From 67fc1074deacad78a8bb2189e3820af194617841 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 8 Jan 2013 13:14:26 +0100 Subject: [PATCH] Change constructor to take ParameterGroup argument. Also make tof to limit against >= 0.0, for case when upstream cell values go below zero. Disabled some debug output. --- examples/compute_tof.cpp | 31 +++++++++--- examples/compute_tof_from_files.cpp | 10 ++-- .../TransportModelTracerTofDiscGal.cpp | 47 +++++++++++++++---- .../TransportModelTracerTofDiscGal.hpp | 4 +- 4 files changed, 67 insertions(+), 25 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index e308782b..1dfc6ec5 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -169,16 +169,19 @@ main(int argc, char** argv) // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); int dg_degree = -1; - bool use_cvi = false; - bool use_limiter = false; bool use_multidim_upwind = false; + // Need to initialize dg solver here, since it uses parameters now. + boost::scoped_ptr dg_solver; if (use_dg) { dg_degree = param.getDefault("dg_degree", 0); - use_cvi = param.getDefault("use_cvi", false); - use_limiter = param.getDefault("use_limiter", false); + dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(*grid->c_grid(), param)); } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } + bool compute_tracer = param.getDefault("compute_tracer", false); + if (use_dg && compute_tracer) { + THROW("DG for tracer not yet implemented."); + } // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -232,12 +235,16 @@ main(int argc, char** argv) // Solve time-of-flight. transport_timer.start(); std::vector tof; + std::vector tracer; if (use_dg) { - Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi, use_limiter); - tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof); + dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof); } else { Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind); - tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); + if (compute_tracer) { + tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tof, tracer); + } else { + tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); + } } transport_timer.stop(); double tt = transport_timer.secsSinceStart(); @@ -249,7 +256,17 @@ main(int argc, char** argv) if (output) { std::string tof_filename = output_dir + "/tof.txt"; std::ofstream tof_stream(tof_filename.c_str()); + tof_stream.precision(16); std::copy(tof.begin(), tof.end(), std::ostream_iterator(tof_stream, "\n")); + if (compute_tracer) { + std::string tracer_filename = output_dir + "/tracer.txt"; + std::ofstream tracer_stream(tracer_filename.c_str()); + tracer_stream.precision(16); + const int nt = tracer.size()/num_cells; + for (int i = 0; i < nt*num_cells; ++i) { + tracer_stream << tracer[i] << (((i + 1) % nt == 0) ? '\n' : ' '); + } + } } std::cout << "\n\n================ End of simulation ===============\n" diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index 4b9e20ce..bf5c077d 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -123,13 +123,12 @@ main(int argc, char** argv) // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); int dg_degree = -1; - bool use_cvi = false; - bool use_limiter = false; bool use_multidim_upwind = false; + // Need to initialize dg solver here, since it uses parameters now. + boost::scoped_ptr dg_solver; if (use_dg) { dg_degree = param.getDefault("dg_degree", 0); - use_cvi = param.getDefault("use_cvi", false); - use_limiter = param.getDefault("use_limiter", false); + dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(grid, param)); } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } @@ -164,8 +163,7 @@ main(int argc, char** argv) std::vector tof; std::vector tracer; if (use_dg) { - Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi, use_limiter); - tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof); + dg_solver->solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof); } else { Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind); if (compute_tracer) { diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index de73f47b..02403ba1 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -128,27 +129,50 @@ namespace Opm /// \param[in] use_cvi If true, use corner point velocity interpolation. /// Otherwise, use the basic constant interpolation. TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, - const bool use_cvi, - const bool use_limiter) + const parameter::ParameterGroup& param) : grid_(grid), - use_cvi_(use_cvi), - use_limiter_(use_limiter), + use_cvi_(false), + use_limiter_(false), limiter_relative_flux_threshold_(1e-3), limiter_method_(MinUpwindFace), limiter_usage_(DuringComputations), coord_(grid.dimensions), velocity_(grid.dimensions) { + use_cvi_ = param.getDefault("use_cvi", use_cvi_); + use_limiter_ = param.getDefault("use_limiter", use_limiter_); + if (use_limiter_) { + limiter_relative_flux_threshold_ = param.getDefault("limiter_relative_flux_threshold", + limiter_relative_flux_threshold_); + const std::string limiter_method_str = param.getDefault("limiter_method", "MinUpwindFace"); + if (limiter_method_str == "MinUpwindFace") { + limiter_method_ = MinUpwindFace; + } else if (limiter_method_str == "MinUpwindAverage") { + limiter_method_ = MinUpwindAverage; + } else { + THROW("Unknown limiter method: " << limiter_method_str); + } + const std::string limiter_usage_str = param.getDefault("limiter_usage", "DuringComputations"); + if (limiter_usage_str == "DuringComputations") { + limiter_usage_ = DuringComputations; + } else if (limiter_usage_str == "AsPostProcess") { + limiter_usage_ = AsPostProcess; + } else if (limiter_usage_str == "AsSimultaneousPostProcess") { + limiter_usage_ = AsSimultaneousPostProcess; + } else { + THROW("Unknown limiter usage spec: " << limiter_usage_str); + } + } // A note about the use_cvi_ member variable: // In principle, we should not need it, since the choice of velocity // interpolation is made below, but we may need to use higher order // quadrature to exploit CVI, so we store the choice. // An alternative would be to add a virtual method isConstant() to // the VelocityInterpolationInterface. - if (use_cvi) { - velocity_interpolation_.reset(new VelocityInterpolationECVI(grid)); + if (use_cvi_) { + velocity_interpolation_.reset(new VelocityInterpolationECVI(grid_)); } else { - velocity_interpolation_.reset(new VelocityInterpolationConstant(grid)); + velocity_interpolation_.reset(new VelocityInterpolationConstant(grid_)); } } @@ -456,7 +480,7 @@ namespace Opm // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. const double total_flux = std::max(-upstream_flux, downstream_flux); - // Find minimum tof on upstream faces. + // Find minimum tof on upstream faces and for this cell. const int dim = grid_.dimensions; const int num_basis = DGBasis::numBasisFunc(dim, degree_); double min_upstream_tof = 1e100; @@ -506,6 +530,9 @@ namespace Opm min_upstream_tof = 0.0; min_here_tof = 0.0; } + if (min_upstream_tof < 0.0) { + min_upstream_tof = 0.0; + } const double tof_c = tof_coeff_[num_basis*cell]; double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof); if (tof_c < min_upstream_tof) { @@ -518,12 +545,12 @@ namespace Opm // Actually do the limiting (if applicable). if (limiter < 1.0) { - std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl; + // std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl; for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) { tof[i] *= limiter; } } else { - std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; + // std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; } } diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 24ad712d..f37cd453 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -33,6 +33,7 @@ namespace Opm class IncompPropertiesInterface; class VelocityInterpolationInterface; + namespace parameter { class ParameterGroup; } /// Implements a discontinuous Galerkin solver for /// (single-phase) time-of-flight using reordering. @@ -51,8 +52,7 @@ namespace Opm /// \param[in] use_cvi If true, use corner point velocity interpolation. /// Otherwise, use the basic constant interpolation. TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, - const bool use_cvi, - const bool use_limiter = false); + const parameter::ParameterGroup& param); /// Solve for time-of-flight.