Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

The mixed approach for the reference solution time integration #50

Closed
wants to merge 2 commits into from
Closed

The mixed approach for the reference solution time integration #50

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6666525
Implemented the HO time-integrated LO solution with the mixed approac…
najlkin Dec 9, 2024
b8d0c18
Added a switch for the HO time integration of the LO solution.
najlkin Dec 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 76 additions & 22 deletions remhos.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool> &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)
{
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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)
{
Expand Down Expand Up @@ -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); }
Expand All @@ -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<Vector&>(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<bool> &mask) const
{
const int NE = Mbf.FESpace()->GetNE();
Expand Down Expand Up @@ -1478,26 +1539,25 @@ void AdvectionOperator::ComputeMask(const Vector &x, Array<bool> &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<Vector&>(X), block_offsets);
const BlockVector block_Y_lo(const_cast<Vector&>(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);
Expand All @@ -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();
Expand Down
40 changes: 38 additions & 2 deletions remhos_solvers.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
{

Expand Down Expand Up @@ -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)
Expand All @@ -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.);
Expand Down Expand Up @@ -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);
Expand All @@ -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.);
Expand Down
26 changes: 25 additions & 1 deletion remhos_solvers.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool> &mask) const
{ MFEM_ABORT("Mask computation not implemented!"); }
Expand All @@ -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
Expand Down Expand Up @@ -97,6 +120,7 @@ class RKIDPSolver : public IDPODESolver
const real_t *a, *b, *c;
real_t *d;
Vector *dxs;
Vector dx_lo;
Array<bool> mask;

// This function constructs coefficients that transform eq. (2.16) from
Expand Down
Loading