From 59441f188a09f329449cd5b7470808d9978311b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 4 Jan 2013 19:59:05 +0100 Subject: [PATCH 1/8] Bugfix: take minimum of allowed slopes, not maximum. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 0dad656e..26076e52 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -404,7 +404,7 @@ namespace Opm const int dim = grid_.dimensions; const int num_basis = DGBasis::numBasisFunc(dim, degree_); - double max_slope_mult = 0.0; + double max_slope_mult = 1e100; int num_upstream_faces = 0; // For inflow faces, ensure that cell tof does not dip below // the minimum value from upstream (for all faces). @@ -453,7 +453,7 @@ namespace Opm break; } const double face_mult = (tof_c - min_upstream)/(tof_c - min_here); - max_slope_mult = std::max(max_slope_mult, face_mult); + max_slope_mult = std::min(max_slope_mult, face_mult); } ASSERT(max_slope_mult >= 0.0); From 7e49c1a475fef96283d1ee54ef6acbc36f7fa07f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 7 Jan 2013 15:10:43 +0100 Subject: [PATCH 2/8] Bugfix, again: it was correct to take the maximum, reinstating. There is a different problem that needs fixing, however: Flux inaccuracies (for example on boundaries) may tag some face as inflow face for a cell, even if it should have been no-flow. This may let the cell avoid limiting, even though it should have been limited according to the proper inflow faces. --- .../transport/reorder/TransportModelTracerTofDiscGal.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 26076e52..70b9b81a 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -404,7 +404,8 @@ namespace Opm const int dim = grid_.dimensions; const int num_basis = DGBasis::numBasisFunc(dim, degree_); - double max_slope_mult = 1e100; + // double max_slope_mult = 1e100; + double max_slope_mult = 0.0; int num_upstream_faces = 0; // For inflow faces, ensure that cell tof does not dip below // the minimum value from upstream (for all faces). @@ -453,7 +454,8 @@ namespace Opm break; } const double face_mult = (tof_c - min_upstream)/(tof_c - min_here); - max_slope_mult = std::min(max_slope_mult, face_mult); + // max_slope_mult = std::min(max_slope_mult, face_mult); + max_slope_mult = std::max(max_slope_mult, face_mult); } ASSERT(max_slope_mult >= 0.0); From cf9fdbbd76954d2424175d0594af6b57c381c938 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 7 Jan 2013 15:48:47 +0100 Subject: [PATCH 3/8] Simplify and correct implementation of limiter. Now we check all corners' tof values for the cell under consideration, not just the inflow face corners'. --- .../TransportModelTracerTofDiscGal.cpp | 68 ++++++++++--------- 1 file changed, 35 insertions(+), 33 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 70b9b81a..2ec842f6 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -404,11 +404,10 @@ namespace Opm const int dim = grid_.dimensions; const int num_basis = DGBasis::numBasisFunc(dim, degree_); - // double max_slope_mult = 1e100; - double max_slope_mult = 0.0; + double min_upstream_tof = 1e100; + double min_here_tof = 1e100; int num_upstream_faces = 0; - // For inflow faces, ensure that cell tof does not dip below - // the minimum value from upstream (for all faces). + // Find minimum tof on upstream faces. for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { const int face = grid_.cell_faces[hface]; double flux = 0.0; @@ -420,48 +419,51 @@ namespace Opm flux = -darcyflux_[face]; upstream_cell = grid_.face_cells[2*face]; } - if (flux >= 0.0) { - // This is a downstream face. - continue; + const bool upstream = (flux < 0.0); + if (upstream) { + ++num_upstream_faces; } - ++num_upstream_faces; + bool interior = (upstream_cell >= 0); - // Evaluate the solution in all corners, and find the appropriate limiter. - bool upstream = (upstream_cell >= 0 && flux < 0.0); - double min_upstream = upstream ? 1e100 : 0.0; - double min_here = 1e100; + // Evaluate the solution in all corners. for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; DGBasis::eval(grid_, cell, degree_, nc, &basis_[0]); const double tof_here = std::inner_product(basis_.begin(), basis_.end(), tof_coeff_ + num_basis*cell, 0.0); - min_here = std::min(min_here, tof_here); + min_here_tof = std::min(min_here_tof, tof_here); if (upstream) { - DGBasis::eval(grid_, upstream_cell, degree_, nc, &basis_nb_[0]); - const double tof_upstream - = std::inner_product(basis_nb_.begin(), basis_nb_.end(), - tof_coeff_ + num_basis*upstream_cell, 0.0); - min_upstream = std::min(min_upstream, tof_upstream); + if (interior) { + DGBasis::eval(grid_, upstream_cell, degree_, nc, &basis_nb_[0]); + const double tof_upstream + = std::inner_product(basis_nb_.begin(), basis_nb_.end(), + tof_coeff_ + num_basis*upstream_cell, 0.0); + min_upstream_tof = std::min(min_upstream_tof, tof_upstream); + } else { + // Allow tof down to 0 on inflow boundaries. + min_upstream_tof = std::min(min_upstream_tof, 0.0); + } } } - // Compute maximum slope multiplier. - const double tof_c = tof_coeff_[num_basis*cell]; - if (tof_c < min_upstream) { - // Handle by setting a flat solution. - std::cout << "Trouble in cell " << cell << std::endl; - max_slope_mult = 0.0; - tof_coeff_[num_basis*cell] = min_upstream; - break; - } - const double face_mult = (tof_c - min_upstream)/(tof_c - min_here); - // max_slope_mult = std::min(max_slope_mult, face_mult); - max_slope_mult = std::max(max_slope_mult, face_mult); } - ASSERT(max_slope_mult >= 0.0); + + // Compute slope multiplier (limiter). + if (num_upstream_faces == 0) { + min_upstream_tof = 0.0; + min_here_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) { + // Handle by setting a flat solution. + std::cout << "Trouble in cell " << cell << std::endl; + limiter = 0.0; + tof_coeff_[num_basis*cell] = min_upstream_tof; + } + ASSERT(limiter >= 0.0); // Actually do the limiting (if applicable). - const double limiter = max_slope_mult; - if (num_upstream_faces > 0 && limiter < 1.0) { + if (limiter < 1.0) { 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_coeff_[i] *= limiter; From 86eb45bd88e3c256b811cf1fa4089bb563b4b95f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 8 Jan 2013 11:04:07 +0100 Subject: [PATCH 4/8] Add methods sequence() and components(). --- .../reorder/TransportModelInterface.cpp | 27 ++++++++++++++----- .../reorder/TransportModelInterface.hpp | 10 ++++++- 2 files changed, 30 insertions(+), 7 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp index 8a90ae34..0a542b88 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/opm/core/transport/reorder/TransportModelInterface.cpp @@ -29,15 +29,18 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux) { // Compute reordered sequence of single-cell problems - std::vector sequence(grid.number_of_cells); - std::vector components(grid.number_of_cells + 1); + sequence_.resize(grid.number_of_cells); + components_.resize(grid.number_of_cells + 1); int ncomponents; time::StopWatch clock; clock.start(); - compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents); + compute_sequence(&grid, darcyflux, &sequence_[0], &components_[0], &ncomponents); clock.stop(); std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl; + // Make vector's size match actual used data. + components_.resize(ncomponents + 1); + // Invoke appropriate solve method for each interdependent component. for (int comp = 0; comp < ncomponents; ++comp) { #if 0 @@ -50,11 +53,23 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g } #endif #endif - const int comp_size = components[comp + 1] - components[comp]; + const int comp_size = components_[comp + 1] - components_[comp]; if (comp_size == 1) { - solveSingleCell(sequence[components[comp]]); + solveSingleCell(sequence_[components_[comp]]); } else { - solveMultiCell(comp_size, &sequence[components[comp]]); + solveMultiCell(comp_size, &sequence_[components_[comp]]); } } } + + +const std::vector& Opm::TransportModelInterface::sequence() const +{ + return sequence_; +} + + +const std::vector& Opm::TransportModelInterface::components() const +{ + return components_; +} diff --git a/opm/core/transport/reorder/TransportModelInterface.hpp b/opm/core/transport/reorder/TransportModelInterface.hpp index 1a93ff36..763ffd2e 100644 --- a/opm/core/transport/reorder/TransportModelInterface.hpp +++ b/opm/core/transport/reorder/TransportModelInterface.hpp @@ -20,6 +20,8 @@ #ifndef OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED #define OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED +#include + struct UnstructuredGrid; namespace Opm @@ -31,7 +33,8 @@ namespace Opm /// method that will have an interface geared to the model's /// needs. (The solve() method is therefore not virtual in this /// class.) The reorderAndTransport() method is provided as an aid - /// to implementing solve() in subclasses. + /// to implementing solve() in subclasses, together with the + /// sequence() and components() methods for accessing the ordering. class TransportModelInterface { public: @@ -41,6 +44,11 @@ namespace Opm virtual void solveMultiCell(const int num_cells, const int* cells) = 0; protected: void reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux); + const std::vector& sequence() const; + const std::vector& components() const; + private: + std::vector sequence_; + std::vector components_; }; From d6596a87afc15c00fa6d22b0208780e06ec97f46 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 8 Jan 2013 11:04:43 +0100 Subject: [PATCH 5/8] Improvements to DG1 flux limiter. Skeleton in place for increased flexibility in methods and usage. (So far behaviour choices are hardcoded, though.) Added relative flux thresholding to existing limiter to avoid flux noise strongly affecting solution. For example no-flow boundaries could be treated as inflow boundaries and make the minimum upwind face limiter meaningless. --- .../TransportModelTracerTofDiscGal.cpp | 104 ++++++++++++++++-- .../TransportModelTracerTofDiscGal.hpp | 14 ++- 2 files changed, 109 insertions(+), 9 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 2ec842f6..de73f47b 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -133,6 +133,9 @@ namespace Opm : grid_(grid), use_cvi_(use_cvi), use_limiter_(use_limiter), + limiter_relative_flux_threshold_(1e-3), + limiter_method_(MinUpwindFace), + limiter_usage_(DuringComputations), coord_(grid.dimensions), velocity_(grid.dimensions) { @@ -194,6 +197,19 @@ namespace Opm grad_basis_.resize(num_basis*grid_.dimensions); velocity_interpolation_->setupFluxes(darcyflux); reorderAndTransport(grid_, darcyflux); + switch (limiter_usage_) { + case AsPostProcess: + applyLimiterAsPostProcess(); + break; + case AsSimultaneousPostProcess: + applyLimiterAsSimultaneousPostProcess(); + break; + case DuringComputations: + // Do nothing. + break; + default: + THROW("Unknown limiter usage choice: " << limiter_usage_); + } } @@ -370,8 +386,8 @@ namespace Opm std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); // Apply limiter. - if (degree_ > 0 && use_limiter_) { - useLimiter(cell); + if (degree_ > 0 && use_limiter_ && limiter_usage_ == DuringComputations) { + applyLimiter(cell, tof_coeff_); } } @@ -389,7 +405,22 @@ namespace Opm - void TransportModelTracerTofDiscGal::useLimiter(const int cell) + void TransportModelTracerTofDiscGal::applyLimiter(const int cell, double* tof) + { + switch (limiter_method_) { + case MinUpwindFace: + applyMinUpwindFaceLimiter(cell, tof); + break; + case MinUpwindAverage: + default: + THROW("Limiter type not implemented: " << limiter_method_); + } + } + + + + + void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof) { if (degree_ != 1) { THROW("This limiter only makes sense for our DG1 implementation."); @@ -399,15 +430,38 @@ namespace Opm // 1. Let M be the minimum TOF value on the upstream faces, // evaluated in the upstream cells. Then the value at all // points in this cell shall be at least M. + // Upstream faces whose flux does not exceed the relative + // flux threshold are not considered for this minimum. // 2. The TOF shall not be below zero in any point. + // Find total upstream/downstream fluxes. + double upstream_flux = 0.0; + double downstream_flux = 0.0; + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + } else { + flux = -darcyflux_[face]; + } + if (flux < 0.0) { + upstream_flux += flux; + } else { + downstream_flux += flux; + } + } + // In the presence of sources, significant fluxes may be missing from the computed fluxes, + // setting the total flux to the (positive) maximum avoids this: since source is either + // 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. const int dim = grid_.dimensions; const int num_basis = DGBasis::numBasisFunc(dim, degree_); - double min_upstream_tof = 1e100; double min_here_tof = 1e100; int num_upstream_faces = 0; - // Find minimum tof on upstream faces. for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { const int face = grid_.cell_faces[hface]; double flux = 0.0; @@ -419,7 +473,7 @@ namespace Opm flux = -darcyflux_[face]; upstream_cell = grid_.face_cells[2*face]; } - const bool upstream = (flux < 0.0); + const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_); if (upstream) { ++num_upstream_faces; } @@ -458,7 +512,7 @@ namespace Opm // Handle by setting a flat solution. std::cout << "Trouble in cell " << cell << std::endl; limiter = 0.0; - tof_coeff_[num_basis*cell] = min_upstream_tof; + tof[num_basis*cell] = min_upstream_tof; } ASSERT(limiter >= 0.0); @@ -466,7 +520,7 @@ namespace Opm if (limiter < 1.0) { 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_coeff_[i] *= limiter; + tof[i] *= limiter; } } else { std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; @@ -476,4 +530,38 @@ namespace Opm + void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess() + { + // Apply the limiter sequentially to all cells. + // This means that a cell's limiting behaviour may be affected by + // any limiting applied to its upstream cells. + const std::vector& seq = TransportModelInterface::sequence(); + const int nc = seq.size(); + ASSERT(nc == grid_.number_of_cells); + for (int i = 0; i < nc; ++i) { + const int cell = seq[i]; + applyLimiter(cell, tof_coeff_); + } + } + + + + + void TransportModelTracerTofDiscGal::applyLimiterAsSimultaneousPostProcess() + { + // Apply the limiter simultaneously to all cells. + // This means that each cell is limited independently from all other cells, + // we write the resulting dofs to a new array instead of writing to tof_coeff_. + // Afterwards we copy the results back to tof_coeff_. + const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_); + std::vector tof_coeffs_new(tof_coeff_, tof_coeff_ + num_basis*grid_.number_of_cells); + for (int c = 0; c < grid_.number_of_cells; ++c) { + applyLimiter(c, &tof_coeffs_new[0]); + } + std::copy(tof_coeffs_new.begin(), tof_coeffs_new.end(), tof_coeff_); + } + + + + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 3188b371..24ad712d 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -88,6 +88,11 @@ namespace Opm boost::shared_ptr velocity_interpolation_; bool use_cvi_; bool use_limiter_; + double limiter_relative_flux_threshold_; + enum LimiterMethod { MinUpwindFace, MinUpwindAverage }; + LimiterMethod limiter_method_; + enum LimiterUsage { DuringComputations, AsPostProcess, AsSimultaneousPostProcess }; + LimiterUsage limiter_usage_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell @@ -105,7 +110,14 @@ namespace Opm std::vector velocity_; // Private methods - void useLimiter(const int cell); + + // Apply some limiter, writing to array tof + // (will read data from tof_coeff_, it is ok to call + // with tof_coeff as tof argument. + void applyLimiter(const int cell, double* tof); + void applyMinUpwindFaceLimiter(const int cell, double* tof); + void applyLimiterAsPostProcess(); + void applyLimiterAsSimultaneousPostProcess(); }; } // namespace Opm 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 6/8] 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. From 9df613913e82d7a2dd442a7d482bdd1e45a07cb2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 8 Jan 2013 13:34:50 +0100 Subject: [PATCH 7/8] Documented new available limiting options in constructor doc. --- .../reorder/TransportModelTracerTofDiscGal.cpp | 17 +++++++++++++++-- .../reorder/TransportModelTracerTofDiscGal.hpp | 17 +++++++++++++++-- 2 files changed, 30 insertions(+), 4 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 02403ba1..82f96585 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -126,8 +126,21 @@ namespace Opm /// Construct solver. /// \param[in] grid A 2d or 3d grid. - /// \param[in] use_cvi If true, use corner point velocity interpolation. - /// Otherwise, use the basic constant interpolation. + /// \param[in] param Parameters for the solver. + /// The following parameters are accepted (defaults): + /// use_cvi (false) Use ECVI velocity interpolation. + /// use_limiter (false) Use a slope limiter. If true, the next three parameters are used. + /// limiter_relative_flux_threshold (1e-3) Ignore upstream fluxes below this threshold, relative to total cell flux. + /// limiter_method ("MinUpwindFace") Limiter method used. Accepted methods are: + /// MinUpwindFace Limit cell tof to >= inflow face tofs. + /// limiter_usage ("DuringComputations") Usage pattern for limiter. Accepted choices are: + /// DuringComputations Apply limiter to cells as they are computed, + /// so downstream cells' solutions may be affected + /// by limiting in upstream cells. + /// AsPostProcess Apply in dependency order, but only after + /// computing (unlimited) solution. + /// AsSimultaneousPostProcess Apply to each cell independently, using un- + /// limited solution in neighbouring cells. TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, const parameter::ParameterGroup& param) : grid_(grid), diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index f37cd453..a92c8fd4 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -49,8 +49,21 @@ namespace Opm public: /// Construct solver. /// \param[in] grid A 2d or 3d grid. - /// \param[in] use_cvi If true, use corner point velocity interpolation. - /// Otherwise, use the basic constant interpolation. + /// \param[in] param Parameters for the solver. + /// The following parameters are accepted (defaults): + /// use_cvi (false) Use ECVI velocity interpolation. + /// use_limiter (false) Use a slope limiter. If true, the next three parameters are used. + /// limiter_relative_flux_threshold (1e-3) Ignore upstream fluxes below this threshold, relative to total cell flux. + /// limiter_method ("MinUpwindFace") Limiter method used. Accepted methods are: + /// MinUpwindFace Limit cell tof to >= inflow face tofs. + /// limiter_usage ("DuringComputations") Usage pattern for limiter. Accepted choices are: + /// DuringComputations Apply limiter to cells as they are computed, + /// so downstream cells' solutions may be affected + /// by limiting in upstream cells. + /// AsPostProcess Apply in dependency order, but only after + /// computing (unlimited) solution. + /// AsSimultaneousPostProcess Apply to each cell independently, using un- + /// limited solution in neighbouring cells. TransportModelTracerTofDiscGal(const UnstructuredGrid& grid, const parameter::ParameterGroup& param); From 33e99520749e5a6570504987fb184e45309843bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 8 Jan 2013 16:00:17 +0100 Subject: [PATCH 8/8] Added MinUpwindAverage limiter method, made it default. --- .../TransportModelTracerTofDiscGal.cpp | 114 +++++++++++++++++- .../TransportModelTracerTofDiscGal.hpp | 1 + 2 files changed, 113 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 82f96585..ae6ab28d 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -147,7 +147,7 @@ namespace Opm use_cvi_(false), use_limiter_(false), limiter_relative_flux_threshold_(1e-3), - limiter_method_(MinUpwindFace), + limiter_method_(MinUpwindAverage), limiter_usage_(DuringComputations), coord_(grid.dimensions), velocity_(grid.dimensions) @@ -157,7 +157,7 @@ namespace Opm 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"); + const std::string limiter_method_str = param.getDefault("limiter_method", "MinUpwindAverage"); if (limiter_method_str == "MinUpwindFace") { limiter_method_ = MinUpwindFace; } else if (limiter_method_str == "MinUpwindAverage") { @@ -449,6 +449,8 @@ namespace Opm applyMinUpwindFaceLimiter(cell, tof); break; case MinUpwindAverage: + applyMinUpwindAverageLimiter(cell, tof); + break; default: THROW("Limiter type not implemented: " << limiter_method_); } @@ -570,6 +572,114 @@ namespace Opm + void TransportModelTracerTofDiscGal::applyMinUpwindAverageLimiter(const int cell, double* tof) + { + if (degree_ != 1) { + THROW("This limiter only makes sense for our DG1 implementation."); + } + + // Limiter principles: + // 1. Let M be the average TOF value of the upstream cells. + /// Then the value at all points in this cell shall be at least M. + // Upstream faces whose flux does not exceed the relative + // flux threshold are not considered for this minimum. + // 2. The TOF shall not be below zero in any point. + + // Find total upstream/downstream fluxes. + double upstream_flux = 0.0; + double downstream_flux = 0.0; + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + } else { + flux = -darcyflux_[face]; + } + if (flux < 0.0) { + upstream_flux += flux; + } else { + downstream_flux += flux; + } + } + // In the presence of sources, significant fluxes may be missing from the computed fluxes, + // setting the total flux to the (positive) maximum avoids this: since source is either + // 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 and for this cell. + const int dim = grid_.dimensions; + const int num_basis = DGBasis::numBasisFunc(dim, degree_); + double min_upstream_tof = 1e100; + double min_here_tof = 1e100; + int num_upstream_faces = 0; + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + int upstream_cell = -1; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face+1]; + } else { + flux = -darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face]; + } + const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_); + if (upstream) { + ++num_upstream_faces; + } + bool interior = (upstream_cell >= 0); + + // Evaluate the solution in all corners. + for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { + const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; + DGBasis::eval(grid_, cell, degree_, nc, &basis_[0]); + const double tof_here = std::inner_product(basis_.begin(), basis_.end(), + tof_coeff_ + num_basis*cell, 0.0); + min_here_tof = std::min(min_here_tof, tof_here); + if (upstream) { + if (interior) { + min_upstream_tof = std::min(min_upstream_tof, tof_coeff_[num_basis*upstream_cell]); + } else { + // Allow tof down to 0 on inflow boundaries. + min_upstream_tof = std::min(min_upstream_tof, 0.0); + } + } + } + } + + // Compute slope multiplier (limiter). + if (num_upstream_faces == 0) { + 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) { + // Handle by setting a flat solution. + std::cout << "Trouble in cell " << cell << std::endl; + limiter = 0.0; + tof[num_basis*cell] = min_upstream_tof; + } + ASSERT(limiter >= 0.0); + + // Actually do the limiting (if applicable). + if (limiter < 1.0) { + // 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; + } + } + + + + void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess() { // Apply the limiter sequentially to all cells. diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index a92c8fd4..529f75b1 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -129,6 +129,7 @@ namespace Opm // with tof_coeff as tof argument. void applyLimiter(const int cell, double* tof); void applyMinUpwindFaceLimiter(const int cell, double* tof); + void applyMinUpwindAverageLimiter(const int cell, double* tof); void applyLimiterAsPostProcess(); void applyLimiterAsSimultaneousPostProcess(); };