From bcbc0213b69b9ae17ea1b1fa25d07d9d2e8b54a9 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 7 Nov 2024 22:17:05 -0800 Subject: [PATCH 01/39] run lines. --- remhos.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 1bda644..3f8fb1f 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -27,7 +27,11 @@ // as part of the Eulerian phase in Arbitrary-Lagrangian Eulerian (ALE) // simulations. // -// Sample runs: see README.md, section 'Verification of Results'. +// single step: +// mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 14 -rs 3 -dt 0.002 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 1 -vs 20 +// RK2: +// mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 14 -rs 3 -dt 0.002 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 2 -vs 20 +// #include "mfem.hpp" #include @@ -984,7 +988,7 @@ int main(int argc, char *argv[]) if (active_dofs[i] == false) { continue; } double us_min = u(i) * s_min_glob; - if (us(i) < us_min) { us(i) = us_min; } + //if (us(i) < us_min) { us(i) = us_min; } } #ifdef REMHOS_FCT_PRODUCT_DEBUG From 2580131cf74f8db3026aad19c7a3c3ad999d6291 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 14 Nov 2024 16:05:00 -0800 Subject: [PATCH 02/39] Make style. --- remhos.cpp | 6 ++++-- remhos_fct.cpp | 13 +++++++------ remhos_lo.hpp | 8 ++++---- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 3f8fb1f..6f61900 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -48,9 +48,11 @@ using namespace mfem; enum class HOSolverType {None, Neumann, CG, LocalInverse}; enum class FCTSolverType {None, FluxBased, ClipScale, - NonlinearPenalty, FCTProject}; + NonlinearPenalty, FCTProject + }; enum class LOSolverType {None, DiscrUpwind, DiscrUpwindPrec, - ResDist, ResDistSubcell, MassBased}; + ResDist, ResDistSubcell, MassBased + }; enum class MonolithicSolverType {None, ResDistMono, ResDistMonoSubcell}; enum class TimeStepControl {FixedTimeStep, LOBoundsError}; diff --git a/remhos_fct.cpp b/remhos_fct.cpp index be2454b..a8e6d96 100644 --- a/remhos_fct.cpp +++ b/remhos_fct.cpp @@ -745,12 +745,13 @@ void ElementFCTProjection::CalcFCTSolution(const ParGridFunction &u, } // element loop } -void ElementFCTProjection::CalcFCTProduct(const ParGridFunction &us, const Vector &m, - const Vector &d_us_HO, const Vector &d_us_LO, - Vector &s_min, Vector &s_max, - const Vector &u_new, - const Array &active_el, - const Array &active_dofs, Vector &d_us) +void ElementFCTProjection::CalcFCTProduct(const ParGridFunction &us, + const Vector &m, + const Vector &d_us_HO, const Vector &d_us_LO, + Vector &s_min, Vector &s_max, + const Vector &u_new, + const Array &active_el, + const Array &active_dofs, Vector &d_us) { us.HostRead(); s_min.HostReadWrite(); diff --git a/remhos_lo.hpp b/remhos_lo.hpp index 29d59a7..83a3616 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -91,11 +91,11 @@ class MassBasedAvg : public LOSolver Vector &el_mass, Vector &el_vol) const; public: - MassBasedAvg(ParFiniteElementSpace &space, HOSolver &hos, - const GridFunction *mesh_vel) - : LOSolver(space), ho_solver(hos), mesh_v(mesh_vel) { } + MassBasedAvg(ParFiniteElementSpace &space, HOSolver &hos, + const GridFunction *mesh_vel) + : LOSolver(space), ho_solver(hos), mesh_v(mesh_vel) { } - virtual void CalcLOSolution(const Vector &u, Vector &du) const; + virtual void CalcLOSolution(const Vector &u, Vector &du) const; }; //PA based Residual Distribution From 71e31a9d5c40fc525960cee497c0a344d527cceb Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 14 Nov 2024 16:06:10 -0800 Subject: [PATCH 03/39] Turned off hard clipping properly. --- remhos.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/remhos.cpp b/remhos.cpp index 6f61900..88873f7 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -982,6 +982,7 @@ int main(int argc, char *argv[]) // now we have implemented only the minimum global bound. u.HostRead(); us.HostReadWrite(); +#if 0 const int s = u.Size(); Array active_elem, active_dofs; ComputeBoolIndicators(NE, u, active_elem, active_dofs); @@ -990,8 +991,9 @@ int main(int argc, char *argv[]) if (active_dofs[i] == false) { continue; } double us_min = u(i) * s_min_glob; - //if (us(i) < us_min) { us(i) = us_min; } + if (us(i) < us_min) { us(i) = us_min; } } +#endif #ifdef REMHOS_FCT_PRODUCT_DEBUG ComputeMinMaxS(NE, us, u, s_min_glob, s_max_glob); From a47b017cafc110addbeabb4bca2257149e84f242 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 14 Nov 2024 16:08:18 -0800 Subject: [PATCH 04/39] Fixed new min/max calculation. --- remhos.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 88873f7..c62090a 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -1436,10 +1436,11 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const s_bool_el_new, s_bool_dofs_new, d_us); #ifdef REMHOS_FCT_PRODUCT_DEBUG - Vector us_new(size); + Vector us_new(size), s_new(size); add(1.0, us, dt, d_us, us_new); std::cout << " out: "; - ComputeMinMaxS(NE, us_new, u_new, myid); + ComputeRatio(NE, us_new, u_new, s_new, s_bool_el_new, s_bool_dofs_new); + ComputeMinMaxS(s_new, s_bool_dofs_new, myid); #endif } else if (lo_solver) { lo_solver->CalcLOSolution(us, d_us); } From 50171f07774824de94e3b816274f76ee828e81c3 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Thu, 14 Nov 2024 16:51:03 -0800 Subject: [PATCH 05/39] Added RK2 IDP mockup. --- remhos.cpp | 37 ++++++++++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/remhos.cpp b/remhos.cpp index c62090a..6fe2ef9 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -78,6 +78,16 @@ double inflow_function(const Vector &x); // Mesh bounding box Vector bb_min, bb_max; +class RK2IDPSolver : public ODESolver +{ + Vector dx12, dx; + +public: + RK2IDPSolver(); + void Init(TimeDependentOperator &f) override; + void Step(Vector &x, double &t, double &dt) override; +}; + class AdvectionOperator : public TimeDependentOperator { private: @@ -294,7 +304,7 @@ int main(int argc, char *argv[]) switch (ode_solver_type) { case 1: ode_solver = new ForwardEulerSolver; break; - case 2: ode_solver = new RK2Solver(1.0); break; + case 2: ode_solver = new RK2IDPSolver(); break; case 3: ode_solver = new RK3SSPSolver; break; case 4: if (myid == 0) { MFEM_WARNING("RK4 may violate the bounds."); } @@ -1863,3 +1873,28 @@ double inflow_function(const Vector &x) } else { return 0.0; } } + +RK2IDPSolver::RK2IDPSolver() +{ +} + +void RK2IDPSolver::Init(TimeDependentOperator &f) +{ + ODESolver::Init(f); + dx12.SetSize(f.Height()); + dx.SetSize(f.Height()); +} + +void RK2IDPSolver::Step(Vector &x, double &t, double &dt) +{ + f->SetTime(t); + f->Mult(x, dx12); + + x.Add(dt/2., dx12); + f->SetTime(t+dt/2.); + f->Mult(x, dx); + + add(2., dx, -1., dx12, dx); + + x.Add(dt/2., dx); +} From ec14c9d92f0d45e09fd7eb8f82618c9a64811a20 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 10:19:21 -0800 Subject: [PATCH 06/39] Separated IDP solvers from remhos.cpp. --- makefile | 4 ++-- remhos.cpp | 36 +---------------------------------- remhos_solvers.cpp | 47 ++++++++++++++++++++++++++++++++++++++++++++++ remhos_solvers.hpp | 37 ++++++++++++++++++++++++++++++++++++ 4 files changed, 87 insertions(+), 37 deletions(-) create mode 100644 remhos_solvers.cpp create mode 100644 remhos_solvers.hpp diff --git a/makefile b/makefile index 81957b7..bb395cf 100644 --- a/makefile +++ b/makefile @@ -106,11 +106,11 @@ CCC = $(strip $(CXX) $(REMHOS_FLAGS)) Ccc = $(strip $(CC) $(CFLAGS) $(GL_OPTS)) SOURCE_FILES = remhos.cpp remhos_tools.cpp remhos_lo.cpp remhos_ho.cpp \ - remhos_fct.cpp remhos_mono.cpp remhos_sync.cpp + remhos_fct.cpp remhos_mono.cpp remhos_sync.cpp remhos_solvers.cpp OBJECT_FILES1 = $(SOURCE_FILES:.cpp=.o) OBJECT_FILES = $(OBJECT_FILES1:.c=.o) HEADER_FILES = remhos_tools.hpp remhos_lo.hpp remhos_ho.hpp remhos_fct.hpp \ - remhos_mono.hpp remhos_sync.hpp + remhos_mono.hpp remhos_sync.hpp remhos_solvers.hpp # Targets diff --git a/remhos.cpp b/remhos.cpp index 6fe2ef9..90b3f7f 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -42,6 +42,7 @@ #include "remhos_mono.hpp" #include "remhos_tools.hpp" #include "remhos_sync.hpp" +#include "remhos_solvers.hpp" using namespace std; using namespace mfem; @@ -78,16 +79,6 @@ double inflow_function(const Vector &x); // Mesh bounding box Vector bb_min, bb_max; -class RK2IDPSolver : public ODESolver -{ - Vector dx12, dx; - -public: - RK2IDPSolver(); - void Init(TimeDependentOperator &f) override; - void Step(Vector &x, double &t, double &dt) override; -}; - class AdvectionOperator : public TimeDependentOperator { private: @@ -1873,28 +1864,3 @@ double inflow_function(const Vector &x) } else { return 0.0; } } - -RK2IDPSolver::RK2IDPSolver() -{ -} - -void RK2IDPSolver::Init(TimeDependentOperator &f) -{ - ODESolver::Init(f); - dx12.SetSize(f.Height()); - dx.SetSize(f.Height()); -} - -void RK2IDPSolver::Step(Vector &x, double &t, double &dt) -{ - f->SetTime(t); - f->Mult(x, dx12); - - x.Add(dt/2., dx12); - f->SetTime(t+dt/2.); - f->Mult(x, dx); - - add(2., dx, -1., dx12, dx); - - x.Add(dt/2., dx); -} diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp new file mode 100644 index 0000000..8a20f12 --- /dev/null +++ b/remhos_solvers.cpp @@ -0,0 +1,47 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at +// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights +// reserved. See files LICENSE and NOTICE for details. +// +// This file is part of CEED, a collection of benchmarks, miniapps, software +// libraries and APIs for efficient high-order finite element and spectral +// element discretizations for exascale applications. For more information and +// source code availability see http://github.com/ceed. +// +// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, +// a collaborative effort of two U.S. Department of Energy organizations (Office +// of Science and the National Nuclear Security Administration) responsible for +// the planning and preparation of a capable exascale ecosystem, including +// software, applications, hardware, advanced system engineering and early +// testbed platforms, in support of the nation's exascale computing imperative. + +#include "remhos_solvers.hpp" + +namespace mfem +{ + +RK2IDPSolver::RK2IDPSolver() +{ +} + +void RK2IDPSolver::Init(TimeDependentOperator &f) +{ + ODESolver::Init(f); + dx12.SetSize(f.Height()); + dx.SetSize(f.Height()); +} + +void RK2IDPSolver::Step(Vector &x, double &t, double &dt) +{ + f->SetTime(t); + f->Mult(x, dx12); + + x.Add(dt/2., dx12); + f->SetTime(t+dt/2.); + f->Mult(x, dx); + + add(2., dx, -1., dx12, dx); + + x.Add(dt/2., dx); +} + +} // namespace mfem diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp new file mode 100644 index 0000000..89041b2 --- /dev/null +++ b/remhos_solvers.hpp @@ -0,0 +1,37 @@ +// Copyright (c) 2024, Lawrence Livermore National Security, LLC. Produced at +// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights +// reserved. See files LICENSE and NOTICE for details. +// +// This file is part of CEED, a collection of benchmarks, miniapps, software +// libraries and APIs for efficient high-order finite element and spectral +// element discretizations for exascale applications. For more information and +// source code availability see http://github.com/ceed. +// +// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, +// a collaborative effort of two U.S. Department of Energy organizations (Office +// of Science and the National Nuclear Security Administration) responsible for +// the planning and preparation of a capable exascale ecosystem, including +// software, applications, hardware, advanced system engineering and early +// testbed platforms, in support of the nation's exascale computing imperative. + +#ifndef MFEM_REMHOS_SOLVERS +#define MFEM_REMHOS_SOLVERS + +#include "mfem.hpp" + +namespace mfem +{ + +class RK2IDPSolver : public ODESolver +{ + Vector dx12, dx; + +public: + RK2IDPSolver(); + void Init(TimeDependentOperator &f) override; + void Step(Vector &x, double &t, double &dt) override; +}; + +} // namespace mfem + +#endif // MFEM_REMHOS_SOLVERS From 3bd9e021b70809f509bc879cb6b073bc26b982c4 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 10:28:48 -0800 Subject: [PATCH 07/39] Fixed time update in RK2. --- remhos_solvers.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 8a20f12..5150749 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -42,6 +42,7 @@ void RK2IDPSolver::Step(Vector &x, double &t, double &dt) add(2., dx, -1., dx12, dx); x.Add(dt/2., dx); + t += dt; } } // namespace mfem From 8836ce993517355bcd4656179f3797b0cbab29c0 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 11:24:31 -0800 Subject: [PATCH 08/39] Added IDPODESolver and LimitedTimeDependentOperator. --- remhos.cpp | 19 +++++++++++------- remhos_solvers.cpp | 4 ++-- remhos_solvers.hpp | 48 ++++++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 60 insertions(+), 11 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 90b3f7f..3b28a5f 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -79,7 +79,7 @@ double inflow_function(const Vector &x); // Mesh bounding box Vector bb_min, bb_max; -class AdvectionOperator : public TimeDependentOperator +class AdvectionOperator : public LimitedTimeDependentOperator { private: BilinearForm &Mbf, &ml; @@ -119,7 +119,12 @@ class AdvectionOperator : public TimeDependentOperator HOSolver *hos, LOSolver *los, FCTSolver *fct, MonolithicSolver *mos); - virtual void Mult(const Vector &x, Vector &y) const; + void Mult(const Vector &x, Vector &y) const override; + + void MultUnlimited(const Vector &x, Vector &y) const override + { Mult(x, y); } + + void LimitMult(const Vector &x, Vector &y) const override { } void SetTimeStepControl(TimeStepControl tsc) { @@ -291,18 +296,18 @@ int main(int argc, char *argv[]) // Define the ODE solver used for time integration. Several explicit // Runge-Kutta methods are available. - ODESolver *ode_solver = NULL; + IDPODESolver *ode_solver = NULL; switch (ode_solver_type) { - case 1: ode_solver = new ForwardEulerSolver; break; + //case 1: ode_solver = new ForwardEulerSolver; break; case 2: ode_solver = new RK2IDPSolver(); break; - case 3: ode_solver = new RK3SSPSolver; break; + /*case 3: ode_solver = new RK3SSPSolver; break; case 4: if (myid == 0) { MFEM_WARNING("RK4 may violate the bounds."); } ode_solver = new RK4Solver; break; case 6: if (myid == 0) { MFEM_WARNING("RK6 may violate the bounds."); } - ode_solver = new RK6Solver; break; + ode_solver = new RK6Solver; break;*/ default: cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; return 3; @@ -1251,7 +1256,7 @@ AdvectionOperator::AdvectionOperator(int size, BilinearForm &Mbf_, LowOrderMethod &_lom, DofInfo &_dofs, HOSolver *hos, LOSolver *los, FCTSolver *fct, MonolithicSolver *mos) : - TimeDependentOperator(size), Mbf(Mbf_), ml(_ml), Kbf(Kbf_), + LimitedTimeDependentOperator(size), Mbf(Mbf_), ml(_ml), Kbf(Kbf_), M_HO(M_HO_), K_HO(K_HO_), lumpedM(_lumpedM), start_mesh_pos(pos.Size()), start_submesh_pos(sub_vel.Size()), diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 5150749..9d441b5 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -23,9 +23,9 @@ RK2IDPSolver::RK2IDPSolver() { } -void RK2IDPSolver::Init(TimeDependentOperator &f) +void RK2IDPSolver::Init(LimitedTimeDependentOperator &f) { - ODESolver::Init(f); + IDPODESolver::Init(f); dx12.SetSize(f.Height()); dx.SetSize(f.Height()); } diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index 89041b2..a3b008e 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -22,13 +22,57 @@ namespace mfem { -class RK2IDPSolver : public ODESolver +class LimitedTimeDependentOperator : public TimeDependentOperator +{ +public: + /** @brief Construct a "square" LimitedTimeDependentOperator (u,t) -> k(u,t), + where u and k have the same dimension @a n. */ + LimitedTimeDependentOperator(int n = 0, real_t t = 0.0) + : TimeDependentOperator(n, t) { } + + /** @brief Construct a LimitedTimeDependentOperator (u,t) -> k(u,t), where + u and k have dimensions @a w and @a h, respectively. */ + LimitedTimeDependentOperator(int h, int w, real_t t = 0.0) + : TimeDependentOperator(h, w, t) { } + + virtual ~LimitedTimeDependentOperator() { } + + void Mult(const Vector &u, Vector &k) const override + { + MultUnlimited(u, k); + LimitMult(u, k); + } + + /// Perform the unlimited action of the operator + virtual void MultUnlimited(const Vector &u, Vector &k) const = 0; + + /// Limit the action vector @a k + virtual void LimitMult(const Vector &u, Vector &k) const = 0; +}; + +class IDPODESolver : public ODESolver +{ +protected: + /// Pointer to the associated LimitedTimeDependentOperator. + LimitedTimeDependentOperator *f; // f(.,t) : R^n --> R^n + + void Init(TimeDependentOperator &f_) override + { MFEM_ABORT("Limited time-dependent operator must be assigned!"); } +public: + IDPODESolver() : ODESolver(), f(NULL) { } + virtual ~IDPODESolver() { } + + virtual void Init(LimitedTimeDependentOperator &f_) + { ODESolver::Init(f_); f = &f_; } +}; + +class RK2IDPSolver : public IDPODESolver { Vector dx12, dx; public: RK2IDPSolver(); - void Init(TimeDependentOperator &f) override; + void Init(LimitedTimeDependentOperator &f) override; void Step(Vector &x, double &t, double &dt) override; }; From afdc7408371d8a05714a06253eb9acb5baa27d4e Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 11:45:25 -0800 Subject: [PATCH 09/39] Cleaned access to blocks. --- remhos.cpp | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 3b28a5f..7c4c248 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -82,6 +82,8 @@ Vector bb_min, bb_max; class AdvectionOperator : public LimitedTimeDependentOperator { private: + const Array &block_offsets; + BilinearForm &Mbf, &ml; ParBilinearForm &Kbf; ParBilinearForm &M_HO, &K_HO; @@ -109,8 +111,8 @@ class AdvectionOperator : public LimitedTimeDependentOperator const Vector &x_min, const Vector &x_max) const; public: - AdvectionOperator(int size, BilinearForm &Mbf_, BilinearForm &_ml, - Vector &_lumpedM, + AdvectionOperator(const Array &offsets, BilinearForm &Mbf_, + BilinearForm &_ml, Vector &_lumpedM, ParBilinearForm &Kbf_, ParBilinearForm &M_HO_, ParBilinearForm &K_HO_, GridFunction &pos, GridFunction *sub_pos, @@ -896,7 +898,7 @@ int main(int argc, char *argv[]) fct_solver = new ElementFCTProjection(pfes, dt); } - AdvectionOperator adv(S.Size(), m, ml, lumpedM, k, M_HO, K_HO, + AdvectionOperator adv(offset, m, ml, lumpedM, k, M_HO, K_HO, x, xsub, v_gf, v_sub_gf, asmbl, lom, dofs, ho_solver, lo_solver, fct_solver, mono_solver); @@ -1246,7 +1248,8 @@ int main(int argc, char *argv[]) return 0; } -AdvectionOperator::AdvectionOperator(int size, BilinearForm &Mbf_, +AdvectionOperator::AdvectionOperator(const Array &offsets, + BilinearForm &Mbf_, BilinearForm &_ml, Vector &_lumpedM, ParBilinearForm &Kbf_, ParBilinearForm &M_HO_, ParBilinearForm &K_HO_, @@ -1256,7 +1259,8 @@ AdvectionOperator::AdvectionOperator(int size, BilinearForm &Mbf_, LowOrderMethod &_lom, DofInfo &_dofs, HOSolver *hos, LOSolver *los, FCTSolver *fct, MonolithicSolver *mos) : - LimitedTimeDependentOperator(size), Mbf(Mbf_), ml(_ml), Kbf(Kbf_), + LimitedTimeDependentOperator(offsets.Last()), block_offsets(offsets), + Mbf(Mbf_), ml(_ml), Kbf(Kbf_), M_HO(M_HO_), K_HO(K_HO_), lumpedM(_lumpedM), start_mesh_pos(pos.Size()), start_submesh_pos(sub_vel.Size()), @@ -1339,10 +1343,10 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const // Needed because X and Y are allocated on the host by the ODESolver. X.Read(); Y.Read(); - Vector u, d_u; - Vector* xptr = const_cast(&X); - u.MakeRef(*xptr, 0, size); - d_u.MakeRef(Y, 0, size); + 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); Vector du_HO(u.Size()), du_LO(u.Size()); x_gf = u; @@ -1385,15 +1389,14 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const d_u.SyncAliasMemory(Y); // Remap the product field, if there is a product field. - if (X.Size() > size) + if (block_offsets.Size() > 2) { MFEM_VERIFY(exec_mode == 1, "Products are processed only in remap mode."); MFEM_VERIFY(dt_control == TimeStepControl::FixedTimeStep, "Automatic time step is not implemented for product remap."); - Vector us, d_us; - us.MakeRef(*xptr, size, size); - d_us.MakeRef(Y, size, size); + const Vector &us = block_X.GetBlock(1); + Vector &d_us = block_Y.GetBlock(1); x_gf = us; x_gf.ExchangeFaceNbrData(); From ac929b4d8fc96654f5175414f60f0afe330f0c8b Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 14:08:22 -0800 Subject: [PATCH 10/39] Moved setting of time step to LimitedTimeDependentOperator. Does not work with adaptive time control atm! --- remhos.cpp | 5 +++-- remhos_solvers.cpp | 2 ++ remhos_solvers.hpp | 6 ++++++ 3 files changed, 11 insertions(+), 2 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 7c4c248..4b42131 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -94,7 +94,6 @@ class AdvectionOperator : public LimitedTimeDependentOperator mutable ParGridFunction x_gf; - double dt; TimeStepControl dt_control; mutable double dt_est; Assembly &asmbl; @@ -137,7 +136,9 @@ class AdvectionOperator : public LimitedTimeDependentOperator } dt_control = tsc; } - void SetDt(double _dt) { dt = _dt; dt_est = dt; } + + void SetDt(double _dt) override { LimitedTimeDependentOperator::SetDt(_dt); dt_est = dt; } + double GetTimeStepEstimate() { return dt_est; } void SetRemapStartPos(const Vector &m_pos, const Vector &sm_pos) diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 9d441b5..8a0c445 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -33,10 +33,12 @@ void RK2IDPSolver::Init(LimitedTimeDependentOperator &f) void RK2IDPSolver::Step(Vector &x, double &t, double &dt) { f->SetTime(t); + f->SetDt(dt/2.); f->Mult(x, dx12); x.Add(dt/2., dx12); f->SetTime(t+dt/2.); + //f->SetDt(dt/2.); f->Mult(x, dx); add(2., dx, -1., dx12, dx); diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index a3b008e..8d5f037 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -24,6 +24,9 @@ namespace mfem class LimitedTimeDependentOperator : public TimeDependentOperator { +protected: + real_t dt; + public: /** @brief Construct a "square" LimitedTimeDependentOperator (u,t) -> k(u,t), where u and k have the same dimension @a n. */ @@ -37,6 +40,9 @@ class LimitedTimeDependentOperator : public TimeDependentOperator virtual ~LimitedTimeDependentOperator() { } + virtual void SetDt(double dt_) { dt = dt_; } + virtual real_t GetDt() const { return dt; } + void Mult(const Vector &u, Vector &k) const override { MultUnlimited(u, k); From d6aef34da801ae9cb6e657bec5c40574fee5fce9 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 14:10:06 -0800 Subject: [PATCH 11/39] Added .gitignore. --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4704bd8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*.o +remhos +*.mesh +*.gf From 938c58c14eb54c22f0f8a91c957020b192b86148 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 14:37:39 -0800 Subject: [PATCH 12/39] Split AdvectionOperator::Mult(). Limiting does not work atm. --- remhos.cpp | 105 +++++++++++++++++++++++++++++++-------------- remhos_solvers.cpp | 6 ++- 2 files changed, 77 insertions(+), 34 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 4b42131..b948b7c 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -120,12 +120,9 @@ class AdvectionOperator : public LimitedTimeDependentOperator HOSolver *hos, LOSolver *los, FCTSolver *fct, MonolithicSolver *mos); - void Mult(const Vector &x, Vector &y) const override; + void MultUnlimited(const Vector &x, Vector &y) const override; - void MultUnlimited(const Vector &x, Vector &y) const override - { Mult(x, y); } - - void LimitMult(const Vector &x, Vector &y) const override { } + void LimitMult(const Vector &x, Vector &y) const override; void SetTimeStepControl(TimeStepControl tsc) { @@ -1271,7 +1268,7 @@ AdvectionOperator::AdvectionOperator(const Array &offsets, asmbl(_asmbl), lom(_lom), dofs(_dofs), ho_solver(hos), lo_solver(los), fct_solver(fct), mono_solver(mos) { } -void AdvectionOperator::Mult(const Vector &X, Vector &Y) const +void AdvectionOperator::MultUnlimited(const Vector &X, Vector &Y) const { if (exec_mode == 1) { @@ -1338,9 +1335,6 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const } } - const int size = Kbf.ParFESpace()->GetVSize(); - const int NE = Kbf.ParFESpace()->GetNE(); - // Needed because X and Y are allocated on the host by the ODESolver. X.Read(); Y.Read(); @@ -1348,28 +1342,15 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const BlockVector block_Y(Y, block_offsets); const Vector &u = block_X.GetBlock(0); Vector &d_u = block_Y.GetBlock(0); - Vector du_HO(u.Size()), du_LO(u.Size()); - - x_gf = u; - x_gf.ExchangeFaceNbrData(); 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."); - lo_solver->CalcLOSolution(u, du_LO); - ho_solver->CalcHOSolution(u, du_HO); - - 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); - fct_solver->CalcFCTSolution(x_gf, lumpedM, du_HO, du_LO, - dofs.xi_min, dofs.xi_max, d_u); + ho_solver->CalcHOSolution(u, d_u); - if (dt_control == TimeStepControl::LOBoundsError) - { - UpdateTimeStepEstimate(u, du_LO, dofs.xi_min, dofs.xi_max); - } + // Limiting is deferred to LimitMult() } else if (lo_solver) { @@ -1399,21 +1380,84 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const const Vector &us = block_X.GetBlock(1); Vector &d_us = block_Y.GetBlock(1); - x_gf = us; - x_gf.ExchangeFaceNbrData(); - 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."); - Vector d_us_HO(us.Size()), d_us_LO; + ho_solver->CalcHOSolution(us, d_us); + + // Limiting is deferred to LimitMult() + } + else if (lo_solver) { lo_solver->CalcLOSolution(us, d_us); } + else if (ho_solver) { ho_solver->CalcHOSolution(us, d_us); } + else { MFEM_ABORT("No solver was chosen."); } + + d_us.SyncAliasMemory(Y); + } +} + +void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const +{ + // Needed because X and Y are allocated on the host by the ODESolver. + X.Read(); Y.Read(); + + 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(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); + + 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); + fct_solver->CalcFCTSolution(x_gf, lumpedM, du_HO, du_LO, + dofs.xi_min, dofs.xi_max, d_u); + + if (dt_control == TimeStepControl::LOBoundsError) + { + UpdateTimeStepEstimate(u, du_LO, dofs.xi_min, dofs.xi_max); + } + } + + 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."); + MFEM_VERIFY(dt_control == TimeStepControl::FixedTimeStep, + "Automatic time step is not implemented for product remap."); + + const Vector &us = block_X.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); } - ho_solver->CalcHOSolution(us, d_us_HO); + + const int size = Kbf.ParFESpace()->GetVSize(); + const int NE = Kbf.ParFESpace()->GetNE(); // Compute the ratio s = us_old / u_old, and old active dofs. Vector s(size); @@ -1453,9 +1497,6 @@ void AdvectionOperator::Mult(const Vector &X, Vector &Y) const ComputeMinMaxS(s_new, s_bool_dofs_new, myid); #endif } - else if (lo_solver) { lo_solver->CalcLOSolution(us, d_us); } - else if (ho_solver) { ho_solver->CalcHOSolution(us, d_us); } - else { MFEM_ABORT("No solver was chosen."); } d_us.SyncAliasMemory(Y); } diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 8a0c445..e02d2d8 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -34,14 +34,16 @@ void RK2IDPSolver::Step(Vector &x, double &t, double &dt) { f->SetTime(t); f->SetDt(dt/2.); - f->Mult(x, dx12); + f->MultUnlimited(x, dx12); + f->LimitMult(x, dx12); x.Add(dt/2., dx12); f->SetTime(t+dt/2.); //f->SetDt(dt/2.); - f->Mult(x, dx); + f->MultUnlimited(x, dx); add(2., dx, -1., dx12, dx); + f->LimitMult(x, dx); x.Add(dt/2., dx); t += dt; From 92c599a90ed68c413b5aa4aa6f8bf54970c71eeb Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 14:46:35 -0800 Subject: [PATCH 13/39] Added forward Euler. --- remhos.cpp | 2 +- remhos_solvers.cpp | 15 ++++++++++++++- remhos_solvers.hpp | 9 ++++++++- 3 files changed, 23 insertions(+), 3 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index b948b7c..703f785 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -299,7 +299,7 @@ int main(int argc, char *argv[]) IDPODESolver *ode_solver = NULL; switch (ode_solver_type) { - //case 1: ode_solver = new ForwardEulerSolver; break; + case 1: ode_solver = new ForwardEulerIDPSolver(); break; case 2: ode_solver = new RK2IDPSolver(); break; /*case 3: ode_solver = new RK3SSPSolver; break; case 4: diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index e02d2d8..3edc2f0 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -19,8 +19,21 @@ namespace mfem { -RK2IDPSolver::RK2IDPSolver() +void ForwardEulerIDPSolver::Init(LimitedTimeDependentOperator &f) { + IDPODESolver::Init(f); + dx.SetSize(f.Height()); +} + +void ForwardEulerIDPSolver::Step(Vector &x, double &t, double &dt) +{ + f->SetTime(t); + f->SetDt(dt); + f->MultUnlimited(x, dx); + f->LimitMult(x, dx); + + x.Add(dt, dx); + t += dt; } void RK2IDPSolver::Init(LimitedTimeDependentOperator &f) diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index 8d5f037..f992f1f 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -72,12 +72,19 @@ class IDPODESolver : public ODESolver { ODESolver::Init(f_); f = &f_; } }; +class ForwardEulerIDPSolver : public IDPODESolver +{ + Vector dx; +public: + void Init(LimitedTimeDependentOperator &f) override; + void Step(Vector &x, double &t, double &dt) override; +}; + class RK2IDPSolver : public IDPODESolver { Vector dx12, dx; public: - RK2IDPSolver(); void Init(LimitedTimeDependentOperator &f) override; void Step(Vector &x, double &t, double &dt) override; }; From df2bb5ab3ccc58abf6a42186f83b8e33c728ff83 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 15:02:28 -0800 Subject: [PATCH 14/39] Fixed time step setting for FCT and LO solvers. --- remhos.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 703f785..ec59e3f 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -134,7 +134,13 @@ class AdvectionOperator : public LimitedTimeDependentOperator dt_control = tsc; } - void SetDt(double _dt) override { LimitedTimeDependentOperator::SetDt(_dt); dt_est = dt; } + void SetDt(double _dt) override + { + LimitedTimeDependentOperator::SetDt(_dt); + dt_est = dt; + if (lo_solver) { lo_solver->UpdateTimeStep(dt); } + if (fct_solver) { fct_solver->UpdateTimeStep(dt); } + } double GetTimeStepEstimate() { return dt_est; } @@ -931,8 +937,6 @@ int main(int argc, char *argv[]) // This also resets the time step estimate when automatic dt is on. adv.SetDt(dt_real); - if (lo_solver) { lo_solver->UpdateTimeStep(dt_real); } - if (fct_solver) { fct_solver->UpdateTimeStep(dt_real); } if (product_sync) { From 84074c014d0c7f677911765d75bec2b4cd1c9ccb Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Fri, 15 Nov 2024 15:12:32 -0800 Subject: [PATCH 15/39] Fixed the debugging output. --- remhos.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index ec59e3f..acc13d7 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -1470,7 +1470,7 @@ void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const #ifdef REMHOS_FCT_PRODUCT_DEBUG const int myid = x_gf.ParFESpace()->GetMyRank(); if (myid == 0) { std::cout << " --- RK stage" << std::endl; } - std::cout << " in: "; + if (myid == 0) { std::cout << " in: "; } ComputeMinMaxS(s, s_bool_dofs, myid); #endif @@ -1496,7 +1496,7 @@ void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const #ifdef REMHOS_FCT_PRODUCT_DEBUG Vector us_new(size), s_new(size); add(1.0, us, dt, d_us, us_new); - std::cout << " out: "; + if (myid == 0) { std::cout << " out: "; } ComputeRatio(NE, us_new, u_new, s_new, s_bool_el_new, s_bool_dofs_new); ComputeMinMaxS(s_new, s_bool_dofs_new, myid); #endif From 0ac5b945db5bb2d059faea63b7d8aeae3faddca8 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 18 Nov 2024 11:22:06 -0800 Subject: [PATCH 16/39] Added masking of inactive DOFs. RK2 works now. --- remhos.cpp | 53 ++++++++++++++++++++++++++++++++++++++++++ remhos_solvers.cpp | 58 +++++++++++++++++++++++++++++++++++++++++++++- remhos_solvers.hpp | 11 +++++++++ 3 files changed, 121 insertions(+), 1 deletion(-) diff --git a/remhos.cpp b/remhos.cpp index acc13d7..0b3c350 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -122,6 +122,8 @@ class AdvectionOperator : public LimitedTimeDependentOperator void MultUnlimited(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 SetTimeStepControl(TimeStepControl tsc) @@ -1401,6 +1403,57 @@ void AdvectionOperator::MultUnlimited(const Vector &X, Vector &Y) const } } +void AdvectionOperator::ComputeMask(const Vector &x, Array &mask) const +{ + const int NE = Mbf.FESpace()->GetNE(); + Array bool_els; + + if (block_offsets.Size() <= 2) + { + // Only product fields must be masked + mask.SetSize(x.Size()); + mask = true; + return; + } + + const BlockVector bx(const_cast(x), block_offsets); + Array bool_dofs; + ComputeBoolIndicators(NE, bx.GetBlock(0), bool_els, bool_dofs); + + // All DOFs must be active for consistency of update + const int ndof_el = bx.GetBlock(0).Size() / NE; + for (int k = 0; k < NE; k++) + { + if (!bool_els[k]) { continue; } + + bool dofs_active = true; + for (int j = 0; j < ndof_el; j++) + { + if (!bool_dofs[j+ndof_el*k]) + { + dofs_active = false; + break; + } + } + + if (!dofs_active) + for (int j = 0; j < ndof_el; j++) + { + bool_dofs[j+ndof_el*k] = false; + } + } + + // Apply the u mask to all product fields + const int ndofs = bool_dofs.Size(); + const int ndim = block_offsets.Size() - 1; + mask.SetSize(ndofs * ndim); + for (int i = 0; i < ndofs; i++) + for (int d = 0; d < ndim; d++) + { + mask[i+ndofs*d] = bool_dofs[i]; + } +} + void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const { // Needed because X and Y are allocated on the host by the ODESolver. diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 3edc2f0..fd45288 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -15,10 +15,64 @@ // testbed platforms, in support of the nation's exascale computing imperative. #include "remhos_solvers.hpp" +#include "general/forall.hpp" namespace mfem { +void IDPODESolver::AddMasked(const Array &mask, real_t a, + const Vector &va, + real_t b, const Vector &vb, Vector &vc) +{ + MFEM_ASSERT(va.Size() == vb.Size() && va.Size() == vc.Size(), + "incompatible Vectors!"); + +#if !defined(MFEM_USE_LEGACY_OPENMP) + const bool use_dev = va.UseDevice() || va.UseDevice() || vc.UseDevice(); + const int N = vc.Size(); + // Note: get read access first, in case c is the same as a/b. + auto ad = va.Read(use_dev); + auto bd = vb.Read(use_dev); + auto cd = vc.Write(use_dev); + auto maskd = mask.Read(use_dev); + mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) + { + cd[i] = (maskd[i])?(a * ad[i] + b * bd[i]):(cd[i]); + }); +#else + const real_t *ap = va.GetData(); + const real_t *bp = vb.GetData(); + real_t *cp = vc.GetData(); + const bool *maskp = mask.GetData(); + const int s = vc.Size(); + #pragma omp parallel for + for (int i = 0; i < s; i++) + { + cp[i] = (maskp[i])?(a * ap[i] + b * bp[i]):(cp[i]); + } +#endif +} + +void IDPODESolver::UpdateMask(const Vector &x, const Vector &dx, real_t dt, + Array &mask) +{ + Array mask_new(mask.Size()); + if (dt != 0.) + { + Vector x_new(x.Size()); + add(x, dt, dx, x_new); + f->ComputeMask(x_new, mask_new); + } + else + { + f->ComputeMask(x, mask_new); + } + for (int i = 0; i < mask.Size(); i++) + { + mask[i] = mask[i] && mask_new[i]; + } +} + void ForwardEulerIDPSolver::Init(LimitedTimeDependentOperator &f) { IDPODESolver::Init(f); @@ -51,11 +105,13 @@ void RK2IDPSolver::Step(Vector &x, double &t, double &dt) f->LimitMult(x, dx12); x.Add(dt/2., dx12); + f->ComputeMask(x, mask); f->SetTime(t+dt/2.); //f->SetDt(dt/2.); f->MultUnlimited(x, dx); - add(2., dx, -1., dx12, dx); + UpdateMask(x, dx, dt/2., mask); + AddMasked(mask, 2., dx, -1., dx12, dx); f->LimitMult(x, dx); x.Add(dt/2., dx); diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index f992f1f..b6337ee 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -52,6 +52,10 @@ class LimitedTimeDependentOperator : public TimeDependentOperator /// Perform the unlimited action of the operator virtual void MultUnlimited(const Vector &u, Vector &k) const = 0; + /// Compute mask of the state for the update + virtual void ComputeMask(const Vector &u, Array &mask) const + { MFEM_ABORT("Mask computation not implemented!"); } + /// Limit the action vector @a k virtual void LimitMult(const Vector &u, Vector &k) const = 0; }; @@ -62,6 +66,12 @@ class IDPODESolver : public ODESolver /// Pointer to the associated LimitedTimeDependentOperator. LimitedTimeDependentOperator *f; // f(.,t) : R^n --> R^n + void AddMasked(const Array &mask, real_t a, const Vector &va, real_t b, + const Vector &vb, Vector &vc); + + void UpdateMask(const Vector &x, const Vector &dx, real_t dt, + Array &mask); + void Init(TimeDependentOperator &f_) override { MFEM_ABORT("Limited time-dependent operator must be assigned!"); } public: @@ -83,6 +93,7 @@ class ForwardEulerIDPSolver : public IDPODESolver class RK2IDPSolver : public IDPODESolver { Vector dx12, dx; + Array mask; public: void Init(LimitedTimeDependentOperator &f) override; From 9fb62bc6cbe4a554b3f46c3244f5e5b1d5286f51 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 18 Nov 2024 14:25:35 -0800 Subject: [PATCH 17/39] Added arbitrary order RK IDP solver and changed RK2 to it. --- remhos_solvers.cpp | 181 ++++++++++++++++++++++++++++++++++----------- remhos_solvers.hpp | 31 ++++++-- 2 files changed, 162 insertions(+), 50 deletions(-) diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index fd45288..0445c25 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -20,41 +20,105 @@ namespace mfem { -void IDPODESolver::AddMasked(const Array &mask, real_t a, - const Vector &va, - real_t b, const Vector &vb, Vector &vc) +void ForwardEulerIDPSolver::Init(LimitedTimeDependentOperator &f) +{ + IDPODESolver::Init(f); + dx.SetSize(f.Height()); +} + +void ForwardEulerIDPSolver::Step(Vector &x, double &t, double &dt) +{ + f->SetTime(t); + f->SetDt(dt); + f->MultUnlimited(x, dx); + f->LimitMult(x, dx); + + x.Add(dt, dx); + t += dt; +} + +void RKIDPSolver::ConstructD() +{ + // Convert high-order to Forward Euler factors + d = new real_t[s*(s+1)/2]; + + const real_t *a_n = a; + const real_t *a_o = a; + int i_o = -1; + real_t c_o = 0.; + + for (int i = 0; i < s; i++) + { + const real_t c_n = (i &mask, real_t b, const Vector &vb, + Vector &va) { - MFEM_ASSERT(va.Size() == vb.Size() && va.Size() == vc.Size(), + MFEM_ASSERT(va.Size() == vb.Size(), "incompatible Vectors!"); #if !defined(MFEM_USE_LEGACY_OPENMP) - const bool use_dev = va.UseDevice() || va.UseDevice() || vc.UseDevice(); - const int N = vc.Size(); + const bool use_dev = va.UseDevice() || va.UseDevice(); + const int N = va.Size(); // Note: get read access first, in case c is the same as a/b. - auto ad = va.Read(use_dev); + auto ad = va.ReadWrite(use_dev); auto bd = vb.Read(use_dev); - auto cd = vc.Write(use_dev); auto maskd = mask.Read(use_dev); mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { - cd[i] = (maskd[i])?(a * ad[i] + b * bd[i]):(cd[i]); + ad[i] += (maskd[i])?(b * bd[i]):(0.); }); #else - const real_t *ap = va.GetData(); + real_t *ap = va.GetData(); const real_t *bp = vb.GetData(); - real_t *cp = vc.GetData(); const bool *maskp = mask.GetData(); - const int s = vc.Size(); + const int s = va.Size(); #pragma omp parallel for for (int i = 0; i < s; i++) { - cp[i] = (maskp[i])?(a * ap[i] + b * bp[i]):(cp[i]); + ap[i] += (maskp[i])?(b * bp[i]):(0.); } #endif } -void IDPODESolver::UpdateMask(const Vector &x, const Vector &dx, real_t dt, - Array &mask) +void RKIDPSolver::UpdateMask(const Vector &x, const Vector &dx, real_t dt, + Array &mask) { Array mask_new(mask.Size()); if (dt != 0.) @@ -73,49 +137,82 @@ void IDPODESolver::UpdateMask(const Vector &x, const Vector &dx, real_t dt, } } -void ForwardEulerIDPSolver::Init(LimitedTimeDependentOperator &f) +RKIDPSolver::RKIDPSolver(int s_, const real_t a_[], const real_t b_[], + const real_t c_[]) + : s(s_), a(a_), b(b_), c(c_) { - IDPODESolver::Init(f); - dx.SetSize(f.Height()); + dxs = new Vector[s]; + ConstructD(); } -void ForwardEulerIDPSolver::Step(Vector &x, double &t, double &dt) +RKIDPSolver::~RKIDPSolver() { - f->SetTime(t); - f->SetDt(dt); - f->MultUnlimited(x, dx); - f->LimitMult(x, dx); - - x.Add(dt, dx); - t += dt; + delete[] dxs; } -void RK2IDPSolver::Init(LimitedTimeDependentOperator &f) +void RKIDPSolver::Init(LimitedTimeDependentOperator &f_) { - IDPODESolver::Init(f); - dx12.SetSize(f.Height()); - dx.SetSize(f.Height()); + IDPODESolver::Init(f_); + for (int i = 0; i < s; i++) + { + dxs[i].SetSize(f->Height()); + } } -void RK2IDPSolver::Step(Vector &x, double &t, double &dt) +void RKIDPSolver::Step(Vector &x, double &t, double &dt) { + // Perform the first step f->SetTime(t); - f->SetDt(dt/2.); - f->MultUnlimited(x, dx12); - f->LimitMult(x, dx12); + f->SetDt(c[0] * dt); + f->MultUnlimited(x, dxs[0]); + f->LimitMult(x, dxs[0]); - x.Add(dt/2., dx12); + x.Add(c[0] * dt, dxs[0]); f->ComputeMask(x, mask); - f->SetTime(t+dt/2.); - //f->SetDt(dt/2.); - f->MultUnlimited(x, dx); + f->SetTime(c[0] * dt); - UpdateMask(x, dx, dt/2., mask); - AddMasked(mask, 2., dx, -1., dx12, dx); - f->LimitMult(x, dx); + // Step through higher stages + + const real_t *d_i = d + 1; + + for (int i = 1; i < s; i++) + { + const real_t c_o = c[i-1]; + const real_t c_n = (iSetDt(dct); + f->MultUnlimited(x, dxs[i]); + + // Update mask with the HO update + UpdateMask(x, dxs[i], dct, mask); + + // Convert the HO update to Forward Euler + if (d_i[i] != 1.) + { + AddMasked(mask, d_i[i]-1., dxs[i], dxs[i]); + } + for (int j = 0; j < i; j++) + { + AddMasked(mask, d_i[j], dxs[j], dxs[i]); + } + + // Limit the step + f->LimitMult(x, dxs[i]); + + // Update the state + f->SetTime(t + c_n * dt); + x.Add(dct, dxs[i]); + d_i += i+1; + } - x.Add(dt/2., dx); t += dt; } +const real_t RK2IDPSolver::a[] = {.5}; +const real_t RK2IDPSolver::b[] = {0., 1.}; +const real_t RK2IDPSolver::c[] = {.5}; + } // namespace mfem diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index b6337ee..3a551a9 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -66,12 +66,6 @@ class IDPODESolver : public ODESolver /// Pointer to the associated LimitedTimeDependentOperator. LimitedTimeDependentOperator *f; // f(.,t) : R^n --> R^n - void AddMasked(const Array &mask, real_t a, const Vector &va, real_t b, - const Vector &vb, Vector &vc); - - void UpdateMask(const Vector &x, const Vector &dx, real_t dt, - Array &mask); - void Init(TimeDependentOperator &f_) override { MFEM_ABORT("Limited time-dependent operator must be assigned!"); } public: @@ -90,16 +84,37 @@ class ForwardEulerIDPSolver : public IDPODESolver void Step(Vector &x, double &t, double &dt) override; }; -class RK2IDPSolver : public IDPODESolver +class RKIDPSolver : public IDPODESolver { - Vector dx12, dx; + const int s; + const real_t *a, *b, *c; + real_t *d; + Vector *dxs; Array mask; + void ConstructD(); + + void AddMasked(const Array &mask, real_t b, + const Vector &vb, Vector &va); + + void UpdateMask(const Vector &x, const Vector &dx, real_t dt, + Array &mask); + public: + RKIDPSolver(int s_, const real_t a_[], const real_t b_[], const real_t c_[]); + ~RKIDPSolver(); + void Init(LimitedTimeDependentOperator &f) override; void Step(Vector &x, double &t, double &dt) override; }; +class RK2IDPSolver : public RKIDPSolver +{ + static const real_t a[], b[], c[]; +public: + RK2IDPSolver() : RKIDPSolver(2, a, b, c) { } +}; + } // namespace mfem #endif // MFEM_REMHOS_SOLVERS From 7340dc7862ddaab4d62d953442b799d9efc7b4fa Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 18 Nov 2024 14:35:01 -0800 Subject: [PATCH 18/39] Added RK3 and equidistant RK4. --- remhos.cpp | 8 +++----- remhos_solvers.cpp | 11 +++++++++++ remhos_solvers.hpp | 14 ++++++++++++++ 3 files changed, 28 insertions(+), 5 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 0b3c350..ffec017 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -309,11 +309,9 @@ int main(int argc, char *argv[]) { case 1: ode_solver = new ForwardEulerIDPSolver(); break; case 2: ode_solver = new RK2IDPSolver(); break; - /*case 3: ode_solver = new RK3SSPSolver; break; - case 4: - if (myid == 0) { MFEM_WARNING("RK4 may violate the bounds."); } - ode_solver = new RK4Solver; break; - case 6: + case 3: ode_solver = new RK3IDPSolver(); break; + case 4: ode_solver = new RK4IDPSolver(); break; + /*case 6: if (myid == 0) { MFEM_WARNING("RK6 may violate the bounds."); } ode_solver = new RK6Solver; break;*/ default: diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 0445c25..7f50cea 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -211,8 +211,19 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) t += dt; } +//2nd order const real_t RK2IDPSolver::a[] = {.5}; const real_t RK2IDPSolver::b[] = {0., 1.}; const real_t RK2IDPSolver::c[] = {.5}; +//3rd order +const real_t RK3IDPSolver::a[] = {1./3., 0., 2./3.}; +const real_t RK3IDPSolver::b[] = {.25, 0., .75}; +const real_t RK3IDPSolver::c[] = {1./3., 2./3.}; + +//4th order for linear, 3rd for non-linear +const real_t RK4IDPSolver::a[] = {.25, 0., .5, 0., .25, .5}; +const real_t RK4IDPSolver::b[] = {0., 2./3., -1./3., 2./3.}; +const real_t RK4IDPSolver::c[] = {.25, .5, .75}; + } // namespace mfem diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index 3a551a9..540160e 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -115,6 +115,20 @@ class RK2IDPSolver : public RKIDPSolver RK2IDPSolver() : RKIDPSolver(2, a, b, c) { } }; +class RK3IDPSolver : public RKIDPSolver +{ + static const real_t a[], b[], c[]; +public: + RK3IDPSolver() : RKIDPSolver(3, a, b, c) { } +}; + +class RK4IDPSolver : public RKIDPSolver +{ + static const real_t a[], b[], c[]; +public: + RK4IDPSolver() : RKIDPSolver(4, a, b, c) { } +}; + } // namespace mfem #endif // MFEM_REMHOS_SOLVERS From fe5078dee1989690611bb5a9758d17693771a2fc Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 18 Nov 2024 15:39:36 -0800 Subject: [PATCH 19/39] Added support for equal time levels in RK schemes. Changed RK3 to the 3/8 rule. Added RK6, but explodes numerically. --- remhos.cpp | 4 +--- remhos_solvers.cpp | 48 ++++++++++++++++++++++++++++++++++------------ remhos_solvers.hpp | 7 +++++++ 3 files changed, 44 insertions(+), 15 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index ffec017..256545c 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -311,9 +311,7 @@ int main(int argc, char *argv[]) case 2: ode_solver = new RK2IDPSolver(); break; case 3: ode_solver = new RK3IDPSolver(); break; case 4: ode_solver = new RK4IDPSolver(); break; - /*case 6: - if (myid == 0) { MFEM_WARNING("RK6 may violate the bounds."); } - ode_solver = new RK6Solver; break;*/ + case 6: ode_solver = new RK6IDPSolver(); break;//does not work properly! default: cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; return 3; diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 7f50cea..3911ee7 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -72,9 +72,13 @@ void RKIDPSolver::ConstructD() } di[i] = a_n[i] / dc; - i_o = i; - c_o = c_n; - a_o = a_n; + const double c_next = (i < s-2)?(c[i+1]):(1.); + if (c_next > c_n) + { + i_o = i; + c_o = c_n; + a_o = a_n; + } if (i < s-2) { @@ -174,10 +178,10 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) // Step through higher stages const real_t *d_i = d + 1; + real_t c_o = c[0]; for (int i = 1; i < s; i++) { - const real_t c_o = c[i-1]; const real_t c_n = (iLimitMult(x, dxs[i]); // Update the state - f->SetTime(t + c_n * dt); - x.Add(dct, dxs[i]); + const double c_next = (i < s-2)?(c[i+1]):(infinity()); + if (c_next > c_n)// only when advancing after + { + f->SetTime(t + c_n * dt); + x.Add(dct, dxs[i]); + c_o = c_n; + } d_i += i+1; } t += dt; } -//2nd order +//2-stage, 2nd order const real_t RK2IDPSolver::a[] = {.5}; const real_t RK2IDPSolver::b[] = {0., 1.}; const real_t RK2IDPSolver::c[] = {.5}; -//3rd order +//3-stage, 3rd order const real_t RK3IDPSolver::a[] = {1./3., 0., 2./3.}; const real_t RK3IDPSolver::b[] = {.25, 0., .75}; const real_t RK3IDPSolver::c[] = {1./3., 2./3.}; -//4th order for linear, 3rd for non-linear -const real_t RK4IDPSolver::a[] = {.25, 0., .5, 0., .25, .5}; -const real_t RK4IDPSolver::b[] = {0., 2./3., -1./3., 2./3.}; -const real_t RK4IDPSolver::c[] = {.25, .5, .75}; +//4-stage, 4th order for linear, 3rd for non-linear +//const real_t RK4IDPSolver::a[] = {.25, 0., .5, 0., .25, .5}; +//const real_t RK4IDPSolver::b[] = {0., 2./3., -1./3., 2./3.}; +//const real_t RK4IDPSolver::c[] = {.25, .5, .75}; +//4-stage, 4th order, non-equidistant +//const real_t RK4IDPSolver::a[] = {.5, 0., .5, 0., 0., 1.}; +//const real_t RK4IDPSolver::b[] = {1./6., 2./6., 2./6., 1./6.}; +//const real_t RK4IDPSolver::c[] = {.5, .5, 1.}; +//4-stage, 4th order, equidistant (except the last) +const real_t RK4IDPSolver::a[] = {1./3., -1./3., 1., 1., -1., 1.}; +const real_t RK4IDPSolver::b[] = {1./8., 3./8., 3./8., 1./8.}; +const real_t RK4IDPSolver::c[] = {1./3., 2./3., 1.}; + +//6-stage, 5th order, equidistant (except the last) +const real_t RK6IDPSolver::a[] = {.25, 1./8., 1./8., 0., -.5, 1., 3./16., 0., 0., + 9./16., -3./7., 2./7., 12./7., -12./7., 8./7. + }; +const real_t RK6IDPSolver::b[] = {7./90., 0., 32./90., 12./90., 32./90., 7./90.}; +const real_t RK6IDPSolver::c[] = {.25, .25, .5, .75, 1.}; } // namespace mfem diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index 540160e..90d0b5e 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -129,6 +129,13 @@ class RK4IDPSolver : public RKIDPSolver RK4IDPSolver() : RKIDPSolver(4, a, b, c) { } }; +class RK6IDPSolver : public RKIDPSolver +{ + static const real_t a[], b[], c[]; +public: + RK6IDPSolver() : RKIDPSolver(6, a, b, c) { } +}; + } // namespace mfem #endif // MFEM_REMHOS_SOLVERS From c02ce4025e6665c7c8343d49f61ef02c82d993f0 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 18 Nov 2024 16:05:59 -0800 Subject: [PATCH 20/39] Fixed dt type to real_t. --- remhos.cpp | 6 +++--- remhos_fct.hpp | 14 +++++++------- remhos_lo.hpp | 4 ++-- remhos_solvers.hpp | 2 +- 4 files changed, 13 insertions(+), 13 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 256545c..cc4a514 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -95,7 +95,7 @@ class AdvectionOperator : public LimitedTimeDependentOperator mutable ParGridFunction x_gf; TimeStepControl dt_control; - mutable double dt_est; + mutable real_t dt_est; Assembly &asmbl; LowOrderMethod &lom; @@ -136,7 +136,7 @@ class AdvectionOperator : public LimitedTimeDependentOperator dt_control = tsc; } - void SetDt(double _dt) override + void SetDt(real_t _dt) override { LimitedTimeDependentOperator::SetDt(_dt); dt_est = dt; @@ -144,7 +144,7 @@ class AdvectionOperator : public LimitedTimeDependentOperator if (fct_solver) { fct_solver->UpdateTimeStep(dt); } } - double GetTimeStepEstimate() { return dt_est; } + real_t GetTimeStepEstimate() { return dt_est; } void SetRemapStartPos(const Vector &m_pos, const Vector &sm_pos) { diff --git a/remhos_fct.hpp b/remhos_fct.hpp index e3b9f5e..0994093 100644 --- a/remhos_fct.hpp +++ b/remhos_fct.hpp @@ -32,7 +32,7 @@ class FCTSolver protected: ParFiniteElementSpace &pfes; SmoothnessIndicator *smth_indicator; - double dt; + real_t dt; const bool needs_LO_input_for_products; // Computes a compatible slope (piecewise constan = mass_us / mass_u). @@ -51,13 +51,13 @@ class FCTSolver public: FCTSolver(ParFiniteElementSpace &space, - SmoothnessIndicator *si, double dt_, bool needs_LO_prod) + SmoothnessIndicator *si, real_t dt_, bool needs_LO_prod) : pfes(space), smth_indicator(si), dt(dt_), needs_LO_input_for_products(needs_LO_prod) { } virtual ~FCTSolver() { } - virtual void UpdateTimeStep(double dt_new) { dt = dt_new; } + virtual void UpdateTimeStep(real_t dt_new) { dt = dt_new; } bool NeedsLOProductInput() const { return needs_LO_input_for_products; } @@ -110,7 +110,7 @@ class FluxBasedFCT : public FCTSolver public: FluxBasedFCT(ParFiniteElementSpace &space, - SmoothnessIndicator *si, double delta_t, + SmoothnessIndicator *si, real_t delta_t, const SparseMatrix &adv_mat, const Array &adv_smap, const SparseMatrix &mass_mat, int fct_iterations = 1) : FCTSolver(space, si, delta_t, true), @@ -134,7 +134,7 @@ class ClipScaleSolver : public FCTSolver { public: ClipScaleSolver(ParFiniteElementSpace &space, - SmoothnessIndicator *si, double dt) + SmoothnessIndicator *si, real_t dt) : FCTSolver(space, si, dt, false) { } virtual void CalcFCTSolution(const ParGridFunction &u, const Vector &m, @@ -153,7 +153,7 @@ class ClipScaleSolver : public FCTSolver class ElementFCTProjection : public FCTSolver { public: - ElementFCTProjection(ParFiniteElementSpace &space, double dt) + ElementFCTProjection(ParFiniteElementSpace &space, real_t dt) : FCTSolver(space, NULL, dt, false) { } virtual void CalcFCTSolution(const ParGridFunction &u, const Vector &m, @@ -178,7 +178,7 @@ class NonlinearPenaltySolver : public FCTSolver public: NonlinearPenaltySolver(ParFiniteElementSpace &space, - SmoothnessIndicator *si, double dt_) + SmoothnessIndicator *si, real_t dt_) : FCTSolver(space, si, dt_, false) { } virtual void CalcFCTSolution(const ParGridFunction &u, const Vector &m, diff --git a/remhos_lo.hpp b/remhos_lo.hpp index 83a3616..66a1141 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -27,14 +27,14 @@ class LOSolver { protected: ParFiniteElementSpace &pfes; - double dt = -1.0; // usually not known at creation, updated later. + real_t dt = -1.0; // usually not known at creation, updated later. public: LOSolver(ParFiniteElementSpace &space) : pfes(space) { } virtual ~LOSolver() { } - virtual void UpdateTimeStep(double dt_new) { dt = dt_new; } + virtual void UpdateTimeStep(real_t dt_new) { dt = dt_new; } virtual void CalcLOSolution(const Vector &u, Vector &du) const = 0; }; diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index 90d0b5e..8561301 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -40,7 +40,7 @@ class LimitedTimeDependentOperator : public TimeDependentOperator virtual ~LimitedTimeDependentOperator() { } - virtual void SetDt(double dt_) { dt = dt_; } + virtual void SetDt(real_t dt_) { dt = dt_; } virtual real_t GetDt() const { return dt; } void Mult(const Vector &u, Vector &k) const override From 054a198c49e25ba118b58de1e99224e6eaaff8cd Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 18 Nov 2024 16:17:52 -0800 Subject: [PATCH 21/39] Fixed adaptive time step for IDP RK methods. Allows also prolongation of time step. --- remhos.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index cc4a514..b755b63 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -95,7 +95,7 @@ class AdvectionOperator : public LimitedTimeDependentOperator mutable ParGridFunction x_gf; TimeStepControl dt_control; - mutable real_t dt_est; + mutable real_t dt_est, dt_ratio; Assembly &asmbl; LowOrderMethod &lom; @@ -145,6 +145,8 @@ class AdvectionOperator : public LimitedTimeDependentOperator } real_t GetTimeStepEstimate() { return dt_est; } + void ResetTimeStepRatio() { dt_ratio = infinity(); } + real_t GetTimeStepRatio() { return dt_ratio; } void SetRemapStartPos(const Vector &m_pos, const Vector &sm_pos) { @@ -935,6 +937,7 @@ int main(int argc, char *argv[]) // This also resets the time step estimate when automatic dt is on. adv.SetDt(dt_real); + adv.ResetTimeStepRatio(); if (product_sync) { @@ -958,8 +961,8 @@ int main(int argc, char *argv[]) if (dt_control != TimeStepControl::FixedTimeStep) { - double dt_est = adv.GetTimeStepEstimate(); - if (dt_est < dt_real) + real_t dt_ratio = adv.GetTimeStepRatio(); + if (dt_ratio < 1.) { // Repeat with the proper time step. if (myid == 0) @@ -974,7 +977,7 @@ int main(int argc, char *argv[]) if (dt < 1e-12) { MFEM_ABORT("The time step crashed!"); } continue; } - else if (dt_est > 1.25 * dt_real) { dt *= 1.02; } + else if (dt_ratio > 1.25) { dt *= 1.02; } } // S has been modified, update the alias @@ -1583,6 +1586,7 @@ void AdvectionOperator::UpdateTimeStepEstimate(const Vector &x, Kbf.ParFESpace()->GetComm()); dt_est = fmin(dt_est, dt); + dt_ratio = fmin(dt_ratio, (GetDt() != 0.)?(dt / GetDt()):(0.)); } // Velocity coefficient From 2a96acf1fe3aed54a0174283e0f09725f6734559 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 19 Nov 2024 17:28:56 -0800 Subject: [PATCH 22/39] Fixed RK6 and other solvers with repeating time levels. --- remhos_solvers.cpp | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 3911ee7..5649b61 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -165,20 +165,33 @@ void RKIDPSolver::Init(LimitedTimeDependentOperator &f_) void RKIDPSolver::Step(Vector &x, double &t, double &dt) { + real_t c_o = 0.; + // Perform the first step f->SetTime(t); f->SetDt(c[0] * dt); f->MultUnlimited(x, dxs[0]); f->LimitMult(x, dxs[0]); - x.Add(c[0] * dt, dxs[0]); - f->ComputeMask(x, mask); - f->SetTime(c[0] * dt); + // Update state + const double c_next = (s > 2)?(c[1]):(1.); + if (c_next > c[0])// only when advancing after + { + x.Add(c[0] * dt, dxs[0]); + f->ComputeMask(x, mask); + f->SetTime(c[0] * dt); + c_o = c[0]; + } + else + { + Vector x_new(x.Size()); + add(x, c[0] * dt, dxs[0], x_new); + f->ComputeMask(x_new, mask); + } // Step through higher stages const real_t *d_i = d + 1; - real_t c_o = c[0]; for (int i = 1; i < s; i++) { @@ -207,8 +220,8 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) f->LimitMult(x, dxs[i]); // Update the state - const double c_next = (i < s-2)?(c[i+1]):(infinity()); - if (c_next > c_n)// only when advancing after + const double c_next = (i < s-2)?(c[i+1]):(1.); + if (i == s-1 || c_next > c_n)// only when advancing after { f->SetTime(t + c_n * dt); x.Add(dct, dxs[i]); From fb17293bc3541bbea5e519fdb0b858bf3d9f8891 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Tue, 19 Nov 2024 17:56:55 -0800 Subject: [PATCH 23/39] Minor comments. --- remhos.cpp | 2 +- remhos_solvers.cpp | 22 ++++++++++++++-------- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index b755b63..478b2a1 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -313,7 +313,7 @@ int main(int argc, char *argv[]) case 2: ode_solver = new RK2IDPSolver(); break; case 3: ode_solver = new RK3IDPSolver(); break; case 4: ode_solver = new RK4IDPSolver(); break; - case 6: ode_solver = new RK6IDPSolver(); break;//does not work properly! + case 6: ode_solver = new RK6IDPSolver(); break; default: cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; return 3; diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 5649b61..51d3938 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -42,26 +42,27 @@ void RKIDPSolver::ConstructD() // Convert high-order to Forward Euler factors d = new real_t[s*(s+1)/2]; - const real_t *a_n = a; - const real_t *a_o = a; - int i_o = -1; - real_t c_o = 0.; + const real_t *a_n = a; // new coeff line + const real_t *a_o = a; // old coeff line + int i_o = -1; // old stage + real_t c_o = 0.; // old time fraction for (int i = 0; i < s; i++) { - const real_t c_n = (i c_n) { @@ -135,6 +138,8 @@ void RKIDPSolver::UpdateMask(const Vector &x, const Vector &dx, real_t dt, { f->ComputeMask(x, mask_new); } + // All intermediate updates must be on active DOFs + // for a valid high-order update for (int i = 0; i < mask.Size(); i++) { mask[i] = mask[i] && mask_new[i]; @@ -184,6 +189,7 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) } else { + // Only initialize the mask Vector x_new(x.Size()); add(x, c[0] * dt, dxs[0], x_new); f->ComputeMask(x_new, mask); From 15d61567a8d965feff89549114864298068d3301 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Sun, 24 Nov 2024 16:35:15 -0800 Subject: [PATCH 24/39] minor --- remhos.cpp | 5 ++++- remhos_solvers.hpp | 3 +++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/remhos.cpp b/remhos.cpp index 478b2a1..93ea0f3 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -1436,10 +1436,12 @@ void AdvectionOperator::ComputeMask(const Vector &x, Array &mask) const } if (!dofs_active) + { for (int j = 0; j < ndof_el; j++) { bool_dofs[j+ndof_el*k] = false; } + } } // Apply the u mask to all product fields @@ -1447,10 +1449,12 @@ void AdvectionOperator::ComputeMask(const Vector &x, Array &mask) const const int ndim = block_offsets.Size() - 1; mask.SetSize(ndofs * ndim); for (int i = 0; i < ndofs; i++) + { for (int d = 0; d < ndim; d++) { mask[i+ndofs*d] = bool_dofs[i]; } + } } void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const @@ -1463,7 +1467,6 @@ void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const const Vector &u = block_X.GetBlock(0); Vector &d_u = block_Y.GetBlock(0); - if (fct_solver) { MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers."); diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index 8561301..b37a96f 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -57,6 +57,8 @@ class LimitedTimeDependentOperator : public TimeDependentOperator { MFEM_ABORT("Mask computation not implemented!"); } /// Limit the action vector @a k + /// 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; }; @@ -94,6 +96,7 @@ class RKIDPSolver : public IDPODESolver void ConstructD(); + /// Adds only DOFs that have mask = true. void AddMasked(const Array &mask, real_t b, const Vector &vb, Vector &va); From 9d24077976c17582957bb6ced5e6fedb5f403c07 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Sun, 24 Nov 2024 22:54:11 -0800 Subject: [PATCH 25/39] attempt to fix autotest --- .github/workflows/build-and-test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index 1105b3e..e3c68aa 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -115,7 +115,7 @@ jobs: - name: Archive test results patch if: always() - uses: actions/upload-artifact@v2 + uses: actions/upload-artifact@v3 with: name: baseline-patch path: remhos/baseline.patch From 37b9ce2db391c72d4699b14114b5a27c84f5bf52 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 25 Nov 2024 10:08:46 -0800 Subject: [PATCH 26/39] Fixed typo. --- remhos_solvers.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 51d3938..9c9f4cf 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -184,7 +184,7 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) { x.Add(c[0] * dt, dxs[0]); f->ComputeMask(x, mask); - f->SetTime(c[0] * dt); + f->SetTime(t + c[0] * dt); c_o = c[0]; } else From 6c333402a2bd80afc6b7b934affa8bbe3428574c Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Wed, 4 Dec 2024 18:18:42 -0800 Subject: [PATCH 27/39] Restored the options for the old time integrators, added some comments. --- autotest/out_baseline.dat | 18 ++++++++-------- autotest/test.sh | 6 +++--- remhos.cpp | 43 +++++++++++++++++++++++++-------------- remhos_solvers.cpp | 8 +++++++- remhos_solvers.hpp | 10 ++++++++- 5 files changed, 56 insertions(+), 29 deletions(-) diff --git a/autotest/out_baseline.dat b/autotest/out_baseline.dat index 0928712..8cdebe5 100644 --- a/autotest/out_baseline.dat +++ b/autotest/out_baseline.dat @@ -185,19 +185,19 @@ Final mass u: 0.8143001155 Max value u: 0.9995762608 --- Product remap 2D (FCT) -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -Final mass us: 0.1796082878 -Mass loss us: 4.15262e-07 +mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 13 +Final mass us: 0.1767475452 +Mass loss us: 0.00286033 --- Product remap 2D (ClipScale) -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 1 -Final mass us: 0.1814773886 -Mass loss us: 0.00186952 +mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 13 +Final mass us: 0.1782170448 +Mass loss us: 0.00139083 --- Product remap 2D (FCTProject) -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 1 -Final mass us: 0.1814920761 -Mass loss us: 0.0018842 +mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 13 +Final mass us: 0.1776386793 +Mass loss us: 0.00196919 --- BLAST sharpening test - Pacman remap auto-dt mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt -1 -tf 0.75 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1 diff --git a/autotest/test.sh b/autotest/test.sh index 70cd450..79029fc 100755 --- a/autotest/test.sh +++ b/autotest/test.sh @@ -66,17 +66,17 @@ for method in "${methods[@]}"; do done echo -e '\n'"--- Product remap 2D (FCT)" >> $file -run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps" +run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 13" echo -e $run_line >> $file $run_line | grep -e 'mass us' -e 'loss us'>> $file echo -e '\n'"--- Product remap 2D (ClipScale)" >> $file -run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 1" +run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 13" echo -e $run_line >> $file $run_line | grep -e 'mass us' -e 'loss us'>> $file echo -e '\n'"--- Product remap 2D (FCTProject)" >> $file -run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 1" +run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 13" echo -e $run_line >> $file $run_line | grep -e 'mass us' -e 'loss us'>> $file diff --git a/remhos.cpp b/remhos.cpp index 93ea0f3..054ae73 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -47,16 +47,20 @@ using namespace std; using namespace mfem; -enum class HOSolverType {None, Neumann, CG, LocalInverse}; -enum class FCTSolverType {None, FluxBased, ClipScale, - NonlinearPenalty, FCTProject - }; -enum class LOSolverType {None, DiscrUpwind, DiscrUpwindPrec, - ResDist, ResDistSubcell, MassBased - }; -enum class MonolithicSolverType {None, ResDistMono, ResDistMonoSubcell}; +enum class HOSolverType +{ None, Neumann, CG, LocalInverse }; -enum class TimeStepControl {FixedTimeStep, LOBoundsError}; +enum class LOSolverType +{ None, DiscrUpwind, DiscrUpwindPrec, ResDist, ResDistSubcell, MassBased }; + +enum class FCTSolverType +{ None, FluxBased, ClipScale, NonlinearPenalty, FCTProject }; + +enum class MonolithicSolverType +{ None, ResDistMono, ResDistMonoSubcell }; + +enum class TimeStepControl +{ FixedTimeStep, LOBoundsError }; // Choice for the problem setup. The fluid velocity, initial condition and // inflow boundary condition are chosen based on this parameter. @@ -306,18 +310,25 @@ int main(int argc, char *argv[]) // Define the ODE solver used for time integration. Several explicit // Runge-Kutta methods are available. - IDPODESolver *ode_solver = NULL; + IDPODESolver *idp_ode_solver = nullptr; + ODESolver *ode_solver = nullptr; switch (ode_solver_type) { - case 1: ode_solver = new ForwardEulerIDPSolver(); break; - case 2: ode_solver = new RK2IDPSolver(); break; - case 3: ode_solver = new RK3IDPSolver(); break; - case 4: ode_solver = new RK4IDPSolver(); break; - case 6: ode_solver = new RK6IDPSolver(); break; + case 1: ode_solver = new ForwardEulerSolver(); break; + case 2: ode_solver = new RK2Solver(1.0); break; + case 3: ode_solver = new RK3SSPSolver(); break; + case 4: ode_solver = new RK4Solver(); break; + case 6: ode_solver = new RK6Solver(); break; + 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; default: cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; return 3; } + if (ode_solver_type > 10) { ode_solver = idp_ode_solver; } // Check if the input mesh is periodic. const bool periodic = pmesh.GetNodes() != NULL && @@ -1452,6 +1463,8 @@ void AdvectionOperator::ComputeMask(const Vector &x, Array &mask) const { for (int d = 0; d < ndim; d++) { + // Note that this puts a mask on the first field as well, meaning + // that propagation in new elements will be through Forward Euler. mask[i+ndofs*d] = bool_dofs[i]; } } diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 9c9f4cf..43a40d1 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -212,13 +212,19 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) // Update mask with the HO update UpdateMask(x, dxs[i], dct, mask); - // Convert the HO update to Forward Euler + // Form the unlimited update for the stage. + // Note that it converts eq. (2.16) in JLG's paper into an update using + // the previous limited updates. if (d_i[i] != 1.) { + // for mask = 0, we get dxs (nothing happens). + // the loop below won't change it -> Forward Euler. + // for mask = 1, we scale dxs by d_i[i]. AddMasked(mask, d_i[i]-1., dxs[i], dxs[i]); } for (int j = 0; j < i; j++) { + // Use all previous limited updates. AddMasked(mask, d_i[j], dxs[j], dxs[i]); } diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index b37a96f..d3a983d 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -69,7 +69,12 @@ class IDPODESolver : public ODESolver LimitedTimeDependentOperator *f; // f(.,t) : R^n --> R^n void Init(TimeDependentOperator &f_) override - { MFEM_ABORT("Limited time-dependent operator must be assigned!"); } + { + auto lo = dynamic_cast(&f_); + if (lo) { Init(*lo); } + else { MFEM_ABORT("LimitedTimeDependentOperator must be assigned!"); } + } + public: IDPODESolver() : ODESolver(), f(NULL) { } virtual ~IDPODESolver() { } @@ -94,6 +99,9 @@ class RKIDPSolver : public IDPODESolver Vector *dxs; Array mask; + // This function constructs coefficients that transform eq. (2.16) from + // JLG's paper to an update that only uses the previous limited updates. + // This function does not depend on the Operator f in any way. void ConstructD(); /// Adds only DOFs that have mask = true. From 4803c5b2d3931316376b17a948efd4e61c2af487 Mon Sep 17 00:00:00 2001 From: Jan Nikl Date: Mon, 9 Dec 2024 09:36:54 -0800 Subject: [PATCH 28/39] Fixed the final mesh position in remap. --- remhos.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/remhos.cpp b/remhos.cpp index 054ae73..00cbea9 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -1132,6 +1132,14 @@ int main(int argc, char *argv[]) << " (" << ti_total-ti << " repeated)." << endl; } + // Move to the final mesh position + if (exec_mode == 1) + { + add(x0, t_final, v_gf, x); + // Reset precomputed geometric data. + pmesh.DeleteGeometricFactors(); + } + // Print the final meshes and solution. { ofstream meshHO("meshHO_final.mesh"); From 47829569a402819f9aab74ecf71ec6c8095dcdba Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Tue, 17 Dec 2024 16:54:08 -0800 Subject: [PATCH 29/39] Improved the verify_bounds options, added additional checks. MassBasedAvgLOSolver can work with prescribed HO solution. --- remhos.cpp | 152 +++++++++++++++++++++++++++++++-------------- remhos_fct.cpp | 96 ++++++++++++++-------------- remhos_fct.hpp | 4 +- remhos_lo.cpp | 25 ++++---- remhos_lo.hpp | 6 ++ remhos_solvers.cpp | 2 + remhos_sync.cpp | 3 + 7 files changed, 180 insertions(+), 108 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index 00cbea9..d96e90b 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -158,6 +158,8 @@ class AdvectionOperator : public LimitedTimeDependentOperator start_submesh_pos = sm_pos; } + bool verify_bounds = false; + virtual ~AdvectionOperator() { } }; @@ -917,13 +919,16 @@ int main(int argc, char *argv[]) x, xsub, v_gf, v_sub_gf, asmbl, lom, dofs, ho_solver, lo_solver, fct_solver, mono_solver); + adv.verify_bounds = verify_bounds; + fct_solver->verify_bounds = verify_bounds; + double t = 0.0; adv.SetTime(t); adv.SetTimeStepControl(dt_control); ode_solver->Init(adv); - double umin, umax; - GetMinMax(u, umin, umax); + double u_min, u_max; + GetMinMax(u, u_min, u_max); if (exec_mode == 1) { @@ -935,8 +940,8 @@ int main(int argc, char *argv[]) ParGridFunction res = u; double residual = 0.0; - double s_min_glob = numeric_limits::infinity(), - s_max_glob = -numeric_limits::infinity(); + double s_min = numeric_limits::infinity(), + s_max = -numeric_limits::infinity(); // Time-integration (loop over the time iterations, ti, with a time-step dt). bool done = false; @@ -950,20 +955,21 @@ int main(int argc, char *argv[]) adv.SetDt(dt_real); adv.ResetTimeStepRatio(); +#ifdef REMHOS_FCT_PRODUCT_DEBUG if (product_sync) { - ComputeMinMaxS(pmesh.GetNE(), us, u, s_min_glob, s_max_glob); -#ifdef REMHOS_FCT_PRODUCT_DEBUG + double s_min_debug, s_max_debug; + ComputeMinMaxS(pmesh.GetNE(), us, u, s_min_debug, s_max_debug); if (myid == 0) { std::cout << " --- Full time step" << std::endl; std::cout << " in: "; std::cout << std::scientific << std::setprecision(5); - std::cout << "min_s: " << s_min_glob - << "; max_s: " << s_max_glob << std::endl; + std::cout << "min_s: " << s_min_debug + << "; max_s: " << s_max_debug << std::endl; } -#endif } +#endif Sold = S; ode_solver->Step(S, t, dt_real); @@ -993,38 +999,19 @@ int main(int argc, char *argv[]) // S has been modified, update the alias u.SyncMemory(S); + if (product_sync) { us.SyncMemory(S); - - // It is known that RK time integrators with more than 1 stage may - // cause violation of the lower bounds for us. - // The lower bound is corrected, causing small conservation error. - // Correction can also be done with localized bounds for s, but for - // now we have implemented only the minimum global bound. - u.HostRead(); - us.HostReadWrite(); -#if 0 - const int s = u.Size(); - Array active_elem, active_dofs; - ComputeBoolIndicators(NE, u, active_elem, active_dofs); - for (int i = 0; i < s; i++) - { - if (active_dofs[i] == false) { continue; } - - double us_min = u(i) * s_min_glob; - if (us(i) < us_min) { us(i) = us_min; } - } -#endif - #ifdef REMHOS_FCT_PRODUCT_DEBUG - ComputeMinMaxS(NE, us, u, s_min_glob, s_max_glob); + double s_min_debug, s_max_debug; + ComputeMinMaxS(NE, us, u, s_min_debug, s_max_debug); if (myid == 0) { std::cout << " out: "; std::cout << std::scientific << std::setprecision(5); - std::cout << "min_s: " << s_min_glob - << "; max_s: " << s_max_glob << std::endl; + std::cout << "min_s: " << s_min_debug + << "; max_s: " << s_max_debug << std::endl; } #endif } @@ -1032,28 +1019,42 @@ int main(int argc, char *argv[]) // Monotonicity check for debug purposes mainly. if (verify_bounds && forced_bounds && smth_indicator == NULL) { - double umin_new, umax_new; - GetMinMax(u, umin_new, umax_new); + double u_min_new, u_max_new, + s_min_new = s_min, s_max_new = s_max; + GetMinMax(u, u_min_new, u_max_new); + if (product_sync) + { + ComputeMinMaxS(NE, us, u, s_min_new, s_max_new); + } + if (problem_num % 10 != 6 && problem_num % 10 != 7) { if (myid == 0) { - MFEM_VERIFY(umin_new > umin - 1e-12, - "Undershoot of " << umin - umin_new); - MFEM_VERIFY(umax_new < umax + 1e-12, - "Overshoot of " << umax_new - umax); + MFEM_VERIFY(u_min_new > u_min - 1e-12, + "Undershoot of " << u_min - u_min_new); + MFEM_VERIFY(u_max_new < u_max + 1e-12, + "Overshoot of " << u_max_new - u_max); + MFEM_VERIFY(s_min_new > s_min - 1e-12, + "Undershoot in s of " << s_min - s_min_new); + MFEM_VERIFY(s_max_new < s_max + 1e-12, + "Overshoot in s of " << s_max_new - s_max); } - umin = umin_new; - umax = umax_new; + u_min = u_min_new; + u_max = u_max_new; } else { if (myid == 0) { - MFEM_VERIFY(umin_new > 0.0 - 1e-12, - "Undershoot of " << 0.0 - umin_new); - MFEM_VERIFY(umax_new < 1.0 + 1e-12, - "Overshoot of " << umax_new - 1.0); + MFEM_VERIFY(u_min_new > 0.0 - 1e-12, + "Undershoot of " << 0.0 - u_min_new); + MFEM_VERIFY(u_max_new < 1.0 + 1e-12, + "Overshoot of " << u_max_new - 1.0); + MFEM_VERIFY(s_min_new > 0.0 - 1e-12, + "Undershoot in s of " << 0.0 - s_min_new); + MFEM_VERIFY(s_max_new < 1.0 + 1e-12, + "Overshoot in s of " << s_max_new - 1.0); } } } @@ -1172,10 +1173,10 @@ int main(int argc, char *argv[]) mass_u_loc = masses * u; if (product_sync) { mass_us_loc = masses * us; } } - double mass_u, mass_us = 0.0, s_max = 0.0; + double mass_u, mass_us = 0.0; MPI_Allreduce(&mass_u_loc, &mass_u, 1, MPI_DOUBLE, MPI_SUM, comm); const double umax_loc = u.Max(); - MPI_Allreduce(&umax_loc, &umax, 1, MPI_DOUBLE, MPI_MAX, comm); + MPI_Allreduce(&umax_loc, &u_max, 1, MPI_DOUBLE, MPI_MAX, comm); if (product_sync) { ComputeRatio(pmesh.GetNE(), us, u, s, u_bool_el, u_bool_dofs); @@ -1187,7 +1188,7 @@ int main(int argc, char *argv[]) { cout << setprecision(10) << "Final mass u: " << mass_u << endl - << "Max value u: " << umax << endl << setprecision(6) + << "Max value u: " << u_max << endl << setprecision(6) << "Mass loss u: " << abs(mass0_u - mass_u) << endl; if (product_sync) { @@ -1292,6 +1293,45 @@ AdvectionOperator::AdvectionOperator(const Array &offsets, asmbl(_asmbl), lom(_lom), dofs(_dofs), ho_solver(hos), lo_solver(los), fct_solver(fct), mono_solver(mos) { } +void check_violation(const Vector &u_new, + const Vector &u_min, const Vector &u_max, + string info, double tol, const Array *active_dofs) +{ + const int size = u_new.Size(); + for (int i = 0; i < size; i++) + { + if (active_dofs && (*active_dofs)[i] == false) { continue; } + + if (u_new(i) + tol < u_min(i) || u_new(i) > u_max(i) + tol) + { + cout << info << " bounds violation: " << i << " " + << u_min(i) << " " << u_new(i) << " " << u_max(i) << endl; + cout << u_max(i) - u_new(i) << " " << u_new(i) - u_min(i) << endl; + MFEM_ABORT("Aborted due to bounds violation."); + } + } +} + +void check_violation(const Vector &u, double dt, const Vector &du_new, + const Vector &u_min, const Vector &u_max, + string info, double tol, const Array *active_dofs) +{ + const int size = u.Size(); + for (int i = 0; i < size; i++) + { + if (active_dofs && (*active_dofs)[i] == false) { continue; } + + const double u_new = u(i) + dt * du_new(i); + if (u_new + tol < u_min(i) || u_new > u_max(i) + tol) + { + cout << info << " bounds violation: " << i << " " + << u_min(i) << " " << u_new << " " << u_max(i) << endl; + cout << u_max(i) - u_new << " " << u_new - u_min(i) << endl; + MFEM_ABORT("Aborted due to bounds violation."); + } + } +} + void AdvectionOperator::MultUnlimited(const Vector &X, Vector &Y) const { if (exec_mode == 1) @@ -1497,13 +1537,28 @@ void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const Vector du_HO(d_u), du_LO(u.Size()); + auto mba = dynamic_cast(lo_solver); + if (mba) { mba->SetHOSolution(du_HO); } lo_solver->CalcLOSolution(u, du_LO); 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); + + if (verify_bounds) + { + check_violation(u, dt, du_LO, dofs.xi_min, dofs.xi_max, + "LimitMult LO u", 1e-12, nullptr); + } + fct_solver->CalcFCTSolution(x_gf, lumpedM, du_HO, du_LO, dofs.xi_min, dofs.xi_max, d_u); + if (verify_bounds) + { + check_violation(u, dt, d_u, dofs.xi_min, dofs.xi_max, + "LimitMult FCT solution u", 1e-12, nullptr); + } + if (dt_control == TimeStepControl::LOBoundsError) { UpdateTimeStepEstimate(u, du_LO, dofs.xi_min, dofs.xi_max); @@ -1530,6 +1585,7 @@ void AdvectionOperator::LimitMult(const Vector &X, Vector &Y) const x_gf.ExchangeFaceNbrData(); Vector d_us_HO(d_us), d_us_LO; + if (fct_solver->NeedsLOProductInput()) { d_us_LO.SetSize(us.Size()); diff --git a/remhos_fct.cpp b/remhos_fct.cpp index a8e6d96..46b021c 100644 --- a/remhos_fct.cpp +++ b/remhos_fct.cpp @@ -81,23 +81,24 @@ void FCTSolver::CalcCompatibleLOProduct(const ParGridFunction &us, if (s_avg > smax && mass_us - eps < smax * mass_u) { s_avg = smax; } -#ifdef REMHOS_FCT_PRODUCT_DEBUG - // Check if s_avg = mass_us / mass_u is within the bounds of the full - // stencil of active dofs. - if (mass_us + eps < smin * mass_u || - mass_us - eps > smax * mass_u || - s_avg + eps < smin || - s_avg - eps > smax) + if (verify_bounds) { - std::cout << "---\ns_avg element bounds: " - << smin << " " << s_avg << " " << smax << std::endl; - std::cout << "Element " << k << std::endl; - std::cout << "Masses " << mass_us << " " << mass_u << std::endl; - PrintCellValues(k, NE, u_new, "u_loc: "); + // Check if s_avg = mass_us / mass_u is within the bounds of the + // full stencil of active dofs. + if (mass_us + eps < smin * mass_u || + mass_us - eps > smax * mass_u || + s_avg + eps < smin || + s_avg - eps > smax) + { + std::cout << "---\ns_avg element bounds: " + << smin << " " << s_avg << " " << smax << std::endl; + std::cout << "Element " << k << std::endl; + std::cout << "Masses " << mass_us << " " << mass_u << std::endl; + PrintCellValues(k, NE, u_new, "u_loc: "); - MFEM_ABORT("s_avg is not in the full stencil bounds!"); + MFEM_ABORT("s_avg is not in the full stencil bounds!"); + } } -#endif // When s_avg is not in the local bounds for some dof (it should be // within the full stencil of active dofs), reset the bounds to s_avg. @@ -583,46 +584,49 @@ void ClipScaleSolver::CalcFCTProduct(const ParGridFunction &us, const Vector &m, CalcFCTSolution(us, m, d_us_HO, dus_lo_fct, us_min, us_max, d_us); ZeroOutEmptyDofs(active_el, active_dofs, d_us); -#ifdef REMHOS_FCT_PRODUCT_DEBUG - // Check the bounds of the final solution. - const int NE = pfes.GetNE(); - const int ndofs = u_new.Size() / NE; - int dof_id; - const double eps = 1e-12; - Vector us_new(d_us.Size()); - add(1.0, us, dt, d_us, us_new); - for (int k = 0; k < NE; k++) + if (verify_bounds) { - if (active_el[k] == false) { continue; } - - for (int j = 0; j < ndofs; j++) + // Check the bounds of the final solution. + const int NE = pfes.GetNE(); + const int ndofs = u_new.Size() / NE; + int dof_id; + const double eps = 1e-12; + Vector us_new(d_us.Size()); + add(1.0, us, dt, d_us, us_new); + for (int k = 0; k < NE; k++) { - dof_id = k*ndofs + j; - if (active_dofs[dof_id] == false) { continue; } + if (active_el[k] == false) { continue; } - double s = us_new(dof_id) / u_new(dof_id); - if (s + eps < s_min(dof_id) || - s - eps > s_max(dof_id)) + for (int j = 0; j < ndofs; j++) { - std::cout << "Final s " << j << " " << k << " " - << s_min(dof_id) << " " - << s << " " - << s_max(dof_id) << std::endl; - std::cout << "---\n"; - } + dof_id = k*ndofs + j; + if (active_dofs[dof_id] == false) { continue; } - if (us_new(dof_id) + eps < us_min(dof_id) || - us_new(dof_id) - eps > us_max(dof_id)) - { - std::cout << "Final us " << j << " " << k << " " - << us_min(dof_id) << " " - << us_new(dof_id) << " " - << us_max(dof_id) << std::endl; - std::cout << "---\n"; + /* // this doesn't check round-offs in the division. + double s = us_new(dof_id) / u_new(dof_id); + if (s + eps < s_min(dof_id) || + s - eps > s_max(dof_id)) + { + std::cout << "Final s " << j << " " << k << " " + << s_min(dof_id) << " " + << s << " " + << s_max(dof_id) << std::endl; + std::cout << "---\n"; + }*/ + + if (us_new(dof_id) + eps < us_min(dof_id) || + us_new(dof_id) - eps > us_max(dof_id)) + { + std::cout << "Final us " << j << " " << k << " " + << us_min(dof_id) << " " + << us_new(dof_id) << " " + << us_max(dof_id) << std::endl; + std::cout << "---\n"; + MFEM_ABORT("Bounds violation FCT us."); + } } } } -#endif } void ElementFCTProjection::CalcFCTSolution(const ParGridFunction &u, diff --git a/remhos_fct.hpp b/remhos_fct.hpp index 0994093..8e36438 100644 --- a/remhos_fct.hpp +++ b/remhos_fct.hpp @@ -17,7 +17,7 @@ #ifndef MFEM_REMHOS_FCT #define MFEM_REMHOS_FCT -//#define REMHOS_FCT_PRODUCT_DEBUG +#define REMHOS_FCT_PRODUCT_DEBUG #include "mfem.hpp" @@ -83,6 +83,8 @@ class FCTSolver { MFEM_ABORT("Product remap is not implemented for the chosen solver"); } + + bool verify_bounds = false; }; class FluxBasedFCT : public FCTSolver diff --git a/remhos_lo.cpp b/remhos_lo.cpp index efff529..6150c36 100644 --- a/remhos_lo.cpp +++ b/remhos_lo.cpp @@ -247,25 +247,24 @@ void ResidualDistribution::CalcLOSolution(const Vector &u, Vector &du) const void MassBasedAvg::CalcLOSolution(const Vector &u, Vector &du) const { // Compute the new HO solution. - Vector du_HO(u.Size()); ParGridFunction u_HO_new(&pfes); - ho_solver.CalcHOSolution(u, du_HO); - add(1.0, u, dt, du_HO, u_HO_new); - - // Mesh positions for the new HO solution. - ParMesh *pmesh = pfes.GetParMesh(); - GridFunction x_new(pmesh->GetNodes()->FESpace()); - // Copy the current nodes into x. - pmesh->GetNodes(x_new); - if (mesh_v) + if (du_HO) { - // Remap mode - get the positions of the mesh at time [t + dt]. - x_new.Add(dt, *mesh_v); + add(1.0, u, dt, *du_HO, u_HO_new); + du_HO = nullptr; + } + else + { + Vector du_HO(u.Size()); + ho_solver.CalcHOSolution(u, du_HO); + add(1.0, u, dt, du_HO, u_HO_new); } + // Mesh positions. + ParMesh *pmesh = pfes.GetParMesh(); const int NE = pfes.GetNE(); Vector el_mass(NE), el_vol(NE); - MassesAndVolumesAtPosition(u_HO_new, x_new, el_mass, el_vol); + MassesAndVolumesAtPosition(u_HO_new, *pmesh->GetNodes(), el_mass, el_vol); const int ndofs = u.Size() / NE; for (int k = 0; k < NE; k++) diff --git a/remhos_lo.hpp b/remhos_lo.hpp index 66a1141..76011e9 100644 --- a/remhos_lo.hpp +++ b/remhos_lo.hpp @@ -86,6 +86,9 @@ class MassBasedAvg : public LOSolver HOSolver &ho_solver; const GridFunction *mesh_v; + // Temporary HO solution, used only in the next call to CalcLOSolution(). + mutable const Vector *du_HO = nullptr; + void MassesAndVolumesAtPosition(const ParGridFunction &u, const GridFunction &x, Vector &el_mass, Vector &el_vol) const; @@ -95,6 +98,9 @@ class MassBasedAvg : public LOSolver const GridFunction *mesh_vel) : LOSolver(space), ho_solver(hos), mesh_v(mesh_vel) { } + // Temporary HO solution, used only in the next call to CalcLOSolution(). + void SetHOSolution(Vector &du) { du_HO = &du; } + virtual void CalcLOSolution(const Vector &u, Vector &du) const; }; diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 43a40d1..61a7d0d 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -194,6 +194,7 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) add(x, c[0] * dt, dxs[0], x_new); f->ComputeMask(x_new, mask); } + //mask = 1; // Step through higher stages @@ -211,6 +212,7 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) // Update mask with the HO update UpdateMask(x, dxs[i], dct, mask); + //mask = 1; // Form the unlimited update for the stage. // Note that it converts eq. (2.16) in JLG's paper into an update using diff --git a/remhos_sync.cpp b/remhos_sync.cpp index 30bd605..7ff7039 100644 --- a/remhos_sync.cpp +++ b/remhos_sync.cpp @@ -116,6 +116,9 @@ void ZeroOutEmptyDofs(const Array &ind_elem, void ComputeMinMaxS(int NE, const Vector &us, const Vector &u, double &s_min_glob, double &s_max_glob) { + u.HostRead(); + us.HostRead(); + const int size = u.Size(); Vector s(size); Array bool_el, bool_dofs; From e7a73109eb5b7ad84d257e7383cbd1e6a55a4cca Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Wed, 18 Dec 2024 12:34:09 -0800 Subject: [PATCH 30/39] Option to use boolean masks or not. --- remhos.cpp | 12 +++++++++--- remhos_fct.hpp | 2 +- remhos_solvers.cpp | 14 +++++++------- remhos_solvers.hpp | 3 +++ 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index d96e90b..e9c11e7 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -330,7 +330,12 @@ int main(int argc, char *argv[]) cout << "Unknown ODE solver type: " << ode_solver_type << '\n'; return 3; } - if (ode_solver_type > 10) { ode_solver = idp_ode_solver; } + if (ode_solver_type > 10) + { + auto rk_idp = dynamic_cast(idp_ode_solver); + if (rk_idp) { rk_idp->UseMask(true); } + ode_solver = idp_ode_solver; + } // Check if the input mesh is periodic. const bool periodic = pmesh.GetNodes() != NULL && @@ -929,6 +934,9 @@ int main(int argc, char *argv[]) double u_min, u_max; GetMinMax(u, u_min, u_max); + double s_min = numeric_limits::infinity(), + s_max = -numeric_limits::infinity(); + if (product_sync) { ComputeMinMaxS(NE, us, u, s_min, s_max); } if (exec_mode == 1) { @@ -940,8 +948,6 @@ int main(int argc, char *argv[]) ParGridFunction res = u; double residual = 0.0; - double s_min = numeric_limits::infinity(), - s_max = -numeric_limits::infinity(); // Time-integration (loop over the time iterations, ti, with a time-step dt). bool done = false; diff --git a/remhos_fct.hpp b/remhos_fct.hpp index 8e36438..5b14bf5 100644 --- a/remhos_fct.hpp +++ b/remhos_fct.hpp @@ -17,7 +17,7 @@ #ifndef MFEM_REMHOS_FCT #define MFEM_REMHOS_FCT -#define REMHOS_FCT_PRODUCT_DEBUG +//#define REMHOS_FCT_PRODUCT_DEBUG #include "mfem.hpp" diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 61a7d0d..3fad883 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -183,7 +183,7 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) if (c_next > c[0])// only when advancing after { x.Add(c[0] * dt, dxs[0]); - f->ComputeMask(x, mask); + if (use_masks) { f->ComputeMask(x, mask); } f->SetTime(t + c[0] * dt); c_o = c[0]; } @@ -192,9 +192,8 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) // Only initialize the mask Vector x_new(x.Size()); add(x, c[0] * dt, dxs[0], x_new); - f->ComputeMask(x_new, mask); + if (use_masks) { f->ComputeMask(x_new, mask); } } - //mask = 1; // Step through higher stages @@ -211,8 +210,7 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) f->MultUnlimited(x, dxs[i]); // Update mask with the HO update - UpdateMask(x, dxs[i], dct, mask); - //mask = 1; + if (use_masks) { UpdateMask(x, dxs[i], dct, mask); } // Form the unlimited update for the stage. // Note that it converts eq. (2.16) in JLG's paper into an update using @@ -222,12 +220,14 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) // for mask = 0, we get dxs (nothing happens). // the loop below won't change it -> Forward Euler. // for mask = 1, we scale dxs by d_i[i]. - AddMasked(mask, d_i[i]-1., dxs[i], dxs[i]); + if (use_masks) { AddMasked(mask, d_i[i]-1., dxs[i], dxs[i]); } + else { dxs[i] *= d_i[i]; } } for (int j = 0; j < i; j++) { // Use all previous limited updates. - AddMasked(mask, d_i[j], dxs[j], dxs[i]); + if (use_masks) { AddMasked(mask, d_i[j], dxs[j], dxs[i]); } + else { dxs[i].Add(d_i[j], dxs[j]); } } // Limit the step diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp index d3a983d..9cd4503 100644 --- a/remhos_solvers.hpp +++ b/remhos_solvers.hpp @@ -97,6 +97,7 @@ class RKIDPSolver : public IDPODESolver const real_t *a, *b, *c; real_t *d; Vector *dxs; + bool use_masks = false; Array mask; // This function constructs coefficients that transform eq. (2.16) from @@ -115,6 +116,8 @@ class RKIDPSolver : public IDPODESolver RKIDPSolver(int s_, const real_t a_[], const real_t b_[], const real_t c_[]); ~RKIDPSolver(); + void UseMask(bool mask_on) { use_masks = mask_on; } + void Init(LimitedTimeDependentOperator &f) override; void Step(Vector &x, double &t, double &dt) override; }; From 1b148eeda648284475f0b178e27f266cd0bcb687 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Wed, 18 Dec 2024 15:23:33 -0800 Subject: [PATCH 31/39] comments and other minor. --- remhos.cpp | 22 ++++++------ remhos_fct.cpp | 89 ++++++++++++++++++---------------------------- remhos_solvers.cpp | 7 ++-- 3 files changed, 49 insertions(+), 69 deletions(-) diff --git a/remhos.cpp b/remhos.cpp index e9c11e7..11c5cd0 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -333,7 +333,7 @@ int main(int argc, char *argv[]) if (ode_solver_type > 10) { auto rk_idp = dynamic_cast(idp_ode_solver); - if (rk_idp) { rk_idp->UseMask(true); } + if (rk_idp) { rk_idp->UseMask(false); } ode_solver = idp_ode_solver; } @@ -934,8 +934,7 @@ int main(int argc, char *argv[]) double u_min, u_max; GetMinMax(u, u_min, u_max); - double s_min = numeric_limits::infinity(), - s_max = -numeric_limits::infinity(); + double s_min = -1.0, s_max = -1.0; if (product_sync) { ComputeMinMaxS(NE, us, u, s_min, s_max); } if (exec_mode == 1) @@ -1025,6 +1024,7 @@ int main(int argc, char *argv[]) // Monotonicity check for debug purposes mainly. if (verify_bounds && forced_bounds && smth_indicator == NULL) { + const double eps = 1e-10; double u_min_new, u_max_new, s_min_new = s_min, s_max_new = s_max; GetMinMax(u, u_min_new, u_max_new); @@ -1037,13 +1037,13 @@ int main(int argc, char *argv[]) { if (myid == 0) { - MFEM_VERIFY(u_min_new > u_min - 1e-12, + MFEM_VERIFY(u_min_new > u_min - eps, "Undershoot of " << u_min - u_min_new); - MFEM_VERIFY(u_max_new < u_max + 1e-12, + MFEM_VERIFY(u_max_new < u_max + eps, "Overshoot of " << u_max_new - u_max); - MFEM_VERIFY(s_min_new > s_min - 1e-12, + MFEM_VERIFY(s_min_new > s_min - eps, "Undershoot in s of " << s_min - s_min_new); - MFEM_VERIFY(s_max_new < s_max + 1e-12, + MFEM_VERIFY(s_max_new < s_max + eps, "Overshoot in s of " << s_max_new - s_max); } u_min = u_min_new; @@ -1053,13 +1053,13 @@ int main(int argc, char *argv[]) { if (myid == 0) { - MFEM_VERIFY(u_min_new > 0.0 - 1e-12, + MFEM_VERIFY(u_min_new > 0.0 - eps, "Undershoot of " << 0.0 - u_min_new); - MFEM_VERIFY(u_max_new < 1.0 + 1e-12, + MFEM_VERIFY(u_max_new < 1.0 + eps, "Overshoot of " << u_max_new - 1.0); - MFEM_VERIFY(s_min_new > 0.0 - 1e-12, + MFEM_VERIFY(s_min_new > 0.0 - eps, "Undershoot in s of " << 0.0 - s_min_new); - MFEM_VERIFY(s_max_new < 1.0 + 1e-12, + MFEM_VERIFY(s_max_new < 1.0 + eps, "Overshoot in s of " << s_max_new - 1.0); } } diff --git a/remhos_fct.cpp b/remhos_fct.cpp index 46b021c..d479af4 100644 --- a/remhos_fct.cpp +++ b/remhos_fct.cpp @@ -115,32 +115,6 @@ void FCTSolver::CalcCompatibleLOProduct(const ParGridFunction &us, dof_id = k*ndofs + j; d_us_LO_new(dof_id) = (u_new(dof_id) * s_avg - us(dof_id)) / dt; } - -#ifdef REMHOS_FCT_PRODUCT_DEBUG - // Check the LO product solution. - double us_min, us_max; - for (int j = 0; j < ndofs; j++) - { - dof_id = k*ndofs + j; - if (active_dofs[dof_id] == false) { continue; } - - us_min = s_min_loc(j) * u_new(dof_id); - us_max = s_max_loc(j) * u_new(dof_id); - - if (s_avg * u_new(dof_id) + eps < us_min || - s_avg * u_new(dof_id) - eps > us_max) - { - std::cout << "---\ns_avg * u: " << k << " " - << us_min << " " - << s_avg * u_new(dof_id) << " " - << us_max << endl - << u_new(dof_id) << " " << s_avg << endl - << s_min_loc(j) << " " << s_max_loc(j) << "\n---\n"; - - MFEM_ABORT("s_avg * u not in bounds"); - } - } -#endif } } @@ -287,43 +261,48 @@ void FluxBasedFCT::CalcFCTProduct(const ParGridFunction &us, const Vector &m, dus_lo_fct = d_us; } -#ifdef REMHOS_FCT_PRODUCT_DEBUG - // Check the bounds of the final solution. - const double eps = 1e-12; - Vector us_new(d_us.Size()); - add(1.0, us, dt, d_us, us_new); - for (int k = 0; k < NE; k++) + if (verify_bounds) { - if (active_el[k] == false) { continue; } - - for (int j = 0; j < ndofs; j++) + // Check the bounds of the final solution. + const double eps = 1e-12; + Vector us_new(d_us.Size()); + add(1.0, us, dt, d_us, us_new); + for (int k = 0; k < NE; k++) { - dof_id = k*ndofs + j; - if (active_dofs[dof_id] == false) { continue; } + if (active_el[k] == false) { continue; } - double s = us_new(dof_id) / u_new(dof_id); - if (s + eps < s_min(dof_id) || - s - eps > s_max(dof_id)) + for (int j = 0; j < ndofs; j++) { - std::cout << "Final s " << j << " " << k << " " - << s_min(dof_id) << " " - << s << " " - << s_max(dof_id) << std::endl; - std::cout << "---\n"; - } + dof_id = k*ndofs + j; + if (active_dofs[dof_id] == false) { continue; } - if (us_new(dof_id) + eps < us_min(dof_id) || - us_new(dof_id) - eps > us_max(dof_id)) - { - std::cout << "Final us " << j << " " << k << " " - << us_min(dof_id) << " " - << us_new(dof_id) << " " - << us_max(dof_id) << std::endl; - std::cout << "---\n"; + double s = us_new(dof_id) / u_new(dof_id); + if (s + eps < s_min(dof_id) || + s - eps > s_max(dof_id)) + { + std::cout << "Final s " << j << " " << k << " " + << s_min(dof_id) << " " + << s << " " + << s_max(dof_id) << std::endl; + std::cout << "---\n"; + + MFEM_ABORT("s not in bounds after FCT."); + } + + if (us_new(dof_id) + eps < us_min(dof_id) || + us_new(dof_id) - eps > us_max(dof_id)) + { + std::cout << "Final us " << j << " " << k << " " + << us_min(dof_id) << " " + << us_new(dof_id) << " " + << us_max(dof_id) << std::endl; + std::cout << "---\n"; + + MFEM_ABORT("us not in bounds after FCT."); + } } } } -#endif } void FluxBasedFCT::ComputeFluxMatrix(const ParGridFunction &u, diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp index 3fad883..1d64774 100644 --- a/remhos_solvers.cpp +++ b/remhos_solvers.cpp @@ -212,10 +212,11 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) // Update mask with the HO update if (use_masks) { UpdateMask(x, dxs[i], dct, mask); } + // // Form the unlimited update for the stage. // Note that it converts eq. (2.16) in JLG's paper into an update using // the previous limited updates. - if (d_i[i] != 1.) + // { // for mask = 0, we get dxs (nothing happens). // the loop below won't change it -> Forward Euler. @@ -223,14 +224,14 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt) if (use_masks) { AddMasked(mask, d_i[i]-1., dxs[i], dxs[i]); } else { dxs[i] *= d_i[i]; } } + // Use all previous limited updates. for (int j = 0; j < i; j++) { - // Use all previous limited updates. if (use_masks) { AddMasked(mask, d_i[j], dxs[j], dxs[i]); } else { dxs[i].Add(d_i[j], dxs[j]); } } - // Limit the step + // Limit the step (always a Forward Euler step). f->LimitMult(x, dxs[i]); // Update the state From 0ae34d3ba17d214991e8d1f3b3f380593800ec8b Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 00:58:33 -0800 Subject: [PATCH 32/39] updated the autotest --- autotest/out_baseline.dat | 104 ++++++++++++++++++++------------------ autotest/test.sh | 20 ++++---- remhos.cpp | 2 +- remhos_fct.cpp | 13 ----- 4 files changed, 65 insertions(+), 74 deletions(-) diff --git a/autotest/out_baseline.dat b/autotest/out_baseline.dat index 8cdebe5..dfa08bc 100644 --- a/autotest/out_baseline.dat +++ b/autotest/out_baseline.dat @@ -2,215 +2,219 @@ --- Method -ho 1 -lo 2 -fct 2 - Remap pacman nonper-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 1 -lo 2 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 1 -lo 2 -fct 2 Final mass u: 0.08479546635 Max value u: 0.8262759545 - Remap bump nonper-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 1 -lo 2 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 1 -lo 2 -fct 2 Final mass u: 0.1197299711 Max value u: 0.9998930413 - Transport per-1D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 1 -lo 2 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 1 -lo 2 -fct 2 Final mass u: 0.1401241455 Max value u: 0.9052498201 - Transport bump per-unstruct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 1 -lo 2 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 1 -lo 2 -fct 2 Final mass u: 0.3888354875 Max value u: 0.9854644631 - Transport balls-jacks per-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 1 -lo 2 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 1 -lo 2 -fct 2 Final mass u: 0.1623263888 Max value u: 0.7742737139 - Transport bump per-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 1 -lo 2 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 1 -lo 2 -fct 2 Final mass u: 0.9607429525 Max value u: 0.9724537077 - Transport bump nonper-unstruct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 1 -lo 2 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 1 -lo 2 -fct 2 Final mass u: 0.8075105661 Max value u: 0.9999889304 --- Method -ho 3 -lo 4 -fct 2 - Remap pacman nonper-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 3 -lo 4 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 3 -lo 4 -fct 2 Final mass u: 0.0847954729 Max value u: 0.7581364675 - Remap bump nonper-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 3 -lo 4 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 3 -lo 4 -fct 2 Final mass u: 0.1197299801 Max value u: 0.9997499683 - Transport per-1D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 3 -lo 4 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 3 -lo 4 -fct 2 Final mass u: 0.1401241455 Max value u: 0.9039764015 - Transport bump per-unstruct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 3 -lo 4 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 3 -lo 4 -fct 2 Final mass u: 0.3888354875 Max value u: 0.9850024108 - Transport balls-jacks per-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 3 -lo 4 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 3 -lo 4 -fct 2 Final mass u: 0.1623263888 Max value u: 0.7145371968 - Transport bump per-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 3 -lo 4 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 3 -lo 4 -fct 2 Final mass u: 0.9607429525 Max value u: 0.9334903111 - Transport bump nonper-unstruct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 3 -lo 4 -fct 2 +mpirun -np 2 ./remhos -no-vis -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 3 -lo 4 -fct 2 Final mass u: 0.814998186 Max value u: 0.9999889315 --- Method -ho 2 -lo 3 -fct 2 -pa - Remap pacman nonper-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 2 -lo 3 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 2 -lo 3 -fct 2 -pa Final mass u: 0.08479546775 Max value u: 0.7779015453 - Remap bump nonper-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 2 -lo 3 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 2 -lo 3 -fct 2 -pa Final mass u: 0.1197300033 Max value u: 0.9997879406 - Transport per-1D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 2 -lo 3 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 2 -lo 3 -fct 2 -pa Final mass u: 0.1401241455 Max value u: 0.8781060808 - Transport bump per-unstruct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 2 -lo 3 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 2 -lo 3 -fct 2 -pa Final mass u: 0.3888354875 Max value u: 0.9755502191 - Transport balls-jacks per-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 2 -lo 3 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 2 -lo 3 -fct 2 -pa Final mass u: 0.1623263888 Max value u: 0.6374820899 - Transport bump per-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 2 -lo 3 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 2 -lo 3 -fct 2 -pa Final mass u: 0.9607429525 Max value u: 0.9202929163 - Transport bump nonper-unstruct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 2 -lo 3 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 2 -lo 3 -fct 2 -pa Final mass u: 0.7772459527 Max value u: 0.9999889307 --- Method -ho 2 -lo 4 -fct 2 -pa - Remap pacman nonper-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 2 -lo 4 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 2 -lo 4 -fct 2 -pa Final mass u: 0.0847954729 Max value u: 0.7581364676 - Remap bump nonper-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 2 -lo 4 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 2 -lo 4 -fct 2 -pa Final mass u: 0.1197299801 Max value u: 0.9997499683 - Transport per-1D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 2 -lo 4 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 2 -lo 4 -fct 2 -pa Final mass u: 0.1401241455 Max value u: 0.9039764015 - Transport bump per-unstruct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 2 -lo 4 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 2 -lo 4 -fct 2 -pa Final mass u: 0.3888354875 Max value u: 0.9850024108 - Transport balls-jacks per-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 2 -lo 4 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 2 -lo 4 -fct 2 -pa Final mass u: 0.1623263888 Max value u: 0.7145371968 - Transport bump per-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 2 -lo 4 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 2 -lo 4 -fct 2 -pa Final mass u: 0.9607429525 Max value u: 0.9334903111 - Transport bump nonper-unstruct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 2 -lo 4 -fct 2 -pa +mpirun -np 2 ./remhos -no-vis -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 2 -lo 4 -fct 2 -pa Final mass u: 0.7779917929 Max value u: 0.9999889315 --- Method -ho 3 -lo 1 -fct 1 - Remap pacman nonper-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 3 -lo 1 -fct 1 +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.0015 -tf 0.75 -ho 3 -lo 1 -fct 1 Final mass u: 0.08479546845 Max value u: 0.905654904 - Remap bump nonper-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 3 -lo 1 -fct 1 +mpirun -np 2 ./remhos -no-vis -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.7 -ho 3 -lo 1 -fct 1 Final mass u: 0.11972981 Max value u: 0.996173945 - Transport per-1D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 3 -lo 1 -fct 1 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-segment.mesh -p 0 -rs 3 -dt 0.001 -tf 1 -ho 3 -lo 1 -fct 1 Final mass u: 0.1401241455 Max value u: 0.9071157249 - Transport bump per-unstruct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 3 -lo 1 -fct 1 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-hexagon.mesh -p 0 -rs 2 -dt 0.005 -tf 2.5 -ho 3 -lo 1 -fct 1 Final mass u: 0.3888354875 Max value u: 0.9979069772 - Transport balls-jacks per-struct-2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 3 -lo 1 -fct 1 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.004 -tf 0.8 -ho 3 -lo 1 -fct 1 Final mass u: 0.1623263888 Max value u: 0.787875182 - Transport bump per-struct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 3 -lo 1 -fct 1 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.015 -tf 2 -ho 3 -lo 1 -fct 1 Final mass u: 0.9607429525 Max value u: 0.9984668427 - Transport bump nonper-unstruct-3D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 3 -lo 1 -fct 1 +mpirun -np 2 ./remhos -no-vis -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.035 -tf 3 -ho 3 -lo 1 -fct 1 Final mass u: 0.8143001155 Max value u: 0.9995762608 --- Product remap 2D (FCT) -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 13 -Final mass us: 0.1767475452 -Mass loss us: 0.00286033 +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 1 +Final mass us: 0.1815368098 +Mass loss us: 0.00192894 ---- Product remap 2D (ClipScale) -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 13 -Final mass us: 0.1782170448 -Mass loss us: 0.00139083 +--- Product remap 2D IDP2 (ClipScale) +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 1 -lo 5 -fct 2 -ps -s 12 +Final mass us: 0.1796076412 +Mass loss us: 2.31348e-07 ---- Product remap 2D (FCTProject) -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 13 -Final mass us: 0.1776386793 -Mass loss us: 0.00196919 +--- Product remap 2D IDP3 (FCTProject) +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 5 -fct 4 -ps -s 13 +Final mass us: 0.179607829 +Mass loss us: 4.35885e-08 --- BLAST sharpening test - Pacman remap auto-dt -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt -1 -tf 0.75 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1 +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 1 -dt -1 -tf 0.75 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1 +Final mass u: 0.08479612805 +Mass loss u: 6.61247e-07 --- BLAST sharpening test - Transport balls-jacks auto-dt -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/periodic-square.mesh -p 5 -rs 3 -dt -1 -tf 0.8 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1 +mpirun -np 2 ./remhos -no-vis -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.01 -tf 0.8 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1 +Final mass u: 0.1623263888 +Mass loss u: 1.38778e-15 --- Steady monolithic 2 2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 7 -rs 3 -o 1 -dt 0.01 -tf 20 -mono 1 -si 2 +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 7 -rs 3 -o 1 -dt 0.01 -tf 20 -mono 1 -si 2 Final mass u: 0.1570667907 Max value u: 0.9987771164 --- Steady monolithic 1 2D -mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 6 -rs 2 -o 1 -dt 0.01 -tf 20 -mono 1 -si 1 +mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 6 -rs 2 -o 1 -dt 0.01 -tf 20 -mono 1 -si 1 Final mass u: 0.3182739921 Max value u: 1 diff --git a/autotest/test.sh b/autotest/test.sh index 79029fc..b692b5a 100755 --- a/autotest/test.sh +++ b/autotest/test.sh @@ -9,9 +9,9 @@ file="autotest/out_test.dat" ntask=$1 if [ "$2" = "cuda" ]; then - command="lrun -n "$((ntask))" ./remhos -no-vis --verify-bounds -d cuda" + command="lrun -n "$((ntask))" ./remhos -no-vis -d cuda" else - command="mpirun -np "$((ntask))" ./remhos -no-vis --verify-bounds" + command="mpirun -np "$((ntask))" ./remhos -no-vis" fi methods=( "-ho 1 -lo 2 -fct 2" # Hennes 1 @@ -66,29 +66,29 @@ for method in "${methods[@]}"; do done echo -e '\n'"--- Product remap 2D (FCT)" >> $file -run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 13" +run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 1" echo -e $run_line >> $file $run_line | grep -e 'mass us' -e 'loss us'>> $file -echo -e '\n'"--- Product remap 2D (ClipScale)" >> $file -run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 13" +echo -e '\n'"--- Product remap 2D IDP2 (ClipScale)" >> $file +run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 1 -lo 5 -fct 2 -ps -s 12" echo -e $run_line >> $file $run_line | grep -e 'mass us' -e 'loss us'>> $file -echo -e '\n'"--- Product remap 2D (FCTProject)" >> $file -run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 13" +echo -e '\n'"--- Product remap 2D IDP3 (FCTProject)" >> $file +run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 5 -fct 4 -ps -s 13" echo -e $run_line >> $file $run_line | grep -e 'mass us' -e 'loss us'>> $file echo -e '\n'"--- BLAST sharpening test - Pacman remap auto-dt" >> $file run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 1 -dt -1 -tf 0.75 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1" echo -e $run_line >> $file -$run_line | grep -e 'mass us' -e 'loss us'>> $file +$run_line | grep -e 'mass u' -e 'loss u'>> $file echo -e '\n'"--- BLAST sharpening test - Transport balls-jacks auto-dt" >> $file -run_line=$command" -m ./data/periodic-square.mesh -p 5 -rs 3 -dt -1 -tf 0.8 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1" +run_line=$command" -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.01 -tf 0.8 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1" echo -e $run_line >> $file -$run_line | grep -e 'mass us' -e 'loss us'>> $file +$run_line | grep -e 'mass u' -e 'loss u'>> $file echo -e '\n'"--- Steady monolithic 2 2D" >> $file run_line=$command" -m ./data/inline-quad.mesh -p 7 -rs 3 -o 1 -dt 0.01 -tf 20 -mono 1 -si 2" diff --git a/remhos.cpp b/remhos.cpp index 11c5cd0..f4b824a 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -925,7 +925,7 @@ int main(int argc, char *argv[]) ho_solver, lo_solver, fct_solver, mono_solver); adv.verify_bounds = verify_bounds; - fct_solver->verify_bounds = verify_bounds; + if (fct_solver) { fct_solver->verify_bounds = verify_bounds; } double t = 0.0; adv.SetTime(t); diff --git a/remhos_fct.cpp b/remhos_fct.cpp index d479af4..53dc595 100644 --- a/remhos_fct.cpp +++ b/remhos_fct.cpp @@ -276,19 +276,6 @@ void FluxBasedFCT::CalcFCTProduct(const ParGridFunction &us, const Vector &m, dof_id = k*ndofs + j; if (active_dofs[dof_id] == false) { continue; } - double s = us_new(dof_id) / u_new(dof_id); - if (s + eps < s_min(dof_id) || - s - eps > s_max(dof_id)) - { - std::cout << "Final s " << j << " " << k << " " - << s_min(dof_id) << " " - << s << " " - << s_max(dof_id) << std::endl; - std::cout << "---\n"; - - MFEM_ABORT("s not in bounds after FCT."); - } - if (us_new(dof_id) + eps < us_min(dof_id) || us_new(dof_id) - eps > us_max(dof_id)) { From 2d67a960a3a45a260d28deadbf15703637739d30 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 09:54:14 -0800 Subject: [PATCH 33/39] less sensitivity to small numbers in autotest --- autotest/test.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/autotest/test.sh b/autotest/test.sh index b692b5a..fee9f16 100755 --- a/autotest/test.sh +++ b/autotest/test.sh @@ -78,7 +78,7 @@ $run_line | grep -e 'mass us' -e 'loss us'>> $file echo -e '\n'"--- Product remap 2D IDP3 (FCTProject)" >> $file run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 5 -fct 4 -ps -s 13" echo -e $run_line >> $file -$run_line | grep -e 'mass us' -e 'loss us'>> $file +$run_line | grep -e 'mass u' -e 'mass us'>> $file echo -e '\n'"--- BLAST sharpening test - Pacman remap auto-dt" >> $file run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 1 -dt -1 -tf 0.75 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1" @@ -88,7 +88,7 @@ $run_line | grep -e 'mass u' -e 'loss u'>> $file echo -e '\n'"--- BLAST sharpening test - Transport balls-jacks auto-dt" >> $file run_line=$command" -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.01 -tf 0.8 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1" echo -e $run_line >> $file -$run_line | grep -e 'mass u' -e 'loss u'>> $file +$run_line | grep -e 'mass u' -e 'value u'>> $file echo -e '\n'"--- Steady monolithic 2 2D" >> $file run_line=$command" -m ./data/inline-quad.mesh -p 7 -rs 3 -o 1 -dt 0.01 -tf 20 -mono 1 -si 2" From 3f49ec77c642247fd6136cebef4c43eeae68f59f Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 10:20:18 -0800 Subject: [PATCH 34/39] forgot to update .dat --- autotest/out_baseline.dat | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/autotest/out_baseline.dat b/autotest/out_baseline.dat index dfa08bc..836e9a8 100644 --- a/autotest/out_baseline.dat +++ b/autotest/out_baseline.dat @@ -196,8 +196,8 @@ Mass loss us: 2.31348e-07 --- Product remap 2D IDP3 (FCTProject) mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 5 -fct 4 -ps -s 13 +Final mass u: 0.08980386855 Final mass us: 0.179607829 -Mass loss us: 4.35885e-08 --- BLAST sharpening test - Pacman remap auto-dt mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 14 -rs 1 -dt -1 -tf 0.75 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1 @@ -207,7 +207,7 @@ Mass loss u: 6.61247e-07 --- BLAST sharpening test - Transport balls-jacks auto-dt mpirun -np 2 ./remhos -no-vis -m ./data/periodic-square.mesh -p 5 -rs 3 -dt 0.01 -tf 0.8 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1 Final mass u: 0.1623263888 -Mass loss u: 1.38778e-15 +Max value u: 0.2863317261 --- Steady monolithic 2 2D mpirun -np 2 ./remhos -no-vis -m ./data/inline-quad.mesh -p 7 -rs 3 -o 1 -dt 0.01 -tf 20 -mono 1 -si 2 From 5297ef1b10527f398982afd5753617c2688a06ca Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 12:57:55 -0800 Subject: [PATCH 35/39] update workflow versions --- .github/workflows/build-and-test.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index e3c68aa..7022608 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -35,7 +35,7 @@ jobs: # Install will only run on cache miss. - name: cache hypre id: hypre-cache - uses: actions/cache@v2 + uses: actions/cache@v4 with: path: ${{ env.HYPRE_TOP_DIR }} key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2 @@ -51,7 +51,7 @@ jobs: # Install will only run on cache miss. - name: cache metis id: metis-cache - uses: actions/cache@v2 + uses: actions/cache@v4 with: path: ${{ env.METIS_TOP_DIR }} key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2 @@ -80,7 +80,7 @@ jobs: # Install will only run on cache miss. - name: cache mfem id: mfem-cache - uses: actions/cache@v2 + uses: actions/cache@v4 with: path: ${{ env.MFEM_TOP_DIR }} key: ${{ runner.os }}-build-${{ env.MFEM_TOP_DIR }}-${{ env.MFEM_COMMIT }}-v3 @@ -115,7 +115,7 @@ jobs: - name: Archive test results patch if: always() - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: baseline-patch path: remhos/baseline.patch From 698113df828e47d681cbf990433b984f8da35219 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 17:07:44 -0800 Subject: [PATCH 36/39] reversing workflow runnner image --- .github/workflows/build-and-test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index 7022608..29e7932 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -16,7 +16,7 @@ env: jobs: build-and-test: - runs-on: ubuntu-latest + runs-on: ubuntu-22.04 steps: # Checkout Remhos in "remhos" subdirectory. Final path: /home/runner/work/remhos/remhos/remhos From 6a6ee57543efcf122b6595f240575ab639e23729 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 17:13:39 -0800 Subject: [PATCH 37/39] build library only in workflows --- .github/workflows/build-and-test.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index 29e7932..d660853 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -94,6 +94,7 @@ jobs: metis-dir: ${{ env.METIS_TOP_DIR }} mfem-dir: ${{ env.MFEM_TOP_DIR }} mfem-branch: ${{ env.MFEM_BRANCH }} + library-only: 'true' - name: build Remhos run: | From 5180847aa8fdebc8b1ff8e8a3d09cea20693aa69 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 17:42:37 -0800 Subject: [PATCH 38/39] minor --- README.md | 4 ++-- remhos.cpp | 6 ------ 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 0d97aa4..dcee039 100644 --- a/README.md +++ b/README.md @@ -233,7 +233,7 @@ Alternatively, verify the final mass (`mass`) and maximum value (`max`) for the 7. `mpirun -np 8 remhos -m ./data/periodic-cube.mesh -p 0 -rs 1 -o 2 -dt 0.014 -tf 8 -ho 1 -lo 4 -fct 2` 8. `mpirun -np 8 remhos -m ../mfem/data/ball-nurbs.mesh -p 1 -rs 1 -dt 0.02 -tf 3 -ho 1 -lo 4 -fct 2` 9. `mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 14 -rs 1 -dt 0.001 -tf 0.75 -ho 1 -lo 4 -fct 2` -10. `mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps` +10. `mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 14 -rs 3 -dt 0.005 -tf 0.75 -ho 1 -lo 5 -fct 4 -ps -s 13` 11. `mpirun -np 8 remhos -m ./data/cube01_hex.mesh -p 10 -rs 1 -o 2 -dt 0.02 -tf 0.8 -ho 1 -lo 4 -fct 2` 12. `mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 7 -rs 3 -o 1 -dt 0.01 -tf 20 -mono 1 -si 2` 13. `mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 6 -rs 2 -o 1 -dt 0.01 -tf 20 -mono 1 -si 1` @@ -250,7 +250,7 @@ Alternatively, verify the final mass (`mass`) and maximum value (`max`) for the | 7. | 0.9607429525 | 0.7678305756 | | 8. | 0.8087104604 | 0.9999889315 | | 9. | 0.08479546709| 0.8156091428 | -| 10. | 0.08980397023| 0.9886734209 | +| 10. | 0.09317738757| 0.9994170644 | | 11. | 0.1197294512 | 0.9990312449 | | 12. | 0.1570667907 | 0.9987771164 | | 13. | 0.3182739921 | 1 | diff --git a/remhos.cpp b/remhos.cpp index f4b824a..df487cb 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -26,12 +26,6 @@ // equations that are used to perform discontinuous field interpolation (remap) // as part of the Eulerian phase in Arbitrary-Lagrangian Eulerian (ALE) // simulations. -// -// single step: -// mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 14 -rs 3 -dt 0.002 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 1 -vs 20 -// RK2: -// mpirun -np 8 remhos -m ./data/inline-quad.mesh -p 14 -rs 3 -dt 0.002 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 2 -vs 20 -// #include "mfem.hpp" #include From 37692f967a3abb00d5d9b67385479ca372c8b9ca Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Thu, 19 Dec 2024 17:44:36 -0800 Subject: [PATCH 39/39] minor --- remhos.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/remhos.cpp b/remhos.cpp index df487cb..19dd616 100644 --- a/remhos.cpp +++ b/remhos.cpp @@ -26,6 +26,8 @@ // equations that are used to perform discontinuous field interpolation (remap) // as part of the Eulerian phase in Arbitrary-Lagrangian Eulerian (ALE) // simulations. +// +// Sample runs: see README.md, section 'Verification of Results'. #include "mfem.hpp" #include