diff --git a/remhos.cpp b/remhos.cpp index 00cbea9..f518e33 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -124,11 +124,15 @@ class AdvectionOperator : public LimitedTimeDependentOperator HOSolver *hos, LOSolver *los, FCTSolver *fct, MonolithicSolver *mos); + bool RequiresLOSolution() const override { return (fct_solver != nullptr); } + void MultUnlimited(const Vector &x, Vector &y) const override; + void MultUnlimitedLO(const Vector &x, Vector &y) const override; void ComputeMask(const Vector &x, Array &mask) const override; - void LimitMult(const Vector &x, Vector &y) const override; + void LimitMult(const Vector &x, Vector &y) const override { } + void LimitMult(const Vector &x, const Vector &y_lo, Vector &y) const override; void SetTimeStepControl(TimeStepControl tsc) { @@ -322,8 +326,14 @@ int main(int argc, char *argv[]) case 11: idp_ode_solver = new ForwardEulerIDPSolver(); break; case 12: idp_ode_solver = new RK2IDPSolver(); break; case 13: idp_ode_solver = new RK3IDPSolver(); break; - case 14: idp_ode_solver = new RK4IDPSolver(); break; - case 16: idp_ode_solver = new RK6IDPSolver(); break; + case 14: + if (myid == 0) { MFEM_WARNING("RK4 may violate the bounds."); } + idp_ode_solver = new RK4IDPSolver(); + break; + case 16: + if (myid == 0) { MFEM_WARNING("RK6 may violate the bounds."); } + idp_ode_solver = new RK6IDPSolver(); + break; default: cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; return 3; @@ -1370,11 +1380,12 @@ void AdvectionOperator::MultUnlimited(const Vector &X, Vector &Y) const if (mono_solver) { mono_solver->CalcSolution(u, d_u); } else if (fct_solver) { - MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers."); + MFEM_VERIFY(ho_solver, "FCT requires HO solver."); ho_solver->CalcHOSolution(u, d_u); - // Limiting is deferred to LimitMult() + // Low-order solution is deferred to MultUnlimitedLO() + // while limiting is performed in LimitMult() } else if (lo_solver) { @@ -1407,11 +1418,12 @@ void AdvectionOperator::MultUnlimited(const Vector &X, Vector &Y) const if (mono_solver) { mono_solver->CalcSolution(us, d_us); } else if (fct_solver) { - MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers."); + MFEM_VERIFY(ho_solver, "FCT requires HO solver."); ho_solver->CalcHOSolution(us, d_us); - // Limiting is deferred to LimitMult() + // Low-order solution is deferred to MultUnlimitedLO() + // while limiting is performed in LimitMult() } else if (lo_solver) { lo_solver->CalcLOSolution(us, d_us); } else if (ho_solver) { ho_solver->CalcHOSolution(us, d_us); } @@ -1421,6 +1433,55 @@ void AdvectionOperator::MultUnlimited(const Vector &X, Vector &Y) const } } +void AdvectionOperator::MultUnlimitedLO(const Vector &X, Vector &Y) const +{ + // Needed because X and Y are allocated on the host by the ODESolver. + X.Read(); Y.Write(); + + const BlockVector block_X(const_cast(X), block_offsets); + BlockVector block_Y(Y, block_offsets); + const Vector &u = block_X.GetBlock(0); + Vector &d_u = block_Y.GetBlock(0); + + if (fct_solver) + { + MFEM_VERIFY(lo_solver, "FCT requires LO solver."); + + lo_solver->CalcLOSolution(u, d_u); + + // Limiting is deferred to LimitMult() + } + + d_u.SyncAliasMemory(Y); + + // Remap the product field, if there is a product field. + if (block_offsets.Size() > 2) + { + MFEM_VERIFY(exec_mode == 1, "Products are processed only in remap mode."); + + const Vector &us = block_X.GetBlock(1); + Vector &d_us = block_Y.GetBlock(1); + + if (fct_solver) + { + MFEM_VERIFY(lo_solver, "FCT requires LO solver."); + + if (fct_solver->NeedsLOProductInput()) + { + lo_solver->CalcLOSolution(us, d_us); + } + else + { + d_us = 0.;// dummy + } + + // Limiting is deferred to LimitMult() + } + + d_us.SyncAliasMemory(Y); + } +} + void AdvectionOperator::ComputeMask(const Vector &x, Array &mask) const { const int NE = Mbf.FESpace()->GetNE(); @@ -1478,26 +1539,25 @@ void AdvectionOperator::ComputeMask(const Vector &x, Array &mask) const } } -void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const +void AdvectionOperator::LimitMult(const Vector &X, const Vector &Y_lo, + Vector &Y) const { // Needed because X and Y are allocated on the host by the ODESolver. - X.Read(); Y.Read(); + X.Read(); Y_lo.Read(); Y.Read(); const BlockVector block_X(const_cast(X), block_offsets); + const BlockVector block_Y_lo(const_cast(Y_lo), block_offsets); BlockVector block_Y(Y, block_offsets); const Vector &u = block_X.GetBlock(0); + const Vector &du_LO = block_Y_lo.GetBlock(0); Vector &d_u = block_Y.GetBlock(0); if (fct_solver) { - MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers."); - x_gf = u; x_gf.ExchangeFaceNbrData(); - Vector du_HO(d_u), du_LO(u.Size()); - - lo_solver->CalcLOSolution(u, du_LO); + Vector du_HO(d_u); dofs.ComputeElementsMinMax(u, dofs.xe_min, dofs.xe_max, NULL, NULL); dofs.ComputeBounds(dofs.xe_min, dofs.xe_max, dofs.xi_min, dofs.xi_max); @@ -1520,21 +1580,15 @@ void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const "Automatic time step is not implemented for product remap."); const Vector &us = block_X.GetBlock(1); + const Vector &d_us_LO = block_Y_lo.GetBlock(1); Vector &d_us = block_Y.GetBlock(1); if (fct_solver) { - MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers."); - x_gf = us; x_gf.ExchangeFaceNbrData(); - Vector d_us_HO(d_us), d_us_LO; - if (fct_solver->NeedsLOProductInput()) - { - d_us_LO.SetSize(us.Size()); - lo_solver->CalcLOSolution(us, d_us_LO); - } + Vector d_us_HO(d_us); const int size = Kbf.ParFESpace()->GetVSize(); const int NE = Kbf.ParFESpace()->GetNE(); diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 43a40d1..59f921b 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -17,6 +17,9 @@ #include "remhos_solvers.hpp" #include "general/forall.hpp" +/// Switches on high-order time integration of the low-order solution +#define REMHOS_HO_RK_LO_SOLUTION + namespace mfem { @@ -166,6 +169,10 @@ void RKIDPSolver::Init(LimitedTimeDependentOperator &f_) { dxs[i].SetSize(f->Height()); } + if (f->RequiresLOSolution()) + { + dx_lo.SetSize(f->Height()); + } } void RKIDPSolver::Step(Vector &x, double &t, double &dt) @@ -176,7 +183,15 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) f->SetTime(t); f->SetDt(c[0] * dt); f->MultUnlimited(x, dxs[0]); - f->LimitMult(x, dxs[0]); + if (f->RequiresLOSolution()) + { + f->MultUnlimitedLO(x, dx_lo); + f->LimitMult(x, dx_lo, dxs[0]); + } + else + { + f->LimitMult(x, dxs[0]); + } // Update state const double c_next = (s > 2)?(c[1]):(1.); @@ -208,6 +223,10 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) // Explicit HO step f->SetDt(dct); f->MultUnlimited(x, dxs[i]); + if (f->RequiresLOSolution()) + { + f->MultUnlimitedLO(x, dx_lo); + } // Update mask with the HO update UpdateMask(x, dxs[i], dct, mask); @@ -228,8 +247,25 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) AddMasked(mask, d_i[j], dxs[j], dxs[i]); } +#ifdef REMHOS_HO_RK_LO_SOLUTION + // Combine LO soluion with the previous limited updates to obtain + // a bounds-preserving HO time-integrated reference solution. + if (f->RequiresLOSolution()) + { + if (d_i[i] != 1.) + { + AddMasked(mask, d_i[i]-1., dx_lo, dx_lo); + } + for (int j = 0; j < i; j++) + { + // Use all previous limited updates. + AddMasked(mask, d_i[j], dxs[j], dx_lo); + } + } +#endif // REMHOS_HO_RK_LO_SOLUTION + // Limit the step - f->LimitMult(x, dxs[i]); + f->LimitMult(x, dx_lo, dxs[i]); // Update the state const double c_next = (i < s-2)?(c[i+1]):(1.); diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index d3a983d..7bb43d0 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -43,15 +43,31 @@ class LimitedTimeDependentOperator : public TimeDependentOperator virtual void SetDt(real_t dt_) { dt = dt_; } virtual real_t GetDt() const { return dt; } + /// Checks if the low-order solution is needed + virtual bool RequiresLOSolution() const { return false; } + void Mult(const Vector &u, Vector &k) const override { MultUnlimited(u, k); - LimitMult(u, k); + if (RequiresLOSolution()) + { + Vector k_lo(Height()); + MultUnlimitedLO(u, k_lo); + LimitMult(u, k_lo, k); + } + else + { + LimitMult(u, k); + } } /// Perform the unlimited action of the operator virtual void MultUnlimited(const Vector &u, Vector &k) const = 0; + /// Perform the unlimited low-order action of the operator + virtual void MultUnlimitedLO(const Vector &u, Vector &k) const + { k = 0.; } + /// Compute mask of the state for the update virtual void ComputeMask(const Vector &u, Array &mask) const { MFEM_ABORT("Mask computation not implemented!"); } @@ -60,6 +76,13 @@ class LimitedTimeDependentOperator : public TimeDependentOperator /// Assumes that MultUnlimited(u, k) has been called, which has computed the /// unlimited solution in @a k. virtual void LimitMult(const Vector &u, Vector &k) const = 0; + + /// Limit the action vector @a k with LO solution @a k_lo + /// Assumes that MultUnlimited(u, k) and MultUnlimited(u, k_lo) have been + /// called, which has computed the unlimited solution in @a k and the LO + /// solution in @a k_lo. + virtual void LimitMult(const Vector &u, const Vector &k_lo, Vector &k) const + { LimitMult(u, k); } }; class IDPODESolver : public ODESolver @@ -97,6 +120,7 @@ class RKIDPSolver : public IDPODESolver const real_t *a, *b, *c; real_t *d; Vector *dxs; + Vector dx_lo; Array mask; // This function constructs coefficients that transform eq. (2.16) from