From 2fcdb252792f9565e26b2f3bdc1d0ca324ae7bb5 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 11 Dec 2025 10:05:40 -0800 Subject: [PATCH 01/28] missing flag check for rebalancing --- laghos.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/laghos.cpp b/laghos.cpp index 0996fa22..2d87340f 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -443,7 +443,7 @@ int main(int argc, char *argv[]) // Refine the mesh further in parallel to increase the resolution. for (int lev = 0; lev < rp_levels; lev++) { pmesh->UniformRefinement(); } - if (!cartesian_partitioning && enable_nc && dim > 1) + if (!cartesian_partitioning && enable_nc && dim > 1 && enable_rebalance) { if (myid == 0) { cout << "Rebalancing mesh" << endl; } pmesh->Rebalance(); From 18c064be2ad8144adf6d67bf2402acb9c249c4cd Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 18 Dec 2025 14:45:25 -0800 Subject: [PATCH 02/28] Added ability to generate the triple-pt mesh --- laghos.cpp | 86 ++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 70 insertions(+), 16 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index 1477a35d..a323c225 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -87,6 +87,10 @@ double rho0(const Vector &); double gamma_func(const Vector &); void v0(const Vector &, Vector &); +// for generated meshes +static void AssignMeshBdrAttrs2D(Mesh &, real_t, real_t); +static void AssignMeshBdrAttrs3D(Mesh &, real_t, real_t, real_t, real_t); + static long GetMaxRssMB(); static void display_banner(std::ostream&); static void Checks(const int ti, const double norm, int &checks); @@ -285,26 +289,32 @@ int main(int argc, char *argv[]) } if (dim == 2) { - mesh = new Mesh(Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, - true)); - const int NBE = mesh->GetNBE(); - for (int b = 0; b < NBE; b++) - { - Element *bel = mesh->GetBdrElement(b); - const int attr = (b < NBE/2) ? 2 : 1; - bel->SetAttribute(attr); + switch (problem) { + case 3: + mesh = new Mesh(Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, + true, 7_r, 3_r)); + AssignMeshBdrAttrs2D(*mesh, 0_r, 7_r); + break; + default: + mesh = new Mesh( + Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true)); + AssignMeshBdrAttrs2D(*mesh, 0_r, 1_r); + break; } } if (dim == 3) { - mesh = new Mesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, - true)); - const int NBE = mesh->GetNBE(); - for (int b = 0; b < NBE; b++) - { - Element *bel = mesh->GetBdrElement(b); - const int attr = (b < NBE/3) ? 3 : (b < 2*NBE/3) ? 1 : 2; - bel->SetAttribute(attr); + switch (problem) { + case 3: + mesh = new Mesh(Mesh::MakeCartesian3D( + nx, ny, nz, Element::HEXAHEDRON, 7_r, 3_r, 3_r, true)); + AssignMeshBdrAttrs3D(*mesh, 0_r, 7_r, 0_r, 3_r); + break; + default: + mesh = new Mesh(Mesh::MakeCartesian3D( + nx, ny, nz, Element::HEXAHEDRON, 1_r, 1_r, 1_r, true)); + AssignMeshBdrAttrs3D(*mesh, 0_r, 1_r, 0_r, 1_r); + break; } } } @@ -1213,3 +1223,47 @@ static void Checks(const int ti, const double nrm, int &chk) } } } + +static void AssignMeshBdrAttrs2D(Mesh& mesh, real_t xmin, real_t xmax) +{ + Vector pos(3); + constexpr real_t tol = 1e-6; + const int NBE = mesh.GetNBE(); + IntegrationPoint center; + center.x = 0.5; + center.y = 0.5; + center.z = 0.5; + for (int b = 0; b < NBE; b++) { + Element *bel = mesh.GetBdrElement(b); + auto eltrans = mesh.GetBdrElementTransformation(b); + eltrans->Transform(center, pos); + int attr = 2; + if (pos[0] <= xmin + tol || pos[0] >= xmax - tol) { + attr = 1; + } + bel->SetAttribute(attr); + } +} + +static void AssignMeshBdrAttrs3D(Mesh &mesh, real_t xmin, real_t xmax, + real_t ymin, real_t ymax) { + Vector pos(3); + constexpr real_t tol = 1e-6; + const int NBE = mesh.GetNBE(); + IntegrationPoint center; + center.x = 0.5; + center.y = 0.5; + center.z = 0.5; + for (int b = 0; b < NBE; b++) { + Element *bel = mesh.GetBdrElement(b); + auto eltrans = mesh.GetBdrElementTransformation(b); + eltrans->Transform(center, pos); + int attr = 3; + if (pos[0] <= xmin + tol || pos[0] >= xmax - tol) { + attr = 1; + } else if (pos[1] <= ymin + tol || pos[1] >= ymax - tol) { + attr = 2; + } + bel->SetAttribute(attr); + } +} From b1375b842193444a1229d9e5f070c42b2aacae70 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Mon, 12 Jan 2026 09:22:55 -0800 Subject: [PATCH 03/28] Added a way to compute the exact sedov shock solution --- CMakeLists.txt | 16 ++- adaptive_quad.hpp | 160 ++++++++++++++++++++++++ bisect.hpp | 76 ++++++++++++ laghos.cpp | 84 ++++++++++++- sedov/sedov.cpp | 306 ++++++++++++++++++++++++++++++++++++++++++++++ sedov_sol.cpp | 173 ++++++++++++++++++++++++++ sedov_sol.hpp | 77 ++++++++++++ 7 files changed, 888 insertions(+), 4 deletions(-) create mode 100644 adaptive_quad.hpp create mode 100644 bisect.hpp create mode 100644 sedov/sedov.cpp create mode 100644 sedov_sol.cpp create mode 100644 sedov_sol.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 917dd2a3..a77949da 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,7 +53,7 @@ endif() list(APPEND SOURCES - laghos_assembly.cpp laghos.cpp laghos_solver.cpp + laghos_assembly.cpp laghos.cpp laghos_solver.cpp sedov_sol.cpp ) if (MFEM_USE_CUDA) @@ -69,3 +69,17 @@ target_link_libraries(laghos PUBLIC MPI::MPI_CXX PUBLIC HYPRE::HYPRE ) + +if (MFEM_USE_CUDA) + set_source_files_properties(sedov_sol.cpp sedov/sedov.cpp PROPERTIES LANGUAGE CUDA) +endif() + +add_executable(sedov sedov_sol.cpp sedov/sedov.cpp) + +target_include_directories(sedov PUBLIC "${MFEM_DIR}/../../../include/mfem") + +target_link_libraries(sedov + PUBLIC mfem + PUBLIC MPI::MPI_CXX + PUBLIC HYPRE::HYPRE +) diff --git a/adaptive_quad.hpp b/adaptive_quad.hpp new file mode 100644 index 00000000..d3a27bf1 --- /dev/null +++ b/adaptive_quad.hpp @@ -0,0 +1,160 @@ +// 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. + +#ifndef LAGHOS_ADAPTIVE_QUAD_HPP +#define LAGHOS_ADAPTIVE_QUAD_HPP + +#include +#include + +/// +/// Implements the 21-point adaptive Gauss-Kronrod quadrature method +/// +template struct gk21 { + Fun fun; + Err err_fun; + using res_type = decltype(fun(0.0)); + constexpr static size_t gl_points() { return 10; } + constexpr static size_t gk_points() { return 11; } + +private: + res_type integrate_recurse(double lower, double upper, size_t curr_depth, + size_t max_depth = 20) const { + static constexpr double data[] = { + // gl_abscissa + -1.488743389816312108848260011297200e-01, + -4.333953941292471907992659431657842e-01, + -6.794095682990244062343273651148736e-01, + -8.650633666889845107320966884234930e-01, + -9.739065285171717200779640120844521e-01, + 1.488743389816312108848260011297200e-01, + 4.333953941292471907992659431657842e-01, + 6.794095682990244062343273651148736e-01, + 8.650633666889845107320966884234930e-01, + 9.739065285171717200779640120844521e-01, + // gl_weights + 2.955242247147528701738929946513383e-01, + 2.692667193099963550912269215694694e-01, + 2.190863625159820439955349342281632e-01, + 1.494513491505805931457763396576973e-01, + 6.667134430868813759356880989333179e-02, + 2.955242247147528701738929946513383e-01, + 2.692667193099963550912269215694694e-01, + 2.190863625159820439955349342281632e-01, + 1.494513491505805931457763396576973e-01, + 6.667134430868813759356880989333179e-02, + // glk_weights + 1.477391049013384913748415159720680e-01, + 1.347092173114733259280540017717068e-01, + 1.093871588022976418992105903258050e-01, + 7.503967481091995276704314091619001e-02, + 3.255816230796472747881897245938976e-02, + 1.477391049013384913748415159720680e-01, + 1.347092173114733259280540017717068e-01, + 1.093871588022976418992105903258050e-01, + 7.503967481091995276704314091619001e-02, + 3.255816230796472747881897245938976e-02, + // gk_abscissa + 0.000000000000000000000000000000000e00, + -2.943928627014601981311266031038656e-01, + -5.627571346686046833390000992726941e-01, + -7.808177265864168970637175783450424e-01, + -9.301574913557082260012071800595083e-01, + -9.956571630258080807355272806890028e-01, + 2.943928627014601981311266031038656e-01, + 5.627571346686046833390000992726941e-01, + 7.808177265864168970637175783450424e-01, + 9.301574913557082260012071800595083e-01, + 9.956571630258080807355272806890028e-01, + // gk_weights + 1.494455540029169056649364683898212e-01, + 1.427759385770600807970942731387171e-01, + 1.234919762620658510779581098310742e-01, + 9.312545458369760553506546508336634e-02, + 5.475589657435199603138130024458018e-02, + 1.169463886737187427806439606219205e-02, + 1.427759385770600807970942731387171e-01, + 1.234919762620658510779581098310742e-01, + 9.312545458369760553506546508336634e-02, + 5.475589657435199603138130024458018e-02, + 1.169463886737187427806439606219205e-02, + }; + // TODO: where to copy gk21_base data to scratch memory? + res_type gl_sum = 0; + res_type gk_sum = 0; + double jac = (upper - lower) * 0.5; + for (int i = 0; i < gl_points(); ++i) { + res_type f_eval = fun((data[i] + 1) * jac + lower); + gl_sum += f_eval * data[gl_points() + i]; + gk_sum += f_eval * data[2 * gl_points() + i]; + } + for (int i = 0; i < gk_points(); ++i) { + res_type f_eval = fun((data[3 * gl_points() + i] + 1) * jac + lower); + gk_sum += f_eval * data[3 * gl_points() + gk_points() + i]; + } + gk_sum *= jac; + gl_sum *= jac; + if (curr_depth < max_depth && !err_fun(gk_sum, gl_sum)) { + gk_sum = integrate_recurse(lower, lower + jac, curr_depth + 1, max_depth); + gk_sum += + integrate_recurse(lower + jac, upper, curr_depth + 1, max_depth); + } + + return gk_sum; + } + +public: + res_type integrate(double lower, double upper, size_t start_segs = 1, + size_t max_depth = 20) const { + double dx = (upper - lower) / start_segs; + res_type res = 0; + double curr = lower; + for (size_t i = 0; i < start_segs; ++i) { + double next = lower + (i + 1) * dx; + res += integrate_recurse(curr, next, 1, max_depth); + curr = next; + } + return res; + } +}; + +template +auto gk21_integrate(F &&f, E &&e, double lower, double upper, + size_t start_segs = 1, size_t max_depth = 20) { + gk21 integrator{f, e}; + return integrator.integrate(lower, upper, start_segs, max_depth); +} + +struct scalar_error_functor { + double eps_abs; + double eps_rel; + template bool operator()(const T &ho, const T &lo) const { + if (!std::isfinite(ho)) { + return true; + } + double delta = std::fabs(ho - lo); + if (delta < eps_abs) { + return true; + } + double denom = std::max(std::fabs(ho), std::fabs(lo)); + if (delta < eps_rel * denom) { + return true; + } + return false; + } +}; + +#endif diff --git a/bisect.hpp b/bisect.hpp new file mode 100644 index 00000000..2259d9c4 --- /dev/null +++ b/bisect.hpp @@ -0,0 +1,76 @@ +// 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. + +#ifndef LAGHOS_BISECT_HPP +#define LAGHOS_BISECT_HPP + +#include +#include + +#include + +/// Bisection root finder +template double bisection(Fun &&fun, double lower, double upper) { + double lv = fun(lower); + constexpr double tol = 1e-20; + if (std::fabs(lv) < tol) { + return lower; + } + double rv = fun(upper); + if (std::fabs(rv) < tol) { + return upper; + } + if (std::copysign(1., lv) * std::copysign(1., rv) > 0) { + throw std::runtime_error("bisection: no sign change"); + } + auto dx_init = upper - lower; + auto dx_last = dx_init; + while (true) { + double mid = 0.5 * (lower + upper); + auto dx = mid - lower; + double mv = fun(mid); + if (dx < dx_init * 1e-16 || dx >= dx_last) { + if (fabs(mv) < fabs(lv)) { + if (fabs(mv) < fabs(rv)) { + return mid; + } else if (fabs(rv) < fabs(lv)) { + return upper; + } else { + return lower; + } + } else if (fabs(rv) < fabs(lv)) { + return upper; + } else { + return lower; + } + } + if (std::fabs(mv) < tol) { + return mid; + } + if (std::copysign(1., lv) != std::copysign(1., mv)) { + upper = mid; + rv = mv; + } else if (std::copysign(1., rv) != std::copysign(1., mv)) { + lower = mid; + lv = mv; + } else { + throw std::runtime_error("bisection: no sign change"); + } + dx_last = dx; + } +} + +#endif diff --git a/laghos.cpp b/laghos.cpp index a323c225..0da0c719 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -63,6 +63,7 @@ #include #include #include "laghos_solver.hpp" +#include "sedov_sol.hpp" #ifdef USE_CALIPER #include #include @@ -135,13 +136,12 @@ int main(int argc, char *argv[]) int partition_type = 0; const char *device = "cpu"; bool check = false; + bool check_exact = false; bool mem_usage = false; bool fom = false; bool gpu_aware_mpi = false; int dev = 0; int dev_pool_size = 4; - double blast_energy = 0.25; - double blast_position[] = {0.0, 0.0, 0.0}; bool enable_nc = true; bool enable_rebalance = true; @@ -213,6 +213,10 @@ int main(int argc, char *argv[]) "Device configuration string, see Device::Configure()."); args.AddOption(&check, "-chk", "--checks", "-no-chk", "--no-checks", "Enable 2D checks."); + args.AddOption(&check_exact, "-err", "--exact-error", "-no-err", + "--no-exact-error", + "Enable comparing the Sedov problem (problem 1) against the " + "exact solution."); args.AddOption(&mem_usage, "-mb", "--mem", "-no-mem", "--no-mem", "Enable memory usage."); args.AddOption(&fom, "-f", "--fom", "-no-fom", "--no-fom", @@ -237,6 +241,13 @@ int main(int argc, char *argv[]) } if (Mpi::Root()) { args.PrintOptions(cout); } + if (check_exact) { + MFEM_VERIFY( + problem == 1, + "Can only compare problem 1 (Sedov) against the exact solution"); + MFEM_VERIFY(strncmp(mesh_file, "default", 7) == 0, "check: mesh_file"); + } + #ifdef USE_CALIPER cali_config_set("CALI_CALIPER_ATTRIBUTE_DEFAULT_SCOPE", "process"); CALI_CXX_MARK_FUNCTION; @@ -598,11 +609,16 @@ int main(int argc, char *argv[]) ParGridFunction l2_rho0_gf(&l2_fes), l2_e(&l2_fes); l2_rho0_gf.ProjectCoefficient(rho0_coeff); rho0_gf.ProjectGridFunction(l2_rho0_gf); + + double blast_energy = 1; + double blast_position[] = {0.0, 0.0, 0.0}; if (problem == 1) { // For the Sedov test, we use a delta function at the origin. + // divide amount of blast energy by 2^d due to simulating only a portion + // of the symmetric blast. DeltaCoefficient e_coeff(blast_position[0], blast_position[1], - blast_position[2], blast_energy); + blast_position[2], blast_energy / pow(2, dim)); l2_e.ProjectCoefficient(e_coeff); } else @@ -958,6 +974,68 @@ int main(int argc, char *argv[]) adiak::fini(); #endif + if (check_exact) { + // compare against the exact Sedov solution + double gamma = 1.4; + double rho0 = 1; + double omega = 0; + + SedovSol asol(dim, gamma, rho0, blast_energy, omega); + + asol.SetTime(t_final); + + MFEM_VERIFY(asol.r2 <= 0.9, + "Solution reflections off boundaries detected, cannot compare " + "against exact solution."); + + int err_order = std::max((std::max(order_v, order_e) + 1) * 2, order_q) * 2; + const IntegrationRule &irule = + IntRules.Get(pmesh->GetTypicalElementGeometry(), err_order); + + QuadratureSpace qspace(*pmesh, irule); + // only compare density + QuadratureFunction sim_qfunc(qspace, 1); + QuadratureFunction err_qfunc(qspace, 1); + + hydro.ComputeDensity(rho_gf); + + sim_qfunc.ProjectGridFunction(rho_gf); + + auto slambda = [&](const Vector &x, Vector &res) { + real_t tmp[3]; + Vector dr(tmp, dim); + double r = 0; + + for (int i = 0; i < dim; ++i) { + dr[i] = x[i] - blast_position[i]; + r += dr[i] * dr[i]; + } + r = sqrt(r); + if (r) { + for (int i = 0; i < dim; ++i) { + dr[i] /= r; + } + } else { + dr = 0_r; + } + double rho, v, P; + asol.EvalSol(r, rho, v, P); + res[0] = rho; + }; + VectorFunctionCoefficient asol_coeff(1, slambda); + asol_coeff.Project(err_qfunc); + err_qfunc -= sim_qfunc; + + err_qfunc.HostReadWrite(); + for (int i = 0; i < err_qfunc.Size(); ++i) { + err_qfunc[i] = pow(err_qfunc[i], 2); + } + real_t lrho_err = err_qfunc.Integrate(); + if (Mpi::Root()) { + cout << "Density L2 error: " << sqrt(lrho_err) << endl; + } + } + // Free the used memory. delete ode_solver; delete pmesh; diff --git a/sedov/sedov.cpp b/sedov/sedov.cpp new file mode 100644 index 00000000..a93b58a2 --- /dev/null +++ b/sedov/sedov.cpp @@ -0,0 +1,306 @@ +// 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. + +// Computes and writes out the Sedov shock solution +// +// See: +// James R. Kamm, Evaluation of the Sedov-von Neumann-Taylor Blast Wave Solution +// LA-UR-00-6055 + +#include "../sedov_sol.hpp" + +#include + +#include +#include + +using namespace mfem; + +// static void ProjectCoeff(ParGridFunction &u, VectorCoefficient &coeff, +// const mfem::IntegrationRule *ir); + +int main(int argc, char *argv[]) { + // Initialize MPI. + Mpi::Init(); + int myid = Mpi::WorldRank(); + Hypre::Init(); + + // Parse command-line options. + int dim = 3; + int rs_levels = 2; + int rp_levels = 0; + int nx = 2; + int ny = 2; + int nz = 2; + int order_v = 2; + int order_e = 1; + int order_q = -1; + double t_final = 0.6; + const char *basename = "results/Sedov"; + + OptionsParser args(argc, argv); + args.AddOption(&dim, "-dim", "--dimension", "Dimension of the problem."); + args.AddOption(&nx, "-nx", "--xelems", "Elements in x-dimension"); + args.AddOption(&ny, "-ny", "--yelems", "Elements in y-dimension"); + args.AddOption(&nz, "-nz", "--zelems", "Elements in z-dimension"); + args.AddOption(&rs_levels, "-rs", "--refine-serial", + "Number of times to refine the mesh uniformly in serial."); + args.AddOption(&rp_levels, "-rp", "--refine-parallel", + "Number of times to refine the mesh uniformly in parallel."); + args.AddOption(&order_v, "-ok", "--order-kinematic", + "Order (degree) of the kinematic finite element space."); + args.AddOption(&order_e, "-ot", "--order-thermo", + "Order (degree) of the thermodynamic finite element space."); + args.AddOption(&order_q, "-oq", "--order-intrule", + "Order of the integration rule."); + args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0."); + args.AddOption(&basename, "-k", "--outputfilename", + "Name of the visit dump files"); + args.Parse(); + if (!args.Good()) { + if (Mpi::Root()) { + args.PrintUsage(std::cout); + } + return 1; + } + if (Mpi::Root()) { + args.PrintOptions(std::cout); + } + + Device::SetMemoryTypes(MemoryType::HOST, MemoryType::DEVICE); + + // On all processors, use the default builtin 1D/2D/3D mesh or read the + // serial one given on the command line. + std::unique_ptr mesh; + switch (dim) { + case 1: + mesh.reset(new Mesh(Mesh::MakeCartesian1D(nx))); + break; + case 2: + mesh.reset(new Mesh( + Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true))); + break; + case 3: + mesh.reset(new Mesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, + 1_r, 1_r, 1_r, true))); + break; + default: + if (Mpi::Root()) { + std::cout << "Invalid number of dims" << std::endl; + } + return -1; + } + + // Refine the mesh in serial to increase the resolution. + for (int lev = 0; lev < rs_levels; lev++) { + mesh->UniformRefinement(); + } + const int mesh_NE = mesh->GetNE(); + if (Mpi::Root()) { + std::cout << "Number of zones in the serial mesh: " << mesh_NE << std::endl; + } + + // Parallel partitioning of the mesh. + std::unique_ptr pmesh; + const int num_tasks = Mpi::WorldSize(); + + if (myid == 0) { + std::cout << "Non-Cartesian partitioning through METIS will be used.\n"; +#ifndef MFEM_USE_METIS + std::cout << "MFEM was built without METIS. " + << "Adjust the number of tasks to use a Cartesian split." + << std::endl; +#endif + } +#ifndef MFEM_USE_METIS + return 1; +#endif + pmesh.reset(new ParMesh(MPI_COMM_WORLD, *mesh)); + + // Refine the mesh further in parallel to increase the resolution. + for (int lev = 0; lev < rp_levels; lev++) { + pmesh->UniformRefinement(); + } + + int NE = pmesh->GetNE(), ne_min, ne_max; + MPI_Reduce(&NE, &ne_min, 1, MPI_INT, MPI_MIN, 0, pmesh->GetComm()); + MPI_Reduce(&NE, &ne_max, 1, MPI_INT, MPI_MAX, 0, pmesh->GetComm()); + if (myid == 0) { + std::cout << "Zones min/max: " << ne_min << " " << ne_max << std::endl; + } + + if (order_q <= 0) { + order_q = (std::max(order_v, order_e) + 1) * 2; + } + + const IntegrationRule &irule = + IntRules.Get(pmesh->GetTypicalElementGeometry(), order_q); + + QuadratureSpace qspace(*pmesh, irule); + QuadratureFunction qfunc(qspace, 2 + dim); + + double gamma = 1.4; + double blast_energy = 1; + double blast_position[] = {0.0, 0.0, 0.0}; + double rho0 = 1; + double omega = 0; + { + SedovSol asol(dim, gamma, rho0, blast_energy, omega); + asol.SetTime(t_final); + if (myid == 0) { + std::cout << "a = " << asol.a << std::endl; + std::cout << "b = " << asol.b << std::endl; + std::cout << "c = " << asol.c << std::endl; + std::cout << "d = " << asol.d << std::endl; + std::cout << "e = " << asol.e << std::endl; + + std::cout << "alpha0 = " << asol.alpha0 << std::endl; + std::cout << "alpha1 = " << asol.alpha1 << std::endl; + std::cout << "alpha2 = " << asol.alpha2 << std::endl; + std::cout << "alpha3 = " << asol.alpha3 << std::endl; + std::cout << "alpha4 = " << asol.alpha4 << std::endl; + std::cout << "alpha5 = " << asol.alpha5 << std::endl; + + std::cout << "V0 = " << asol.V0 << std::endl; + std::cout << "Vv = " << asol.Vv << std::endl; + std::cout << "V2 = " << asol.V2 << std::endl; + std::cout << "Vs = " << asol.Vs << std::endl; + std::cout << "alpha = " << asol.alpha << std::endl; + + std::cout << "r2 (shock position) = " << asol.r2 << std::endl; + std::cout << "U (shock speed) = " << asol.U << std::endl; + std::cout << "rho1 (pre-shock density) = " << asol.rho1 << std::endl; + std::cout << "rho2 (post-shock density) = " << asol.rho2 << std::endl; + std::cout << "v2 (post-shock velocity) = " << asol.v2 << std::endl; + std::cout << "p2 (post-shock pressure) = " << asol.p2 << std::endl; + } + auto slambda = [&](const Vector &x, Vector &res) { + real_t tmp[3]; + Vector dr(tmp, dim); + double r = 0; + + for (int i = 0; i < dim; ++i) { + dr[i] = x[i] - blast_position[i]; + r += dr[i] * dr[i]; + } + r = sqrt(r); + if (r) { + for (int i = 0; i < dim; ++i) { + dr[i] /= r; + } + } + else + { + dr = 0_r; + } + double rho, v, P; + asol.EvalSol(r, rho, v, P); + res[0] = rho; + for (int i = 0; i < dim; ++i) { + res[1 + i] = v * dr[i]; + } + // internal energy + res[1 + dim] = P / (gamma - 1); + }; + VectorFunctionCoefficient asol_coeff(2 + dim, slambda); + asol_coeff.Project(qfunc); + } + + // Define the parallel finite element spaces. We use: + // - H1 (Gauss-Lobatto, continuous) for position and velocity. + // - L2 (Bernstein, discontinuous) for specific internal energy. + L2_FECollection L2FEC(order_e, dim, BasisType::Positive); + H1_FECollection H1FEC(order_v, dim); + ParFiniteElementSpace L2FESpace(pmesh.get(), &L2FEC); + ParFiniteElementSpace H1FESpace(pmesh.get(), &H1FEC, pmesh->Dimension()); + + ParGridFunction rho_gf(&L2FESpace); + ParGridFunction v_gf(&H1FESpace); + ParGridFunction energy_gf(&L2FESpace); +#if 0 + // TODO: need to allow vector LF integrator to specify intrule + VectorQuadratureFunctionCoefficient qcoeff(qfunc); + { + qcoeff.SetComponent(0, 1); + ProjectCoeff(rho_gf, qcoeff, &irule); + } + { + qcoeff.SetComponent(1, dim); + ProjectCoeff(v_gf, qcoeff, &irule); + } + { + qcoeff.SetComponent(1 + dim, 1); + ProjectCoeff(energy_gf, qcoeff, &irule); + } +#endif + + { + std::stringstream fname; + fname << basename << "_mesh"; + pmesh->Save(fname.str().c_str()); + } + { + std::stringstream fname; + fname << basename << "_qfunc"; + std::ofstream out(fname.str()); + qfunc.Save(out); + } +#if 0 + { + std::stringstream fname; + fname << basename << "_rho"; + rho_gf.Save(fname.str().c_str()); + } + { + std::stringstream fname; + fname << basename << "_v"; + v_gf.Save(fname.str().c_str()); + } + { + std::stringstream fname; + fname << basename << "_energy"; + energy_gf.Save(fname.str().c_str()); + } +#endif + + return 0; +} + +// static void ProjectCoeff(ParGridFunction &u, VectorCoefficient &coeff, +// const IntegrationRule *ir) { +// LinearForm b(u.FESpace()); +// b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(coeff, ir)); +// b.UseFastAssembly(true); +// b.Assemble(); + +// BilinearForm a(u.FESpace()); +// a.SetAssemblyLevel(AssemblyLevel::FULL); +// a.AddDomainIntegrator(new VectorFEMassIntegrator()); +// a.Assemble(); +// // Set solver and preconditioner +// SparseMatrix A(a.SpMat()); +// GSSmoother prec(A); +// CGSolver cg; +// cg.SetPreconditioner(prec); +// cg.SetOperator(A); +// cg.SetRelTol(1e-12); +// cg.SetMaxIter(1000); +// cg.SetPrintLevel(0); + +// // Solve and get solution +// u = 0.0; +// cg.Mult(b, u); +// } diff --git a/sedov_sol.cpp b/sedov_sol.cpp new file mode 100644 index 00000000..c4f9cafd --- /dev/null +++ b/sedov_sol.cpp @@ -0,0 +1,173 @@ +// 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 "sedov_sol.hpp" + +#include "adaptive_quad.hpp" +#include "bisect.hpp" + +#include +#include + +#include + +SedovSol::SedovSol(int dim_, double gamma_, double rho_0_, double blast_energy_, + double omega_) + : dim(dim_), gamma(gamma_), rho_0(rho_0_), omega(omega_), + blast_energy(blast_energy_) { + a = (dim + 2 - omega) * (gamma + 1) * 0.25; + b = (gamma + 1) / (gamma - 1); + c = (dim + 2 - omega) * gamma * 0.5; + d = ((dim + 2 - omega) * (gamma + 1) / + ((dim + 2 - omega) * (gamma + 1) - 2 * (2 + dim * (gamma - 1)))); + e = (2 + dim * (gamma - 1)) * 0.5; + + alpha0 = 2. / (dim + 2 - omega); + alpha2 = -(gamma - 1) / (2 * (gamma - 1) + dim - gamma * omega); + alpha1 = + ((dim + 2 - omega) * gamma / (2 + dim * (gamma - 1)) * + (2 * (dim * (2 - gamma) - omega) / pow(gamma * (dim + 2 - omega), 2) - + alpha2)); + alpha3 = (dim - omega) / (2 * (gamma - 1) + dim - dim * omega); + alpha4 = + (dim + 2 - omega) * (dim - omega) * alpha1 / (dim * (2 - gamma) - omega); + alpha5 = (omega * (1 + gamma) - 2 * dim) / (dim * (2 - gamma) - omega); + + V0 = 2. / ((dim + 2 - omega) * gamma); + Vv = 2. / (dim + 2 - omega); + V2 = 4. / ((dim + 2 - omega) * (gamma + 1)); + Vs = 2. / ((gamma - 1) * dim + 2); + + if (V2 == Vs) { + // singular + alpha = (gamma + 1) / (gamma - 1) * pow(2, dim) / + pow(dim * ((gamma - 1) * dim + 2), 2); + if (dim > 1) { + alpha *= M_PI; + } + } else { + // standard or vacuum + auto Vmin = std::min(V0, Vv); + auto J1_integrand = [dim = dim, gamma = gamma, alpha0 = alpha0, + alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, + alpha4 = alpha4, alpha5 = alpha5, a = a, b = b, c = c, + d = d, e = e, omega = omega](double V) { + return -(gamma + 1) / (gamma - 1) * pow(V, 2) * + (alpha0 / V + alpha2 * c / (c * V - 1) - + alpha1 * e / (1 - e * V)) * + pow((pow((a * V), alpha0) * pow((b * (c * V - 1)), alpha2) * + pow((d * (1 - e * V)), alpha1)), + (-(dim + 2 - omega))) * + pow((b * (c * V - 1)), alpha3) * pow((d * (1 - e * V)), alpha4) * + pow((b * (1 - c * V / gamma)), alpha5); + }; + scalar_error_functor err_fun; + err_fun.eps_abs = 1.49e-8; + err_fun.eps_rel = 1.49e-8; + auto J1 = gk21_integrate(J1_integrand, err_fun, Vmin, V2); + + auto J2_integrand = [dim = dim, gamma = gamma, alpha0 = alpha0, + alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, + alpha4 = alpha4, alpha5 = alpha5, a = a, b = b, c = c, + d = d, e = e, omega = omega](double V) { + return -(gamma + 1) / (2 * gamma) * pow(V, 2) * (c * V - gamma) / + (1 - c * V) * + (alpha0 / V + alpha2 * c / (c * V - 1) - + alpha1 * e / (1 - e * V)) * + pow((pow((a * V), alpha0) * pow((b * (c * V - 1)), alpha2) * + pow((d * (1 - e * V)), alpha1)), + (-(dim + 2 - omega))) * + pow((b * (c * V - 1)), alpha3) * pow((d * (1 - e * V)), alpha4) * + pow((b * (1 - c * V / gamma)), alpha5); + }; + auto J2 = gk21_integrate(J2_integrand, err_fun, Vmin, V2); + double I1 = pow(2, dim - 2) * J1; + double I2 = pow(2, (dim - 1)) / (gamma - 1) * J2; + if (dim > 1) { + I1 *= M_PI; + I2 *= M_PI; + } + alpha = I1 + I2; + } +} + +void SedovSol::SetTime(double t_) { + t = t_; + + r2 = pow((blast_energy / (alpha * rho_0)), (1. / (dim + 2 - omega))) * + pow(t, (2. / (dim + 2 - omega))); + U = (2 / (dim + 2 - omega)) * (r2 / t); + rho1 = rho_0 * pow(r2, -omega); + rho2 = ((gamma + 1) / (gamma - 1)) * rho1; + v2 = (2 / (gamma + 1)) * U; + p2 = (2 / (gamma + 1)) * rho1 * U * U; +} + +void SedovSol::EvalSol(double r, double &rho, double &v, double &P) const { + if (r >= r2) { + // pre-shock state + rho = rho_0 * pow(r, -omega); + v = 0; + P = 0; + return; + } + // post-shock state + if (V2 == Vs) { + // singular + rho = rho2 * pow((r / r2), (dim - 2)); + v = v2 * r / r2; + P = p2 * pow((r / r2), dim); + } else { + // find V(r) + auto x1 = [&](double V) { return a * V; }; + auto x2 = [&](double V) { return b * (c * V - 1); }; + auto x3 = [&](double V) { return d * (1 - e * V); }; + auto x4 = [&](double V) { return b * (1 - c * V / gamma); }; + auto lmbda = [&](double V) { + return pow(x1(V), -alpha0) * pow(x2(V), -alpha2) * pow(x3(V), -alpha1); + }; + auto f = [&](double V) { return x1(V) * lmbda(V); }; + auto g = [&](double V) { + return pow(x1(V), alpha0 * omega) * + pow(x2(V), (alpha3 + alpha2 * omega)) * + pow(x3(V), (alpha4 + alpha1 * omega)) * pow(x4(V), alpha5); + }; + auto h = [&](double V) { + return pow(x1(V), (alpha0 * dim)) * + pow(x3(V), (alpha4 + alpha1 * (omega - 2))) * + pow(x4(V), (1 + alpha5)); + }; + double V; + if (V2 < Vs) { + // standard + V = bisection([&](double V_) { return r2 * lmbda(V_) - r; }, V0, V2); + } else { + // vacuum + V = bisection([&](double V_) { return r2 * lmbda(V_) - r; }, Vv, V2); + double r_vacuum = r2 * lmbda(Vv); + if (r <= r_vacuum) { + // vacuum part + rho = 0; + v = 0; + P = 0; + return; + } + } + rho = rho2 * g(V); + v = v2 * f(V); + P = p2 * h(V); + } +} diff --git a/sedov_sol.hpp b/sedov_sol.hpp new file mode 100644 index 00000000..29aa3df7 --- /dev/null +++ b/sedov_sol.hpp @@ -0,0 +1,77 @@ +// 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. + +#ifndef LAGHOS_SEDOV_SOL_HPP +#define LAGHOS_SEDOV_SOL_HPP + +/// Taylor-von Neumann-Sedov blast wave solution +struct SedovSol { + /// 1 for plane wave, 2 for cylinder, 3 for sphere + int dim; + /// time to compute the solution at + double t = 0; + /// ideal gas gamma + double gamma; + /// initial density = rho_0 * pow(r, -omega) + double rho_0; + double omega; + /// initial blast energy + double blast_energy; + + /// currently only supports uniform initial density + /// computed quantities used for computing the solution + /// these values don't depend on time + double a; + double b; + double c; + double d; + double e; + + double alpha0; + double alpha1; + double alpha2; + double alpha3; + double alpha4; + double alpha5; + + double V0; + double Vv; + double V2; + double Vs; + + double alpha; + + /// these values depend on time + /// shock position + double r2; + /// shock speed + double U; + /// pre-shock density + double rho1; + /// post-shock state + double rho2; + double v2; + double p2; + + void SetTime(double t); + + void EvalSol(double r, double &rho, double &v, double &P) const; + + SedovSol(int dim, double gamma, double rho_0, double blast_energy, + double omega = 0); +}; + +#endif From 9cebc3b2c6607a38ae521bee3440f0931ae45ad3 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Mon, 12 Jan 2026 19:53:43 -0800 Subject: [PATCH 04/28] hopefully fix error in computing error --- laghos.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index 0da0c719..1a5384ee 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -1024,11 +1024,11 @@ int main(int argc, char *argv[]) }; VectorFunctionCoefficient asol_coeff(1, slambda); asol_coeff.Project(err_qfunc); - err_qfunc -= sim_qfunc; + sim_qfunc.HostRead(); err_qfunc.HostReadWrite(); for (int i = 0; i < err_qfunc.Size(); ++i) { - err_qfunc[i] = pow(err_qfunc[i], 2); + err_qfunc[i] = pow(err_qfunc[i] - AsConst(sim_qfunc)[i], 2); } real_t lrho_err = err_qfunc.Integrate(); if (Mpi::Root()) { From d44a6b89dcc2f87885ab3f5020144bd424de2f6e Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Mon, 12 Jan 2026 20:19:43 -0800 Subject: [PATCH 05/28] always use the CPU version to project seems to be some sort of bug in ProjectGridFunction on GPUs? --- laghos.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/laghos.cpp b/laghos.cpp index 1a5384ee..d2f09bc2 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -999,7 +999,12 @@ int main(int argc, char *argv[]) hydro.ComputeDensity(rho_gf); - sim_qfunc.ProjectGridFunction(rho_gf); + rho_gf.HostReadWrite(); + + { + GridFunctionCoefficient ctmp(&rho_gf); + ctmp.Coefficient::Project(sim_qfunc); + } auto slambda = [&](const Vector &x, Vector &res) { real_t tmp[3]; From c29e9f5c938d84f9e0456e69c0ba0b8493d1f632 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Tue, 13 Jan 2026 11:27:47 -0800 Subject: [PATCH 06/28] Added ability to specify generated domain size --- laghos.cpp | 92 ++++++++++++++++++++++++++++++++++--------------- sedov/sedov.cpp | 10 ++++-- 2 files changed, 71 insertions(+), 31 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index d2f09bc2..f21d6550 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -143,6 +143,8 @@ int main(int argc, char *argv[]) int dev = 0; int dev_pool_size = 4; + real_t Sx = 1, Sy = 1, Sz = 1; + bool enable_nc = true; bool enable_rebalance = true; @@ -155,6 +157,12 @@ int main(int argc, char *argv[]) "Elements in y-dimension (do not specify mesh_file)"); args.AddOption(&nz, "-nz", "--zelems", "Elements in z-dimension (do not specify mesh_file)"); + args.AddOption(&Sx, "-Sx", "--xwidth", + "Domain width in x-dimension (do not specify mesh_file)"); + args.AddOption(&Sy, "-Sy", "--ywidth", + "Domain width in y-dimension (do not specify mesh_file)"); + args.AddOption(&Sz, "-Sz", "--zwidth", + "Domain width in z-dimension (do not specify mesh_file)"); args.AddOption(&rs_levels, "-rs", "--refine-serial", "Number of times to refine the mesh uniformly in serial."); args.AddOption(&rp_levels, "-rp", "--refine-parallel", @@ -292,41 +300,65 @@ int main(int argc, char *argv[]) } else { + if (Mpi::Root()) { + // generated domain checks + if (problem == 3) { + switch (dim) { + case 1: + if (Sx <= 0.5) { + cout << "WARNING: The triple point is initialized at x=0.5. Sx = " + << Sx + << " puts the triple point outside the simulation " + "domain." + << endl; + } + break; + case 3: + if (Sz <= 1.5) { + cout << "WARNING: The triple point is initialized at z=1.5. Sz = " + << Sz + << " puts the triple point outside the simulation " + "domain." + << endl; + } + case 2: + if (Sx <= 1) { + cout << "WARNING: The triple point is initialized at x=1. Sx = " + << Sx + << " puts the triple point outside the simulation " + "domain." + << endl; + } + + if (Sy <= 1.5) { + cout << "WARNING: The triple point is initialized at y=1.5. Sy = " + << Sy + << " puts the triple point outside the simulation " + "domain." + << endl; + } + break; + } + } + } + if (dim == 1) { - mesh = new Mesh(Mesh::MakeCartesian1D(2)); + mesh = new Mesh(Mesh::MakeCartesian1D(nx, Sx)); mesh->GetBdrElement(0)->SetAttribute(1); mesh->GetBdrElement(1)->SetAttribute(1); } if (dim == 2) { - switch (problem) { - case 3: - mesh = new Mesh(Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, - true, 7_r, 3_r)); - AssignMeshBdrAttrs2D(*mesh, 0_r, 7_r); - break; - default: - mesh = new Mesh( - Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true)); - AssignMeshBdrAttrs2D(*mesh, 0_r, 1_r); - break; - } + mesh = new Mesh(Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, + true, Sx, Sy)); + AssignMeshBdrAttrs2D(*mesh, 0_r, Sx); } if (dim == 3) { - switch (problem) { - case 3: - mesh = new Mesh(Mesh::MakeCartesian3D( - nx, ny, nz, Element::HEXAHEDRON, 7_r, 3_r, 3_r, true)); - AssignMeshBdrAttrs3D(*mesh, 0_r, 7_r, 0_r, 3_r); - break; - default: - mesh = new Mesh(Mesh::MakeCartesian3D( - nx, ny, nz, Element::HEXAHEDRON, 1_r, 1_r, 1_r, true)); - AssignMeshBdrAttrs3D(*mesh, 0_r, 1_r, 0_r, 1_r); - break; - } + mesh = new Mesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, + Sx, Sy, Sz, true)); + AssignMeshBdrAttrs3D(*mesh, 0_r, Sx, 0_r, Sy); } } dim = mesh->Dimension(); @@ -984,9 +1016,13 @@ int main(int argc, char *argv[]) asol.SetTime(t_final); - MFEM_VERIFY(asol.r2 <= 0.9, - "Solution reflections off boundaries detected, cannot compare " - "against exact solution."); + if (strncmp(mesh_file, "default", 7) == 0) { + real_t min_r = std::min(std::min(Sx, Sy), Sz); + MFEM_VERIFY( + asol.r2 <= min_r, + "Solution reflections off boundaries detected, cannot compare " + "against exact solution."); + } int err_order = std::max((std::max(order_v, order_e) + 1) * 2, order_q) * 2; const IntegrationRule &irule = diff --git a/sedov/sedov.cpp b/sedov/sedov.cpp index a93b58a2..7baf96a9 100644 --- a/sedov/sedov.cpp +++ b/sedov/sedov.cpp @@ -50,12 +50,16 @@ int main(int argc, char *argv[]) { int order_q = -1; double t_final = 0.6; const char *basename = "results/Sedov"; + real_t Sx = 1, Sy = 1, Sz = 1; OptionsParser args(argc, argv); args.AddOption(&dim, "-dim", "--dimension", "Dimension of the problem."); args.AddOption(&nx, "-nx", "--xelems", "Elements in x-dimension"); args.AddOption(&ny, "-ny", "--yelems", "Elements in y-dimension"); args.AddOption(&nz, "-nz", "--zelems", "Elements in z-dimension"); + args.AddOption(&Sx, "-Sx", "--xwidth", "Domain width in x-dimension"); + args.AddOption(&Sy, "-Sy", "--ywidth", "Domain width in y-dimension"); + args.AddOption(&Sz, "-Sz", "--zwidth", "Domain width in z-dimension"); args.AddOption(&rs_levels, "-rs", "--refine-serial", "Number of times to refine the mesh uniformly in serial."); args.AddOption(&rp_levels, "-rp", "--refine-parallel", @@ -87,15 +91,15 @@ int main(int argc, char *argv[]) { std::unique_ptr mesh; switch (dim) { case 1: - mesh.reset(new Mesh(Mesh::MakeCartesian1D(nx))); + mesh.reset(new Mesh(Mesh::MakeCartesian1D(nx, Sx))); break; case 2: mesh.reset(new Mesh( - Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true))); + Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true, Sx, Sy))); break; case 3: mesh.reset(new Mesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, - 1_r, 1_r, 1_r, true))); + Sx, Sy, Sz, true))); break; default: if (Mpi::Root()) { From 67eaf51c0a59cee60283b7b3c923f14235927c67 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Wed, 14 Jan 2026 20:09:19 -0800 Subject: [PATCH 07/28] fixed exact solution --- sedov/sedov.cpp | 1 + sedov_sol.cpp | 27 +++++++++++++++------------ 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/sedov/sedov.cpp b/sedov/sedov.cpp index 7baf96a9..ab6ee063 100644 --- a/sedov/sedov.cpp +++ b/sedov/sedov.cpp @@ -162,6 +162,7 @@ int main(int argc, char *argv[]) { double rho0 = 1; double omega = 0; { + std::cout << std::setprecision(16); SedovSol asol(dim, gamma, rho0, blast_energy, omega); asol.SetTime(t_final); if (myid == 0) { diff --git a/sedov_sol.cpp b/sedov_sol.cpp index c4f9cafd..b9fcea2b 100644 --- a/sedov_sol.cpp +++ b/sedov_sol.cpp @@ -39,7 +39,7 @@ SedovSol::SedovSol(int dim_, double gamma_, double rho_0_, double blast_energy_, alpha2 = -(gamma - 1) / (2 * (gamma - 1) + dim - gamma * omega); alpha1 = ((dim + 2 - omega) * gamma / (2 + dim * (gamma - 1)) * - (2 * (dim * (2 - gamma) - omega) / pow(gamma * (dim + 2 - omega), 2) - + (2 * (dim * (2 - gamma) - omega) / (gamma * pow((dim + 2 - omega), 2)) - alpha2)); alpha3 = (dim - omega) / (2 * (gamma - 1) + dim - dim * omega); alpha4 = @@ -75,25 +75,28 @@ SedovSol::SedovSol(int dim_, double gamma_, double rho_0_, double blast_energy_, pow((b * (1 - c * V / gamma)), alpha5); }; scalar_error_functor err_fun; - err_fun.eps_abs = 1.49e-8; - err_fun.eps_rel = 1.49e-8; - auto J1 = gk21_integrate(J1_integrand, err_fun, Vmin, V2); + err_fun.eps_abs = 1.49e-15; + err_fun.eps_rel = 1.49e-15; + auto J1 = gk21_integrate(J1_integrand, err_fun, Vmin, V2, 20, 64); auto J2_integrand = [dim = dim, gamma = gamma, alpha0 = alpha0, alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, alpha4 = alpha4, alpha5 = alpha5, a = a, b = b, c = c, d = d, e = e, omega = omega](double V) { - return -(gamma + 1) / (2 * gamma) * pow(V, 2) * (c * V - gamma) / - (1 - c * V) * - (alpha0 / V + alpha2 * c / (c * V - 1) - - alpha1 * e / (1 - e * V)) * - pow((pow((a * V), alpha0) * pow((b * (c * V - 1)), alpha2) * - pow((d * (1 - e * V)), alpha1)), - (-(dim + 2 - omega))) * + double denom = 1 - c * V; + if (fabs(denom) <= 1e-15) { + denom = std::copysign(1e-15, denom); + } + return -(gamma + 1) / (2 * gamma) * pow(V, 2) * (c * V - gamma) / denom * + (alpha0 / V + alpha2 * c / -denom - alpha1 * e / (1 - e * V)) * + pow(pow(a * V, alpha0) * pow(b * (c * V - 1), alpha2) * + pow(d * (1 - e * V), alpha1), + -(dim + 2 - omega)) * pow((b * (c * V - 1)), alpha3) * pow((d * (1 - e * V)), alpha4) * pow((b * (1 - c * V / gamma)), alpha5); + }; - auto J2 = gk21_integrate(J2_integrand, err_fun, Vmin, V2); + auto J2 = gk21_integrate(J2_integrand, err_fun, Vmin, V2, 20, 64); double I1 = pow(2, dim - 2) * J1; double I2 = pow(2, (dim - 1)) / (gamma - 1) * J2; if (dim > 1) { From a0c1071570e7a689a158830694d44683db98404a Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 13:50:21 -0800 Subject: [PATCH 08/28] updated baseline to match reduced 3D sedov shock E0 --- makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/makefile b/makefile index b6cf3601..de0712b6 100644 --- a/makefile +++ b/makefile @@ -265,7 +265,7 @@ tests: $(shell echo 'step = 0339, dt = 0.000702, |e| = 4.9695537349e+01' >> BASELINE.dat) $(shell echo 'step = 1041, dt = 0.000121, |e| = 3.3909635545e+03' >> BASELINE.dat) $(shell echo 'step = 1154, dt = 0.001655, |e| = 4.6303396053e+01' >> BASELINE.dat) - $(shell echo 'step = 0560, dt = 0.002449, |e| = 1.3408616722e+02' >> BASELINE.dat) + $(shell echo 'step = 0495, dt = 0.002645, |e| = 7.5737286941e+01' >> BASELINE.dat) $(shell echo 'step = 0413, dt = 0.000470, |e| = 3.2012077410e+01' >> BASELINE.dat) $(shell echo 'step = 2872, dt = 0.000064, |e| = 5.6547039096e+01' >> BASELINE.dat) $(shell echo 'step = 0858, dt = 0.000474, |e| = 5.6691500623e+01' >> BASELINE.dat) From 0bce1479150e246e89f68c928ffcbdbc26d89761 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 13:54:15 -0800 Subject: [PATCH 09/28] last digit different --- makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/makefile b/makefile index de0712b6..2d14f4ff 100644 --- a/makefile +++ b/makefile @@ -265,7 +265,7 @@ tests: $(shell echo 'step = 0339, dt = 0.000702, |e| = 4.9695537349e+01' >> BASELINE.dat) $(shell echo 'step = 1041, dt = 0.000121, |e| = 3.3909635545e+03' >> BASELINE.dat) $(shell echo 'step = 1154, dt = 0.001655, |e| = 4.6303396053e+01' >> BASELINE.dat) - $(shell echo 'step = 0495, dt = 0.002645, |e| = 7.5737286941e+01' >> BASELINE.dat) + $(shell echo 'step = 0495, dt = 0.002645, |e| = 7.5737286940e+01' >> BASELINE.dat) $(shell echo 'step = 0413, dt = 0.000470, |e| = 3.2012077410e+01' >> BASELINE.dat) $(shell echo 'step = 2872, dt = 0.000064, |e| = 5.6547039096e+01' >> BASELINE.dat) $(shell echo 'step = 0858, dt = 0.000474, |e| = 5.6691500623e+01' >> BASELINE.dat) From c22b55f000edb6702e81faaf99e84c98620d04ee Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 14:26:26 -0800 Subject: [PATCH 10/28] should be 1 --- makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/makefile b/makefile index 2d14f4ff..de0712b6 100644 --- a/makefile +++ b/makefile @@ -265,7 +265,7 @@ tests: $(shell echo 'step = 0339, dt = 0.000702, |e| = 4.9695537349e+01' >> BASELINE.dat) $(shell echo 'step = 1041, dt = 0.000121, |e| = 3.3909635545e+03' >> BASELINE.dat) $(shell echo 'step = 1154, dt = 0.001655, |e| = 4.6303396053e+01' >> BASELINE.dat) - $(shell echo 'step = 0495, dt = 0.002645, |e| = 7.5737286940e+01' >> BASELINE.dat) + $(shell echo 'step = 0495, dt = 0.002645, |e| = 7.5737286941e+01' >> BASELINE.dat) $(shell echo 'step = 0413, dt = 0.000470, |e| = 3.2012077410e+01' >> BASELINE.dat) $(shell echo 'step = 2872, dt = 0.000064, |e| = 5.6547039096e+01' >> BASELINE.dat) $(shell echo 'step = 0858, dt = 0.000474, |e| = 5.6691500623e+01' >> BASELINE.dat) From 528466b049c346a907e30d11445e14b7e7070251 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 14:42:04 -0800 Subject: [PATCH 11/28] setup the CMake to also be able to use caliper --- CMakeLists.txt | 12 ++++++++++++ laghos.cpp | 12 ++++++------ makefile | 2 +- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a77949da..8b32a6d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,6 +18,11 @@ find_package(HYPRE REQUIRED) find_package(MFEM REQUIRED) find_package(MPI REQUIRED) +option(LAGHOS_USE_CALIPER "Use Caliper" OFF) +if (LAGHOS_USE_CALIPER) + find_package(caliper REQUIRED) +endif() + set(CMAKE_CXX_STANDARD 17 CACHE STRING "C++ standard to use.") set(CMAKE_CXX_STANDARD_REQUIRED ON CACHE BOOL "Force the use of the chosen C++ standard.") @@ -69,6 +74,13 @@ target_link_libraries(laghos PUBLIC MPI::MPI_CXX PUBLIC HYPRE::HYPRE ) +if(LAGHOS_USE_CALIPER) + target_link_libraries(laghos + PUBLIC caliper + ) + target_compile_definitions(laghos + PUBLIC LAGHOS_USE_CALIPER) +endif() if (MFEM_USE_CUDA) set_source_files_properties(sedov_sol.cpp sedov/sedov.cpp PROPERTIES LANGUAGE CUDA) diff --git a/laghos.cpp b/laghos.cpp index f21d6550..889457dc 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -64,7 +64,7 @@ #include #include "laghos_solver.hpp" #include "sedov_sol.hpp" -#ifdef USE_CALIPER +#ifdef LAGHOS_USE_CALIPER #include #include #endif @@ -256,7 +256,7 @@ int main(int argc, char *argv[]) MFEM_VERIFY(strncmp(mesh_file, "default", 7) == 0, "check: mesh_file"); } -#ifdef USE_CALIPER +#ifdef LAGHOS_USE_CALIPER cali_config_set("CALI_CALIPER_ATTRIBUTE_DEFAULT_SCOPE", "process"); CALI_CXX_MARK_FUNCTION; @@ -773,13 +773,13 @@ int main(int argc, char *argv[]) // } // -#ifdef USE_CALIPER +#ifdef LAGHOS_USE_CALIPER CALI_CXX_MARK_LOOP_BEGIN(mainloop_annotation, "timestep loop"); #endif int ti = 1; for (; !last_step; ti++) { -#ifdef USE_CALIPER +#ifdef LAGHOS_USE_CALIPER CALI_CXX_MARK_LOOP_ITERATION(mainloop_annotation, static_cast(ti)); #endif if (t + dt >= t_final) @@ -942,7 +942,7 @@ int main(int argc, char *argv[]) Checks(ti, e_norm, checks); } } -#ifdef USE_CALIPER +#ifdef LAGHOS_USE_CALIPER CALI_CXX_MARK_LOOP_END(mainloop_annotation); adiak::value("steps", ti); #endif @@ -1002,7 +1002,7 @@ int main(int argc, char *argv[]) vis_e.close(); } -#ifdef USE_CALIPER +#ifdef LAGHOS_USE_CALIPER adiak::fini(); #endif diff --git a/makefile b/makefile index de0712b6..95f51371 100644 --- a/makefile +++ b/makefile @@ -81,7 +81,7 @@ ifeq ($(wildcard $(CALIPER_DIR)),) else CALIPER_INCLUDE := -I $(CALIPER_DIR)/include CALIPER_LIBS := -L $(CALIPER_DIR)/lib64 -lcaliper - CALIPER_FLAGS := -DUSE_CALIPER + CALIPER_FLAGS := -DLAGHOS_USE_CALIPER endif # Only configure Adiak if Caliper is enabled From 87ff6c1081baced7ceb47821a768003be6e31012 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 16:52:08 -0800 Subject: [PATCH 12/28] switch back to E0=2 for 3D test since that was consistent --- laghos.cpp | 4 +++- makefile | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index 889457dc..7833b08a 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -143,6 +143,7 @@ int main(int argc, char *argv[]) int dev = 0; int dev_pool_size = 4; + double blast_energy = 1; real_t Sx = 1, Sy = 1, Sz = 1; bool enable_nc = true; @@ -157,6 +158,8 @@ int main(int argc, char *argv[]) "Elements in y-dimension (do not specify mesh_file)"); args.AddOption(&nz, "-nz", "--zelems", "Elements in z-dimension (do not specify mesh_file)"); + args.AddOption(&blast_energy, "-E0", "--blast-energy", + "Sedov initial blast energy (for problem 1)"); args.AddOption(&Sx, "-Sx", "--xwidth", "Domain width in x-dimension (do not specify mesh_file)"); args.AddOption(&Sy, "-Sy", "--ywidth", @@ -642,7 +645,6 @@ int main(int argc, char *argv[]) l2_rho0_gf.ProjectCoefficient(rho0_coeff); rho0_gf.ProjectGridFunction(l2_rho0_gf); - double blast_energy = 1; double blast_position[] = {0.0, 0.0, 0.0}; if (problem == 1) { diff --git a/makefile b/makefile index 95f51371..1423634d 100644 --- a/makefile +++ b/makefile @@ -241,7 +241,7 @@ tests: cat RUN.dat | tail -n 18 | head -n 1 | \ awk '{ printf("step = %04d, dt = %s |e| = %.10e\n", $$2, $$8, $$11); }' >> RESULTS.dat $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP) \ - ./laghos -p 1 -dim 3 -rs 2 -tf 0.6 -pa -vs 100 | tee RUN.dat + ./laghos -p 1 -dim 3 -E0 2 -rs 2 -tf 0.6 -pa -vs 100 | tee RUN.dat cat RUN.dat | tail -n 18 | head -n 1 | \ awk '{ printf("step = %04d, dt = %s |e| = %.10e\n", $$2, $$8, $$11); }' >> RESULTS.dat $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP) \ @@ -265,7 +265,7 @@ tests: $(shell echo 'step = 0339, dt = 0.000702, |e| = 4.9695537349e+01' >> BASELINE.dat) $(shell echo 'step = 1041, dt = 0.000121, |e| = 3.3909635545e+03' >> BASELINE.dat) $(shell echo 'step = 1154, dt = 0.001655, |e| = 4.6303396053e+01' >> BASELINE.dat) - $(shell echo 'step = 0495, dt = 0.002645, |e| = 7.5737286941e+01' >> BASELINE.dat) + $(shell echo 'step = 0560, dt = 0.002449, |e| = 1.3408616722e+02' >> BASELINE.dat) $(shell echo 'step = 0413, dt = 0.000470, |e| = 3.2012077410e+01' >> BASELINE.dat) $(shell echo 'step = 2872, dt = 0.000064, |e| = 5.6547039096e+01' >> BASELINE.dat) $(shell echo 'step = 0858, dt = 0.000474, |e| = 5.6691500623e+01' >> BASELINE.dat) From 91f187c6289518deb11bae92a7ae403196c408ea Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 17:01:09 -0800 Subject: [PATCH 13/28] updated readme --- README.md | 207 ++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 147 insertions(+), 60 deletions(-) diff --git a/README.md b/README.md index 6644435b..7f37d2e9 100644 --- a/README.md +++ b/README.md @@ -128,99 +128,184 @@ Other computational motives in Laghos include the following: Laghos has the following external dependencies: -- *hypre*, used for parallel linear algebra, we recommend version 2.11.2
- https://github.com/hypre-space/hypre/releases/tag/v2.11.2 +- *hypre*, used for parallel linear algebra, we recommend version 2.31.0 or new
+ https://github.com/hypre-space/hypre/releases/tag/v2.31.0 -- METIS, used for parallel domain decomposition (optional), we recommend [version 4.0.3](https://github.com/mfem/tpls/blob/gh-pages/metis-4.0.3.tar.gz)
- https://github.com/mfem/tpls +- METIS, used for parallel domain decomposition (optional) + https://github.com/KarypisLab/METIS.git - MFEM, used for (high-order) finite element discretization, its GitHub master branch
https://github.com/mfem/mfem -To build the miniapp, first download *hypre* and METIS from the links above -and put everything on the same level as the `Laghos` directory: +- Umpire, used for device memory pools in hypre and MFEM. This is only recommended for GPU-accelerated builds. (optional)
+ https://github.com/LLNL/Umpire.git + +- CMake 3.24.0+ +- C and C++17 compiler +- MPI + +This installs built dependencies to `INSTALLDIR` using the `CC` C-compiler and `CXX` C++17-compiler. + +Build METIS (optional): +```sh +git clone https://github.com/KarypisLab/METIS.git +cd METIS +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_C_COMPILER=$CC -DCMAKE_INSTALL_PREFIX=$INSTALLDIR +make -j install +``` +For large runs (problem size above 2 billion unknowns), add `-DMETIS_USE_LONGINDEX=ON` option to the above `cmake` line. If building without METIS only Cartesian partitioning is supported. + +Build Umpire (CUDA, optional): ```sh -~> ls -Laghos/ v2.11.2.tar.gz metis-4.0.3.tar.gz +git clone https://github.com/LLNL/Umpire.git +cd Umpire +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_CUDA_ARCHITECTURES=native -DENABLE_CUDA=ON -DUMPIRE_ENABLE_C=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_CUDA_COMPILER=$CUDACC +make -j install ``` -Build *hypre*: +Build Umpire (HIP, optional): ```sh -~> tar -zxvf v2.11.2.tar.gz -~> cd hypre-2.11.2/src/ -~/hypre-2.11.2/src> ./configure --disable-fortran -~/hypre-2.11.2/src> make -j -~/hypre-2.11.2/src> cd ../.. -~> ln -s hypre-2.11.2 hypre +git clone https://github.com/LLNL/Umpire.git +cd Umpire +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_HIP_ARCHITECTURES=native -DENABLE_HIP=ON -DUMPIRE_ENABLE_C=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_HIP_COMPILER=$HIPCC +make -j install ``` -For large runs (problem size above 2 billion unknowns), add the -`--enable-bigint` option to the above `configure` line. -Build METIS: +Build *hypre* (CPU-only): + +```sh +git clone https://github.com/hypre-space/hypre.git +cd hypre/build +cmake ../src -DCMAKE_BUILD_TYPE=Release -DHYPRE_ENABLE_MIXEDINT=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX +make -j install +``` + +Build *hypre* (CUDA): + +```sh +git clone https://github.com/hypre-space/hypre.git +cd hypre/build +cmake ../src -DCMAKE_BUILD_TYPE=Release -DHYPRE_ENABLE_MIXEDINT=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_ENABLE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=native -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_CUDA_COMPILER=$CUDACC -DHYPRE_ENABLE_GPU_AWARE_MPI=ON -DHYPRE_ENABLE_UMPIRE=ON +make -j install +``` +``HYPRE_ENABLE_GPU_AWARE_MPI`` and ``HYPRE_ENABLE_UMPIRE`` may be optionally turned off. + +Build *hypre* (HIP): + ```sh -~> tar -zxvf metis-4.0.3.tar.gz -~> cd metis-4.0.3 -~/metis-4.0.3> make -~/metis-4.0.3> cd .. -~> ln -s metis-4.0.3 metis-4.0 +git clone https://github.com/hypre-space/hypre.git +cd hypre/build +cmake ../src -DCMAKE_BUILD_TYPE=Release -DHYPRE_ENABLE_MIXEDINT=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_ENABLE_HIP=ON -DCMAKE_HIP_ARCHITECTURES=native -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_HIP_COMPILER=$HIPCC -DHYPRE_ENABLE_GPU_AWARE_MPI=ON -DHYPRE_ENABLE_UMPIRE=ON +make -j install ``` -This build is optional, as MFEM can be build without METIS by specifying -`MFEM_USE_METIS = NO` below. +``HYPRE_ENABLE_GPU_AWARE_MPI`` and ``HYPRE_ENABLE_UMPIRE`` may be optionally turned off. -Clone and build the parallel version of MFEM: +Build MFEM (CPU-only): ```sh -~> git clone https://github.com/mfem/mfem.git ./mfem -~> cd mfem/ -~/mfem> git checkout master -~/mfem> make parallel -j -~/mfem> cd .. +git clone https://github.com/mfem/mfem.git +cd mfem +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_DIR=$INSTALLDIR -DMETIS_DIR=$INSTALLDIR -DMFEM_USE_MPI=ON -DMFEM_USE_METIS=ON -DCMAKE_CXX_COMPILER=$CXX +make -j install ``` -The above uses the `master` branch of MFEM. +`MFEM_USE_METIS` may be optionally disabled. See the [MFEM building page](http://mfem.org/building/) for additional details. -(Optional) Clone and build GLVis: +Build MFEM (CUDA): ```sh -~> git clone https://github.com/GLVis/glvis.git ./glvis -~> cd glvis/ -~/glvis> make -~/glvis> cd .. +git clone https://github.com/mfem/mfem.git +cd mfem +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_DIR=$INSTALLDIR -DMETIS_DIR=$INSTALLDIR -DMFEM_USE_MPI=ON -DMFEM_USE_METIS=ON -DMFEM_USE_CUDA=ON -DMFEM_USE_UMPIRE=ON -DCMAKE_CUDA_ARCHITECTURES=native -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_CUDA_COMPILER=$CUDACC -DUMPIRE_DIR=$INSTALLDIR +make -j install +``` +`MFEM_USE_METIS` and `MFEM_USE_UMPIRE may be optionally disabled. +See the [MFEM building page](http://mfem.org/building/) for additional details. + +Build MFEM (HIP): +```sh +git clone https://github.com/mfem/mfem.git +cd mfem +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_DIR=$INSTALLDIR -DMETIS_DIR=$INSTALLDIR -DMFEM_USE_MPI=ON -DMFEM_USE_METIS=ON -DMFEM_USE_HIP=ON -DMFEM_USE_UMPIRE=ON -DCMAKE_HIP_ARCHITECTURES=native -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_HIP_COMPILER=$HIPCC -DUMPIRE_DIR=$INSTALLDIR +make -j install +``` +`MFEM_USE_METIS` and `MFEM_USE_UMPIRE` may be optionally disabled. +See the [MFEM building page](http://mfem.org/building/) for additional details. + +GLVis (optional): +```sh +git clone https://github.com/GLVis/glvis.git +cd glvis +mkdir build +cd build +cmake .. -DMFEM_DIR=$INSTALLDIR -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_CXX_COMPILER=$CXX ``` The easiest way to visualize Laghos results is to have GLVis running in a separate terminal. Then the `-vis` option in Laghos will stream results directly to the GLVis socket. -(Optional) Build Caliper +Caliper (Optional): 1. Clone and build Adiak: ```sh -~> git clone --recursive https://github.com/LLNL/Adiak.git -~> cd Adiak -~/Adiak> mkdir build && cd build -~/Adiak> cmake -DBUILD_SHARED_LIBS=On -DENABLE_MPI=On \ - -DCMAKE_INSTALL_PREFIX=../../adiak .. -~/Adiak> make && make install -~/Adiak> cd ../.. +git clone --recursive https://github.com/LLNL/Adiak.git +cd Adiak +mkdir build && cd build +cmake -DBUILD_SHARED_LIBS=On -DENABLE_MPI=On \ + -DCMAKE_INSTALL_PREFIX=$INSTALLDIR .. +make -j install ``` 2. Clone and build Caliper: ```sh -~> git clone https://github.com/LLNL/Caliper.git -~> cd Caliper -~/Caliper> mkdir build && cd build -~/Caliper> cmake -DWITH_MPI=True -DWITH_ADIAK=True -Dadiak_ROOT=../../adiak/ \ - -DCMAKE_INSTALL_PREFIX=../../caliper .. -~/Caliper> make && make install -~/Caliper> cd ../.. +git clone https://github.com/LLNL/Caliper.git +cd Caliper +mkdir build && cd build +cmake -DWITH_MPI=True -DWITH_ADIAK=True -Dadiak_ROOT=$INSTALLDIR \ + -DCMAKE_INSTALL_PREFIX=$INSTALLDIR .. +make -j install ``` -Build Laghos +Laghos (CPU-only): ```sh -~> cd Laghos/ -~/Laghos> make -j +git clone https://github.com/CEED/Laghos.git +cd Laghos +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_CXX_COMPILER=$CXX -Dcaliper_ROOT=$INSTALLDIR -DLAGHOS_USE_CALIPER=ON +make -j ``` -This can be followed by `make test` and `make install` to check and install the -build respectively. See `make help` for additional options. +`LAGHOS_USE_CALIPER` may be optionally disabled. -See also the `make setup` target that can be used to automated the -download and building of hypre, METIS and MFEM. +Laghos (CUDA): +```sh +git clone https://github.com/CEED/Laghos.git +cd Laghos +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_CUDA_COMPILER=$CUDACC -DCMAKE_CUDA_ARCHITECTURES=native -Dcaliper_ROOT=$INSTALLDIR -DLAGHOS_USE_CALIPER=ON +make -j +``` +`LAGHOS_USE_CALIPER` may be optionally disabled. + +Laghos (HIP): +```sh +git clone https://github.com/CEED/Laghos.git +cd Laghos +mkdir build +cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_HIP_COMPILER=$HIPCC -DCMAKE_HIP_ARCHITECTURES=native -Dcaliper_ROOT=$INSTALLDIR -DLAGHOS_USE_CALIPER=ON +make -j +``` +`LAGHOS_USE_CALIPER` may be optionally disabled. ## Running @@ -232,13 +317,15 @@ partial assembly option (`-pa`). Some sample runs in 2D and 3D respectively are: ```sh mpirun -np 8 ./laghos -p 1 -dim 2 -rs 3 -tf 0.8 -pa -mpirun -np 8 ./laghos -p 1 -dim 3 -rs 2 -tf 0.6 -pa -vis +mpirun -np 8 ./laghos -p 1 -dim 3 -E0 2 -rs 2 -tf 0.6 -pa -vis ``` The latter produces the following density plot (notice the `-vis` option) [![Sedov blast image](data/sedov.png)](https://glvis.org/live/?stream=../data/laghos.saved) +To compare against the analytical soluton the `-err` option can be used to compute $\int_{\Omega} \|\rho_{sim} - \rho_{exact}\|_2 dV$ of the final solution. + #### Taylor-Green and Gresho vortices Laghos includes also smooth test problems that expose all the principal From 726944ac5f1a391523f1e482d2d80bea7d7b1277 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 19:02:02 -0800 Subject: [PATCH 14/28] debug print out results.dat and baseline.dat --- makefile | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/makefile b/makefile index 1423634d..5bfc4ca0 100644 --- a/makefile +++ b/makefile @@ -270,6 +270,10 @@ tests: $(shell echo 'step = 2872, dt = 0.000064, |e| = 5.6547039096e+01' >> BASELINE.dat) $(shell echo 'step = 0858, dt = 0.000474, |e| = 5.6691500623e+01' >> BASELINE.dat) $(shell echo 'step = 0776, dt = 0.000045, |e| = 4.0982431726e+02' >> BASELINE.dat) + echo "RESULTS.dat:" + cat RESULTS.dat + echo "BASELINE:" + cat BASELINE.dat diff --report-identical-files RESULTS.dat BASELINE.dat # Setup: download & install third party libraries: HYPRE, METIS & MFEM From cd3179b4b6ad5f39417766c8e36373624df673e4 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 19:45:36 -0800 Subject: [PATCH 15/28] makefile bug --- makefile | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/makefile b/makefile index 5bfc4ca0..3b05783e 100644 --- a/makefile +++ b/makefile @@ -261,6 +261,7 @@ tests: -ot 2 -tf 0.62831853 -s 7 -pa -vs 100 | tee RUN.dat cat RUN.dat | tail -n 21 | head -n 1 | \ awk '{ printf("step = %04d, dt = %s |e| = %.10e\n", $$2, $$8, $$11); }' >> RESULTS.dat + cat RESULTS.dat $(shell cat << EOF > BASELINE.dat) $(shell echo 'step = 0339, dt = 0.000702, |e| = 4.9695537349e+01' >> BASELINE.dat) $(shell echo 'step = 1041, dt = 0.000121, |e| = 3.3909635545e+03' >> BASELINE.dat) @@ -270,10 +271,6 @@ tests: $(shell echo 'step = 2872, dt = 0.000064, |e| = 5.6547039096e+01' >> BASELINE.dat) $(shell echo 'step = 0858, dt = 0.000474, |e| = 5.6691500623e+01' >> BASELINE.dat) $(shell echo 'step = 0776, dt = 0.000045, |e| = 4.0982431726e+02' >> BASELINE.dat) - echo "RESULTS.dat:" - cat RESULTS.dat - echo "BASELINE:" - cat BASELINE.dat diff --report-identical-files RESULTS.dat BASELINE.dat # Setup: download & install third party libraries: HYPRE, METIS & MFEM From 6c4f603abf220fdde3de885f00f4a2943b8fc21b Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 20:25:56 -0800 Subject: [PATCH 16/28] now it seems to be passing? --- makefile | 1 - 1 file changed, 1 deletion(-) diff --git a/makefile b/makefile index 3b05783e..1423634d 100644 --- a/makefile +++ b/makefile @@ -261,7 +261,6 @@ tests: -ot 2 -tf 0.62831853 -s 7 -pa -vs 100 | tee RUN.dat cat RUN.dat | tail -n 21 | head -n 1 | \ awk '{ printf("step = %04d, dt = %s |e| = %.10e\n", $$2, $$8, $$11); }' >> RESULTS.dat - cat RESULTS.dat $(shell cat << EOF > BASELINE.dat) $(shell echo 'step = 0339, dt = 0.000702, |e| = 4.9695537349e+01' >> BASELINE.dat) $(shell echo 'step = 1041, dt = 0.000121, |e| = 3.3909635545e+03' >> BASELINE.dat) From 4f15f7b23183da1a27bf8f2a669dbc10d60c93ec Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 20:32:33 -0800 Subject: [PATCH 17/28] mixedint is optional for hypre --- README.md | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 7f37d2e9..38ddf130 100644 --- a/README.md +++ b/README.md @@ -182,7 +182,7 @@ Build *hypre* (CPU-only): ```sh git clone https://github.com/hypre-space/hypre.git cd hypre/build -cmake ../src -DCMAKE_BUILD_TYPE=Release -DHYPRE_ENABLE_MIXEDINT=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX +cmake ../src -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX make -j install ``` @@ -191,7 +191,7 @@ Build *hypre* (CUDA): ```sh git clone https://github.com/hypre-space/hypre.git cd hypre/build -cmake ../src -DCMAKE_BUILD_TYPE=Release -DHYPRE_ENABLE_MIXEDINT=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_ENABLE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=native -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_CUDA_COMPILER=$CUDACC -DHYPRE_ENABLE_GPU_AWARE_MPI=ON -DHYPRE_ENABLE_UMPIRE=ON +cmake ../src -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_ENABLE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=native -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_CUDA_COMPILER=$CUDACC -DHYPRE_ENABLE_GPU_AWARE_MPI=ON -DHYPRE_ENABLE_UMPIRE=ON make -j install ``` ``HYPRE_ENABLE_GPU_AWARE_MPI`` and ``HYPRE_ENABLE_UMPIRE`` may be optionally turned off. @@ -201,10 +201,11 @@ Build *hypre* (HIP): ```sh git clone https://github.com/hypre-space/hypre.git cd hypre/build -cmake ../src -DCMAKE_BUILD_TYPE=Release -DHYPRE_ENABLE_MIXEDINT=ON -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_ENABLE_HIP=ON -DCMAKE_HIP_ARCHITECTURES=native -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_HIP_COMPILER=$HIPCC -DHYPRE_ENABLE_GPU_AWARE_MPI=ON -DHYPRE_ENABLE_UMPIRE=ON +cmake ../src -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=$INSTALLDIR -DHYPRE_ENABLE_HIP=ON -DCMAKE_HIP_ARCHITECTURES=native -DCMAKE_C_COMPILER=$CC -DCMAKE_CXX_COMPILER=$CXX -DCMAKE_HIP_COMPILER=$HIPCC -DHYPRE_ENABLE_GPU_AWARE_MPI=ON -DHYPRE_ENABLE_UMPIRE=ON make -j install ``` ``HYPRE_ENABLE_GPU_AWARE_MPI`` and ``HYPRE_ENABLE_UMPIRE`` may be optionally turned off. +For large runs (problem size above 2 billion unknowns), enable the `HYPRE_ENABLE_MIXEDINT` option. Build MFEM (CPU-only): ```sh From dda1efff0bbe967009ce0fbc24819e3c3371e0fb Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 15 Jan 2026 21:28:19 -0800 Subject: [PATCH 18/28] removed old code --- sedov/sedov.cpp | 87 +------------------------------------------------ 1 file changed, 1 insertion(+), 86 deletions(-) diff --git a/sedov/sedov.cpp b/sedov/sedov.cpp index ab6ee063..d7c4e283 100644 --- a/sedov/sedov.cpp +++ b/sedov/sedov.cpp @@ -29,9 +29,6 @@ using namespace mfem; -// static void ProjectCoeff(ParGridFunction &u, VectorCoefficient &coeff, -// const mfem::IntegrationRule *ir); - int main(int argc, char *argv[]) { // Initialize MPI. Mpi::Init(); @@ -45,9 +42,7 @@ int main(int argc, char *argv[]) { int nx = 2; int ny = 2; int nz = 2; - int order_v = 2; - int order_e = 1; - int order_q = -1; + int order_q = 4; double t_final = 0.6; const char *basename = "results/Sedov"; real_t Sx = 1, Sy = 1, Sz = 1; @@ -64,10 +59,6 @@ int main(int argc, char *argv[]) { "Number of times to refine the mesh uniformly in serial."); args.AddOption(&rp_levels, "-rp", "--refine-parallel", "Number of times to refine the mesh uniformly in parallel."); - args.AddOption(&order_v, "-ok", "--order-kinematic", - "Order (degree) of the kinematic finite element space."); - args.AddOption(&order_e, "-ot", "--order-thermo", - "Order (degree) of the thermodynamic finite element space."); args.AddOption(&order_q, "-oq", "--order-intrule", "Order of the integration rule."); args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0."); @@ -146,10 +137,6 @@ int main(int argc, char *argv[]) { std::cout << "Zones min/max: " << ne_min << " " << ne_max << std::endl; } - if (order_q <= 0) { - order_q = (std::max(order_v, order_e) + 1) * 2; - } - const IntegrationRule &irule = IntRules.Get(pmesh->GetTypicalElementGeometry(), order_q); @@ -224,34 +211,6 @@ int main(int argc, char *argv[]) { asol_coeff.Project(qfunc); } - // Define the parallel finite element spaces. We use: - // - H1 (Gauss-Lobatto, continuous) for position and velocity. - // - L2 (Bernstein, discontinuous) for specific internal energy. - L2_FECollection L2FEC(order_e, dim, BasisType::Positive); - H1_FECollection H1FEC(order_v, dim); - ParFiniteElementSpace L2FESpace(pmesh.get(), &L2FEC); - ParFiniteElementSpace H1FESpace(pmesh.get(), &H1FEC, pmesh->Dimension()); - - ParGridFunction rho_gf(&L2FESpace); - ParGridFunction v_gf(&H1FESpace); - ParGridFunction energy_gf(&L2FESpace); -#if 0 - // TODO: need to allow vector LF integrator to specify intrule - VectorQuadratureFunctionCoefficient qcoeff(qfunc); - { - qcoeff.SetComponent(0, 1); - ProjectCoeff(rho_gf, qcoeff, &irule); - } - { - qcoeff.SetComponent(1, dim); - ProjectCoeff(v_gf, qcoeff, &irule); - } - { - qcoeff.SetComponent(1 + dim, 1); - ProjectCoeff(energy_gf, qcoeff, &irule); - } -#endif - { std::stringstream fname; fname << basename << "_mesh"; @@ -263,49 +222,5 @@ int main(int argc, char *argv[]) { std::ofstream out(fname.str()); qfunc.Save(out); } -#if 0 - { - std::stringstream fname; - fname << basename << "_rho"; - rho_gf.Save(fname.str().c_str()); - } - { - std::stringstream fname; - fname << basename << "_v"; - v_gf.Save(fname.str().c_str()); - } - { - std::stringstream fname; - fname << basename << "_energy"; - energy_gf.Save(fname.str().c_str()); - } -#endif - return 0; } - -// static void ProjectCoeff(ParGridFunction &u, VectorCoefficient &coeff, -// const IntegrationRule *ir) { -// LinearForm b(u.FESpace()); -// b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(coeff, ir)); -// b.UseFastAssembly(true); -// b.Assemble(); - -// BilinearForm a(u.FESpace()); -// a.SetAssemblyLevel(AssemblyLevel::FULL); -// a.AddDomainIntegrator(new VectorFEMassIntegrator()); -// a.Assemble(); -// // Set solver and preconditioner -// SparseMatrix A(a.SpMat()); -// GSSmoother prec(A); -// CGSolver cg; -// cg.SetPreconditioner(prec); -// cg.SetOperator(A); -// cg.SetRelTol(1e-12); -// cg.SetMaxIter(1000); -// cg.SetPrintLevel(0); - -// // Solve and get solution -// u = 0.0; -// cg.Mult(b, u); -// } From d22eda7ef7527ae1f3033e485f0646b493b14d57 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Fri, 16 Jan 2026 00:28:22 -0800 Subject: [PATCH 19/28] formatting, added diagnostic for device memory usage --- adaptive_quad.hpp | 256 ++++++++++++++++---------------- bisect.hpp | 116 +++++++++------ laghos.cpp | 362 ++++++++++++++++++++++++++-------------------- sedov_sol.cpp | 286 +++++++++++++++++++----------------- sedov_sol.hpp | 97 +++++++------ 5 files changed, 614 insertions(+), 503 deletions(-) diff --git a/adaptive_quad.hpp b/adaptive_quad.hpp index d3a27bf1..aae3c007 100644 --- a/adaptive_quad.hpp +++ b/adaptive_quad.hpp @@ -23,138 +23,152 @@ /// /// Implements the 21-point adaptive Gauss-Kronrod quadrature method /// -template struct gk21 { - Fun fun; - Err err_fun; - using res_type = decltype(fun(0.0)); - constexpr static size_t gl_points() { return 10; } - constexpr static size_t gk_points() { return 11; } +template struct gk21 +{ + Fun fun; + Err err_fun; + using res_type = decltype(fun(0.0)); + constexpr static size_t gl_points() { return 10; } + constexpr static size_t gk_points() { return 11; } private: - res_type integrate_recurse(double lower, double upper, size_t curr_depth, - size_t max_depth = 20) const { - static constexpr double data[] = { - // gl_abscissa - -1.488743389816312108848260011297200e-01, - -4.333953941292471907992659431657842e-01, - -6.794095682990244062343273651148736e-01, - -8.650633666889845107320966884234930e-01, - -9.739065285171717200779640120844521e-01, - 1.488743389816312108848260011297200e-01, - 4.333953941292471907992659431657842e-01, - 6.794095682990244062343273651148736e-01, - 8.650633666889845107320966884234930e-01, - 9.739065285171717200779640120844521e-01, - // gl_weights - 2.955242247147528701738929946513383e-01, - 2.692667193099963550912269215694694e-01, - 2.190863625159820439955349342281632e-01, - 1.494513491505805931457763396576973e-01, - 6.667134430868813759356880989333179e-02, - 2.955242247147528701738929946513383e-01, - 2.692667193099963550912269215694694e-01, - 2.190863625159820439955349342281632e-01, - 1.494513491505805931457763396576973e-01, - 6.667134430868813759356880989333179e-02, - // glk_weights - 1.477391049013384913748415159720680e-01, - 1.347092173114733259280540017717068e-01, - 1.093871588022976418992105903258050e-01, - 7.503967481091995276704314091619001e-02, - 3.255816230796472747881897245938976e-02, - 1.477391049013384913748415159720680e-01, - 1.347092173114733259280540017717068e-01, - 1.093871588022976418992105903258050e-01, - 7.503967481091995276704314091619001e-02, - 3.255816230796472747881897245938976e-02, - // gk_abscissa - 0.000000000000000000000000000000000e00, - -2.943928627014601981311266031038656e-01, - -5.627571346686046833390000992726941e-01, - -7.808177265864168970637175783450424e-01, - -9.301574913557082260012071800595083e-01, - -9.956571630258080807355272806890028e-01, - 2.943928627014601981311266031038656e-01, - 5.627571346686046833390000992726941e-01, - 7.808177265864168970637175783450424e-01, - 9.301574913557082260012071800595083e-01, - 9.956571630258080807355272806890028e-01, - // gk_weights - 1.494455540029169056649364683898212e-01, - 1.427759385770600807970942731387171e-01, - 1.234919762620658510779581098310742e-01, - 9.312545458369760553506546508336634e-02, - 5.475589657435199603138130024458018e-02, - 1.169463886737187427806439606219205e-02, - 1.427759385770600807970942731387171e-01, - 1.234919762620658510779581098310742e-01, - 9.312545458369760553506546508336634e-02, - 5.475589657435199603138130024458018e-02, - 1.169463886737187427806439606219205e-02, - }; - // TODO: where to copy gk21_base data to scratch memory? - res_type gl_sum = 0; - res_type gk_sum = 0; - double jac = (upper - lower) * 0.5; - for (int i = 0; i < gl_points(); ++i) { - res_type f_eval = fun((data[i] + 1) * jac + lower); - gl_sum += f_eval * data[gl_points() + i]; - gk_sum += f_eval * data[2 * gl_points() + i]; - } - for (int i = 0; i < gk_points(); ++i) { - res_type f_eval = fun((data[3 * gl_points() + i] + 1) * jac + lower); - gk_sum += f_eval * data[3 * gl_points() + gk_points() + i]; - } - gk_sum *= jac; - gl_sum *= jac; - if (curr_depth < max_depth && !err_fun(gk_sum, gl_sum)) { - gk_sum = integrate_recurse(lower, lower + jac, curr_depth + 1, max_depth); - gk_sum += - integrate_recurse(lower + jac, upper, curr_depth + 1, max_depth); - } + res_type integrate_recurse(double lower, double upper, size_t curr_depth, + size_t max_depth = 20) const + { + static constexpr double data[] = + { + // gl_abscissa + -1.488743389816312108848260011297200e-01, + -4.333953941292471907992659431657842e-01, + -6.794095682990244062343273651148736e-01, + -8.650633666889845107320966884234930e-01, + -9.739065285171717200779640120844521e-01, + 1.488743389816312108848260011297200e-01, + 4.333953941292471907992659431657842e-01, + 6.794095682990244062343273651148736e-01, + 8.650633666889845107320966884234930e-01, + 9.739065285171717200779640120844521e-01, + // gl_weights + 2.955242247147528701738929946513383e-01, + 2.692667193099963550912269215694694e-01, + 2.190863625159820439955349342281632e-01, + 1.494513491505805931457763396576973e-01, + 6.667134430868813759356880989333179e-02, + 2.955242247147528701738929946513383e-01, + 2.692667193099963550912269215694694e-01, + 2.190863625159820439955349342281632e-01, + 1.494513491505805931457763396576973e-01, + 6.667134430868813759356880989333179e-02, + // glk_weights + 1.477391049013384913748415159720680e-01, + 1.347092173114733259280540017717068e-01, + 1.093871588022976418992105903258050e-01, + 7.503967481091995276704314091619001e-02, + 3.255816230796472747881897245938976e-02, + 1.477391049013384913748415159720680e-01, + 1.347092173114733259280540017717068e-01, + 1.093871588022976418992105903258050e-01, + 7.503967481091995276704314091619001e-02, + 3.255816230796472747881897245938976e-02, + // gk_abscissa + 0.000000000000000000000000000000000e00, + -2.943928627014601981311266031038656e-01, + -5.627571346686046833390000992726941e-01, + -7.808177265864168970637175783450424e-01, + -9.301574913557082260012071800595083e-01, + -9.956571630258080807355272806890028e-01, + 2.943928627014601981311266031038656e-01, + 5.627571346686046833390000992726941e-01, + 7.808177265864168970637175783450424e-01, + 9.301574913557082260012071800595083e-01, + 9.956571630258080807355272806890028e-01, + // gk_weights + 1.494455540029169056649364683898212e-01, + 1.427759385770600807970942731387171e-01, + 1.234919762620658510779581098310742e-01, + 9.312545458369760553506546508336634e-02, + 5.475589657435199603138130024458018e-02, + 1.169463886737187427806439606219205e-02, + 1.427759385770600807970942731387171e-01, + 1.234919762620658510779581098310742e-01, + 9.312545458369760553506546508336634e-02, + 5.475589657435199603138130024458018e-02, + 1.169463886737187427806439606219205e-02, + }; + // TODO: where to copy gk21_base data to scratch memory? + res_type gl_sum = 0; + res_type gk_sum = 0; + double jac = (upper - lower) * 0.5; + for (int i = 0; i < gl_points(); ++i) + { + res_type f_eval = fun((data[i] + 1) * jac + lower); + gl_sum += f_eval * data[gl_points() + i]; + gk_sum += f_eval * data[2 * gl_points() + i]; + } + for (int i = 0; i < gk_points(); ++i) + { + res_type f_eval = fun((data[3 * gl_points() + i] + 1) * jac + lower); + gk_sum += f_eval * data[3 * gl_points() + gk_points() + i]; + } + gk_sum *= jac; + gl_sum *= jac; + if (curr_depth < max_depth && !err_fun(gk_sum, gl_sum)) + { + gk_sum = integrate_recurse(lower, lower + jac, curr_depth + 1, max_depth); + gk_sum += + integrate_recurse(lower + jac, upper, curr_depth + 1, max_depth); + } - return gk_sum; - } + return gk_sum; + } public: - res_type integrate(double lower, double upper, size_t start_segs = 1, - size_t max_depth = 20) const { - double dx = (upper - lower) / start_segs; - res_type res = 0; - double curr = lower; - for (size_t i = 0; i < start_segs; ++i) { - double next = lower + (i + 1) * dx; - res += integrate_recurse(curr, next, 1, max_depth); - curr = next; - } - return res; - } + res_type integrate(double lower, double upper, size_t start_segs = 1, + size_t max_depth = 20) const + { + double dx = (upper - lower) / start_segs; + res_type res = 0; + double curr = lower; + for (size_t i = 0; i < start_segs; ++i) + { + double next = lower + (i + 1) * dx; + res += integrate_recurse(curr, next, 1, max_depth); + curr = next; + } + return res; + } }; template auto gk21_integrate(F &&f, E &&e, double lower, double upper, - size_t start_segs = 1, size_t max_depth = 20) { - gk21 integrator{f, e}; - return integrator.integrate(lower, upper, start_segs, max_depth); + size_t start_segs = 1, size_t max_depth = 20) +{ + gk21 integrator{f, e}; + return integrator.integrate(lower, upper, start_segs, max_depth); } -struct scalar_error_functor { - double eps_abs; - double eps_rel; - template bool operator()(const T &ho, const T &lo) const { - if (!std::isfinite(ho)) { - return true; - } - double delta = std::fabs(ho - lo); - if (delta < eps_abs) { - return true; - } - double denom = std::max(std::fabs(ho), std::fabs(lo)); - if (delta < eps_rel * denom) { - return true; - } - return false; - } +struct scalar_error_functor +{ + double eps_abs; + double eps_rel; + template bool operator()(const T &ho, const T &lo) const + { + if (!std::isfinite(ho)) + { + return true; + } + double delta = std::fabs(ho - lo); + if (delta < eps_abs) + { + return true; + } + double denom = std::max(std::fabs(ho), std::fabs(lo)); + if (delta < eps_rel * denom) + { + return true; + } + return false; + } }; #endif diff --git a/bisect.hpp b/bisect.hpp index 2259d9c4..7382164c 100644 --- a/bisect.hpp +++ b/bisect.hpp @@ -23,54 +23,76 @@ #include /// Bisection root finder -template double bisection(Fun &&fun, double lower, double upper) { - double lv = fun(lower); - constexpr double tol = 1e-20; - if (std::fabs(lv) < tol) { - return lower; - } - double rv = fun(upper); - if (std::fabs(rv) < tol) { - return upper; - } - if (std::copysign(1., lv) * std::copysign(1., rv) > 0) { - throw std::runtime_error("bisection: no sign change"); - } - auto dx_init = upper - lower; - auto dx_last = dx_init; - while (true) { - double mid = 0.5 * (lower + upper); - auto dx = mid - lower; - double mv = fun(mid); - if (dx < dx_init * 1e-16 || dx >= dx_last) { - if (fabs(mv) < fabs(lv)) { - if (fabs(mv) < fabs(rv)) { - return mid; - } else if (fabs(rv) < fabs(lv)) { - return upper; - } else { - return lower; - } - } else if (fabs(rv) < fabs(lv)) { - return upper; - } else { - return lower; - } - } - if (std::fabs(mv) < tol) { - return mid; - } - if (std::copysign(1., lv) != std::copysign(1., mv)) { - upper = mid; - rv = mv; - } else if (std::copysign(1., rv) != std::copysign(1., mv)) { - lower = mid; - lv = mv; - } else { +template double bisection(Fun &&fun, double lower, double upper) +{ + double lv = fun(lower); + constexpr double tol = 1e-20; + if (std::fabs(lv) < tol) + { + return lower; + } + double rv = fun(upper); + if (std::fabs(rv) < tol) + { + return upper; + } + if (std::copysign(1., lv) * std::copysign(1., rv) > 0) + { throw std::runtime_error("bisection: no sign change"); - } - dx_last = dx; - } + } + auto dx_init = upper - lower; + auto dx_last = dx_init; + while (true) + { + double mid = 0.5 * (lower + upper); + auto dx = mid - lower; + double mv = fun(mid); + if (dx < dx_init * 1e-16 || dx >= dx_last) + { + if (fabs(mv) < fabs(lv)) + { + if (fabs(mv) < fabs(rv)) + { + return mid; + } + else if (fabs(rv) < fabs(lv)) + { + return upper; + } + else + { + return lower; + } + } + else if (fabs(rv) < fabs(lv)) + { + return upper; + } + else + { + return lower; + } + } + if (std::fabs(mv) < tol) + { + return mid; + } + if (std::copysign(1., lv) != std::copysign(1., mv)) + { + upper = mid; + rv = mv; + } + else if (std::copysign(1., rv) != std::copysign(1., mv)) + { + lower = mid; + lv = mv; + } + else + { + throw std::runtime_error("bisection: no sign change"); + } + dx_last = dx; + } } #endif diff --git a/laghos.cpp b/laghos.cpp index 7833b08a..98c9ae90 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -252,11 +252,12 @@ int main(int argc, char *argv[]) } if (Mpi::Root()) { args.PrintOptions(cout); } - if (check_exact) { - MFEM_VERIFY( + if (check_exact) + { + MFEM_VERIFY( problem == 1, "Can only compare problem 1 (Sedov) against the exact solution"); - MFEM_VERIFY(strncmp(mesh_file, "default", 7) == 0, "check: mesh_file"); + MFEM_VERIFY(strncmp(mesh_file, "default", 7) == 0, "check: mesh_file"); } #ifdef LAGHOS_USE_CALIPER @@ -275,7 +276,8 @@ int main(int argc, char *argv[]) const char * allocator_name = "laghos_device_alloc"; size_t umpire_dev_pool_size = ((size_t) dev_pool_size) * 1024 * 1024 * 1024; size_t umpire_dev_block_size = 512; - rm.makeAllocator(allocator_name, rm.getAllocator("DEVICE"), umpire_dev_pool_size, umpire_dev_block_size); + rm.makeAllocator(allocator_name, + rm.getAllocator("DEVICE"), umpire_dev_pool_size, umpire_dev_block_size); #ifdef HYPRE_USING_UMPIRE HYPRE_SetUmpireDevicePoolName(allocator_name); @@ -303,47 +305,54 @@ int main(int argc, char *argv[]) } else { - if (Mpi::Root()) { - // generated domain checks - if (problem == 3) { - switch (dim) { - case 1: - if (Sx <= 0.5) { - cout << "WARNING: The triple point is initialized at x=0.5. Sx = " - << Sx - << " puts the triple point outside the simulation " - "domain." - << endl; - } - break; - case 3: - if (Sz <= 1.5) { - cout << "WARNING: The triple point is initialized at z=1.5. Sz = " - << Sz - << " puts the triple point outside the simulation " - "domain." - << endl; - } - case 2: - if (Sx <= 1) { - cout << "WARNING: The triple point is initialized at x=1. Sx = " - << Sx - << " puts the triple point outside the simulation " - "domain." - << endl; - } - - if (Sy <= 1.5) { - cout << "WARNING: The triple point is initialized at y=1.5. Sy = " - << Sy - << " puts the triple point outside the simulation " - "domain." - << endl; - } - break; + if (Mpi::Root()) + { + // generated domain checks + if (problem == 3) + { + switch (dim) + { + case 1: + if (Sx <= 0.5) + { + cout << "WARNING: The triple point is initialized at x=0.5. Sx = " + << Sx + << " puts the triple point outside the simulation " + "domain." + << endl; + } + break; + case 3: + if (Sz <= 1.5) + { + cout << "WARNING: The triple point is initialized at z=1.5. Sz = " + << Sz + << " puts the triple point outside the simulation " + "domain." + << endl; + } + case 2: + if (Sx <= 1) + { + cout << "WARNING: The triple point is initialized at x=1. Sx = " + << Sx + << " puts the triple point outside the simulation " + "domain." + << endl; + } + + if (Sy <= 1.5) + { + cout << "WARNING: The triple point is initialized at y=1.5. Sy = " + << Sy + << " puts the triple point outside the simulation " + "domain." + << endl; + } + break; + } } - } - } + } if (dim == 1) { @@ -353,15 +362,15 @@ int main(int argc, char *argv[]) } if (dim == 2) { - mesh = new Mesh(Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, - true, Sx, Sy)); - AssignMeshBdrAttrs2D(*mesh, 0_r, Sx); + mesh = new Mesh(Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, + true, Sx, Sy)); + AssignMeshBdrAttrs2D(*mesh, 0_r, Sx); } if (dim == 3) { - mesh = new Mesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, - Sx, Sy, Sz, true)); - AssignMeshBdrAttrs3D(*mesh, 0_r, Sx, 0_r, Sy); + mesh = new Mesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, + Sx, Sy, Sz, true)); + AssignMeshBdrAttrs3D(*mesh, 0_r, Sx, 0_r, Sy); } } dim = mesh->Dimension(); @@ -644,7 +653,7 @@ int main(int argc, char *argv[]) ParGridFunction l2_rho0_gf(&l2_fes), l2_e(&l2_fes); l2_rho0_gf.ProjectCoefficient(rho0_coeff); rho0_gf.ProjectGridFunction(l2_rho0_gf); - + double blast_position[] = {0.0, 0.0, 0.0}; if (problem == 1) { @@ -751,6 +760,7 @@ int main(int argc, char *argv[]) int steps = 0; BlockVector S_old(S); long mem=0, mmax=0, msum=0; + long dmem = 0, dmmax = 0, dmsum = 0; int checks = 0; // const double internal_energy = hydro.InternalEnergy(e_gf); // const double kinetic_energy = hydro.KineticEnergy(v_gf); @@ -836,6 +846,18 @@ int main(int argc, char *argv[]) if (mem_usage) { mem = GetMaxRssMB(); + size_t mfree, mtot; + if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK)) + { + Device::DeviceMem(&mfree, &mtot); + dmem = mtot - mfree; + MPI_Reduce(&dmem, &dmmax, 1, MPI_LONG, MPI_MAX, 0, + pmesh->GetComm()); + MPI_Reduce(&dmem, &dmsum, 1, MPI_LONG, MPI_SUM, 0, + pmesh->GetComm()); + dmmax /= 1024*1024; + dmsum /= 1024*1024; + } MPI_Reduce(&mem, &mmax, 1, MPI_LONG, MPI_MAX, 0, pmesh->GetComm()); MPI_Reduce(&mem, &msum, 1, MPI_LONG, MPI_SUM, 0, pmesh->GetComm()); } @@ -860,7 +882,8 @@ int main(int argc, char *argv[]) cout << std::fixed; if (mem_usage) { - cout << ", mem: " << mmax << "/" << msum << " MB"; + cout << ", mem: " << mmax << "/" << msum << " MB, " + << dmmax << "/" << dmsum << " MB"; } cout << endl; } @@ -945,8 +968,8 @@ int main(int argc, char *argv[]) } } #ifdef LAGHOS_USE_CALIPER - CALI_CXX_MARK_LOOP_END(mainloop_annotation); - adiak::value("steps", ti); + CALI_CXX_MARK_LOOP_END(mainloop_annotation); + adiak::value("steps", ti); #endif MFEM_VERIFY(!check || checks == 2, "Check error!"); @@ -964,6 +987,16 @@ int main(int argc, char *argv[]) if (mem_usage) { + if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK)) + { + size_t mfree, mtot; + Device::DeviceMem(&mfree, &mtot); + dmem = mtot - mfree; + MPI_Reduce(&dmem, &dmmax, 1, MPI_LONG, MPI_MAX, 0, pmesh->GetComm()); + MPI_Reduce(&dmem, &dmsum, 1, MPI_LONG, MPI_SUM, 0, pmesh->GetComm()); + dmmax /= 1024*1024; + dmsum /= 1024*1024; + } mem = GetMaxRssMB(); MPI_Reduce(&mem, &mmax, 1, MPI_LONG, MPI_MAX, 0, pmesh->GetComm()); MPI_Reduce(&mem, &msum, 1, MPI_LONG, MPI_SUM, 0, pmesh->GetComm()); @@ -978,8 +1011,10 @@ int main(int argc, char *argv[]) << fabs(energy_init - energy_final) << endl; if (mem_usage) { - cout << "Maximum memory resident set size: " - << mmax << "/" << msum << " MB" << endl; + cout << "Maximum memory resident set size: " << mmax << "/" << msum + << "MB , " << dmmax << "/" << dmsum << " MB" << endl; + // cout << "Maximum memory resident set size: " << mmax << "/" << msum + // << " MB" << endl; } } @@ -1008,75 +1043,85 @@ int main(int argc, char *argv[]) adiak::fini(); #endif - if (check_exact) { - // compare against the exact Sedov solution - double gamma = 1.4; - double rho0 = 1; - double omega = 0; + if (check_exact) + { + // compare against the exact Sedov solution + double gamma = 1.4; + double rho0 = 1; + double omega = 0; - SedovSol asol(dim, gamma, rho0, blast_energy, omega); + SedovSol asol(dim, gamma, rho0, blast_energy, omega); - asol.SetTime(t_final); + asol.SetTime(t_final); - if (strncmp(mesh_file, "default", 7) == 0) { - real_t min_r = std::min(std::min(Sx, Sy), Sz); - MFEM_VERIFY( - asol.r2 <= min_r, - "Solution reflections off boundaries detected, cannot compare " - "against exact solution."); - } + if (strncmp(mesh_file, "default", 7) == 0) + { + real_t min_r = std::min(std::min(Sx, Sy), Sz); + MFEM_VERIFY( + asol.r2 <= min_r, + "Solution reflections off boundaries detected, cannot compare " + "against exact solution."); + } - int err_order = std::max((std::max(order_v, order_e) + 1) * 2, order_q) * 2; - const IntegrationRule &irule = + int err_order = std::max((std::max(order_v, order_e) + 1) * 2, order_q) * 2; + const IntegrationRule &irule = IntRules.Get(pmesh->GetTypicalElementGeometry(), err_order); - QuadratureSpace qspace(*pmesh, irule); - // only compare density - QuadratureFunction sim_qfunc(qspace, 1); - QuadratureFunction err_qfunc(qspace, 1); - - hydro.ComputeDensity(rho_gf); - - rho_gf.HostReadWrite(); - - { - GridFunctionCoefficient ctmp(&rho_gf); - ctmp.Coefficient::Project(sim_qfunc); - } - - auto slambda = [&](const Vector &x, Vector &res) { - real_t tmp[3]; - Vector dr(tmp, dim); - double r = 0; - - for (int i = 0; i < dim; ++i) { - dr[i] = x[i] - blast_position[i]; - r += dr[i] * dr[i]; - } - r = sqrt(r); - if (r) { - for (int i = 0; i < dim; ++i) { - dr[i] /= r; + QuadratureSpace qspace(*pmesh, irule); + // only compare density + QuadratureFunction sim_qfunc(qspace, 1); + QuadratureFunction err_qfunc(qspace, 1); + + hydro.ComputeDensity(rho_gf); + + rho_gf.HostReadWrite(); + + { + GridFunctionCoefficient ctmp(&rho_gf); + ctmp.Coefficient::Project(sim_qfunc); + } + + auto slambda = [&](const Vector &x, Vector &res) + { + real_t tmp[3]; + Vector dr(tmp, dim); + double r = 0; + + for (int i = 0; i < dim; ++i) + { + dr[i] = x[i] - blast_position[i]; + r += dr[i] * dr[i]; } - } else { - dr = 0_r; - } - double rho, v, P; - asol.EvalSol(r, rho, v, P); - res[0] = rho; - }; - VectorFunctionCoefficient asol_coeff(1, slambda); - asol_coeff.Project(err_qfunc); - - sim_qfunc.HostRead(); - err_qfunc.HostReadWrite(); - for (int i = 0; i < err_qfunc.Size(); ++i) { - err_qfunc[i] = pow(err_qfunc[i] - AsConst(sim_qfunc)[i], 2); - } - real_t lrho_err = err_qfunc.Integrate(); - if (Mpi::Root()) { - cout << "Density L2 error: " << sqrt(lrho_err) << endl; - } + r = sqrt(r); + if (r) + { + for (int i = 0; i < dim; ++i) + { + dr[i] /= r; + } + } + else + { + dr = 0_r; + } + double rho, v, P; + asol.EvalSol(r, rho, v, P); + res[0] = rho; + }; + VectorFunctionCoefficient asol_coeff(1, slambda); + asol_coeff.Project(err_qfunc); + + sim_qfunc.HostRead(); + err_qfunc.HostReadWrite(); + for (int i = 0; i < err_qfunc.Size(); ++i) + { + err_qfunc[i] = pow(err_qfunc[i] - AsConst(sim_qfunc)[i], 2); + } + real_t lrho_err = err_qfunc.Integrate(); + if (Mpi::Root()) + { + cout << "Density L2 error: " << sqrt(lrho_err) << endl; + } } // Free the used memory. @@ -1347,44 +1392,51 @@ static void Checks(const int ti, const double nrm, int &chk) static void AssignMeshBdrAttrs2D(Mesh& mesh, real_t xmin, real_t xmax) { - Vector pos(3); - constexpr real_t tol = 1e-6; - const int NBE = mesh.GetNBE(); - IntegrationPoint center; - center.x = 0.5; - center.y = 0.5; - center.z = 0.5; - for (int b = 0; b < NBE; b++) { - Element *bel = mesh.GetBdrElement(b); - auto eltrans = mesh.GetBdrElementTransformation(b); - eltrans->Transform(center, pos); - int attr = 2; - if (pos[0] <= xmin + tol || pos[0] >= xmax - tol) { - attr = 1; - } - bel->SetAttribute(attr); - } + Vector pos(3); + constexpr real_t tol = 1e-6; + const int NBE = mesh.GetNBE(); + IntegrationPoint center; + center.x = 0.5; + center.y = 0.5; + center.z = 0.5; + for (int b = 0; b < NBE; b++) + { + Element *bel = mesh.GetBdrElement(b); + auto eltrans = mesh.GetBdrElementTransformation(b); + eltrans->Transform(center, pos); + int attr = 2; + if (pos[0] <= xmin + tol || pos[0] >= xmax - tol) + { + attr = 1; + } + bel->SetAttribute(attr); + } } static void AssignMeshBdrAttrs3D(Mesh &mesh, real_t xmin, real_t xmax, - real_t ymin, real_t ymax) { - Vector pos(3); - constexpr real_t tol = 1e-6; - const int NBE = mesh.GetNBE(); - IntegrationPoint center; - center.x = 0.5; - center.y = 0.5; - center.z = 0.5; - for (int b = 0; b < NBE; b++) { - Element *bel = mesh.GetBdrElement(b); - auto eltrans = mesh.GetBdrElementTransformation(b); - eltrans->Transform(center, pos); - int attr = 3; - if (pos[0] <= xmin + tol || pos[0] >= xmax - tol) { - attr = 1; - } else if (pos[1] <= ymin + tol || pos[1] >= ymax - tol) { - attr = 2; - } - bel->SetAttribute(attr); - } + real_t ymin, real_t ymax) +{ + Vector pos(3); + constexpr real_t tol = 1e-6; + const int NBE = mesh.GetNBE(); + IntegrationPoint center; + center.x = 0.5; + center.y = 0.5; + center.z = 0.5; + for (int b = 0; b < NBE; b++) + { + Element *bel = mesh.GetBdrElement(b); + auto eltrans = mesh.GetBdrElementTransformation(b); + eltrans->Transform(center, pos); + int attr = 3; + if (pos[0] <= xmin + tol || pos[0] >= xmax - tol) + { + attr = 1; + } + else if (pos[1] <= ymin + tol || pos[1] >= ymax - tol) + { + attr = 2; + } + bel->SetAttribute(attr); + } } diff --git a/sedov_sol.cpp b/sedov_sol.cpp index b9fcea2b..7afd99bb 100644 --- a/sedov_sol.cpp +++ b/sedov_sol.cpp @@ -26,151 +26,173 @@ SedovSol::SedovSol(int dim_, double gamma_, double rho_0_, double blast_energy_, double omega_) - : dim(dim_), gamma(gamma_), rho_0(rho_0_), omega(omega_), - blast_energy(blast_energy_) { - a = (dim + 2 - omega) * (gamma + 1) * 0.25; - b = (gamma + 1) / (gamma - 1); - c = (dim + 2 - omega) * gamma * 0.5; - d = ((dim + 2 - omega) * (gamma + 1) / - ((dim + 2 - omega) * (gamma + 1) - 2 * (2 + dim * (gamma - 1)))); - e = (2 + dim * (gamma - 1)) * 0.5; + : dim(dim_), gamma(gamma_), rho_0(rho_0_), omega(omega_), + blast_energy(blast_energy_) +{ + a = (dim + 2 - omega) * (gamma + 1) * 0.25; + b = (gamma + 1) / (gamma - 1); + c = (dim + 2 - omega) * gamma * 0.5; + d = ((dim + 2 - omega) * (gamma + 1) / + ((dim + 2 - omega) * (gamma + 1) - 2 * (2 + dim * (gamma - 1)))); + e = (2 + dim * (gamma - 1)) * 0.5; - alpha0 = 2. / (dim + 2 - omega); - alpha2 = -(gamma - 1) / (2 * (gamma - 1) + dim - gamma * omega); - alpha1 = + alpha0 = 2. / (dim + 2 - omega); + alpha2 = -(gamma - 1) / (2 * (gamma - 1) + dim - gamma * omega); + alpha1 = ((dim + 2 - omega) * gamma / (2 + dim * (gamma - 1)) * (2 * (dim * (2 - gamma) - omega) / (gamma * pow((dim + 2 - omega), 2)) - alpha2)); - alpha3 = (dim - omega) / (2 * (gamma - 1) + dim - dim * omega); - alpha4 = + alpha3 = (dim - omega) / (2 * (gamma - 1) + dim - dim * omega); + alpha4 = (dim + 2 - omega) * (dim - omega) * alpha1 / (dim * (2 - gamma) - omega); - alpha5 = (omega * (1 + gamma) - 2 * dim) / (dim * (2 - gamma) - omega); + alpha5 = (omega * (1 + gamma) - 2 * dim) / (dim * (2 - gamma) - omega); - V0 = 2. / ((dim + 2 - omega) * gamma); - Vv = 2. / (dim + 2 - omega); - V2 = 4. / ((dim + 2 - omega) * (gamma + 1)); - Vs = 2. / ((gamma - 1) * dim + 2); + V0 = 2. / ((dim + 2 - omega) * gamma); + Vv = 2. / (dim + 2 - omega); + V2 = 4. / ((dim + 2 - omega) * (gamma + 1)); + Vs = 2. / ((gamma - 1) * dim + 2); - if (V2 == Vs) { - // singular - alpha = (gamma + 1) / (gamma - 1) * pow(2, dim) / - pow(dim * ((gamma - 1) * dim + 2), 2); - if (dim > 1) { - alpha *= M_PI; - } - } else { - // standard or vacuum - auto Vmin = std::min(V0, Vv); - auto J1_integrand = [dim = dim, gamma = gamma, alpha0 = alpha0, - alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, - alpha4 = alpha4, alpha5 = alpha5, a = a, b = b, c = c, - d = d, e = e, omega = omega](double V) { - return -(gamma + 1) / (gamma - 1) * pow(V, 2) * - (alpha0 / V + alpha2 * c / (c * V - 1) - - alpha1 * e / (1 - e * V)) * - pow((pow((a * V), alpha0) * pow((b * (c * V - 1)), alpha2) * - pow((d * (1 - e * V)), alpha1)), - (-(dim + 2 - omega))) * - pow((b * (c * V - 1)), alpha3) * pow((d * (1 - e * V)), alpha4) * - pow((b * (1 - c * V / gamma)), alpha5); - }; - scalar_error_functor err_fun; - err_fun.eps_abs = 1.49e-15; - err_fun.eps_rel = 1.49e-15; - auto J1 = gk21_integrate(J1_integrand, err_fun, Vmin, V2, 20, 64); + if (V2 == Vs) + { + // singular + alpha = (gamma + 1) / (gamma - 1) * pow(2, dim) / + pow(dim * ((gamma - 1) * dim + 2), 2); + if (dim > 1) + { + alpha *= M_PI; + } + } + else + { + // standard or vacuum + auto Vmin = std::min(V0, Vv); + auto J1_integrand = [dim = dim, gamma = gamma, alpha0 = alpha0, + alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, + alpha4 = alpha4, alpha5 = alpha5, a = a, b = b, c = c, + d = d, e = e, omega = omega](double V) + { + return -(gamma + 1) / (gamma - 1) * pow(V, 2) * + (alpha0 / V + alpha2 * c / (c * V - 1) - + alpha1 * e / (1 - e * V)) * + pow((pow((a * V), alpha0) * pow((b * (c * V - 1)), alpha2) * + pow((d * (1 - e * V)), alpha1)), + (-(dim + 2 - omega))) * + pow((b * (c * V - 1)), alpha3) * pow((d * (1 - e * V)), alpha4) * + pow((b * (1 - c * V / gamma)), alpha5); + }; + scalar_error_functor err_fun; + err_fun.eps_abs = 1.49e-15; + err_fun.eps_rel = 1.49e-15; + auto J1 = gk21_integrate(J1_integrand, err_fun, Vmin, V2, 20, 64); + + auto J2_integrand = [dim = dim, gamma = gamma, alpha0 = alpha0, + alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, + alpha4 = alpha4, alpha5 = alpha5, a = a, b = b, c = c, + d = d, e = e, omega = omega](double V) + { + double denom = 1 - c * V; + if (fabs(denom) <= 1e-15) + { + denom = std::copysign(1e-15, denom); + } + return -(gamma + 1) / (2 * gamma) * pow(V, 2) * (c * V - gamma) / denom * + (alpha0 / V + alpha2 * c / -denom - alpha1 * e / (1 - e * V)) * + pow(pow(a * V, alpha0) * pow(b * (c * V - 1), alpha2) * + pow(d * (1 - e * V), alpha1), + -(dim + 2 - omega)) * + pow((b * (c * V - 1)), alpha3) * pow((d * (1 - e * V)), alpha4) * + pow((b * (1 - c * V / gamma)), alpha5); - auto J2_integrand = [dim = dim, gamma = gamma, alpha0 = alpha0, - alpha1 = alpha1, alpha2 = alpha2, alpha3 = alpha3, - alpha4 = alpha4, alpha5 = alpha5, a = a, b = b, c = c, - d = d, e = e, omega = omega](double V) { - double denom = 1 - c * V; - if (fabs(denom) <= 1e-15) { - denom = std::copysign(1e-15, denom); + }; + auto J2 = gk21_integrate(J2_integrand, err_fun, Vmin, V2, 20, 64); + double I1 = pow(2, dim - 2) * J1; + double I2 = pow(2, (dim - 1)) / (gamma - 1) * J2; + if (dim > 1) + { + I1 *= M_PI; + I2 *= M_PI; } - return -(gamma + 1) / (2 * gamma) * pow(V, 2) * (c * V - gamma) / denom * - (alpha0 / V + alpha2 * c / -denom - alpha1 * e / (1 - e * V)) * - pow(pow(a * V, alpha0) * pow(b * (c * V - 1), alpha2) * - pow(d * (1 - e * V), alpha1), - -(dim + 2 - omega)) * - pow((b * (c * V - 1)), alpha3) * pow((d * (1 - e * V)), alpha4) * - pow((b * (1 - c * V / gamma)), alpha5); - - }; - auto J2 = gk21_integrate(J2_integrand, err_fun, Vmin, V2, 20, 64); - double I1 = pow(2, dim - 2) * J1; - double I2 = pow(2, (dim - 1)) / (gamma - 1) * J2; - if (dim > 1) { - I1 *= M_PI; - I2 *= M_PI; - } - alpha = I1 + I2; - } + alpha = I1 + I2; + } } -void SedovSol::SetTime(double t_) { - t = t_; +void SedovSol::SetTime(double t_) +{ + t = t_; - r2 = pow((blast_energy / (alpha * rho_0)), (1. / (dim + 2 - omega))) * - pow(t, (2. / (dim + 2 - omega))); - U = (2 / (dim + 2 - omega)) * (r2 / t); - rho1 = rho_0 * pow(r2, -omega); - rho2 = ((gamma + 1) / (gamma - 1)) * rho1; - v2 = (2 / (gamma + 1)) * U; - p2 = (2 / (gamma + 1)) * rho1 * U * U; + r2 = pow((blast_energy / (alpha * rho_0)), (1. / (dim + 2 - omega))) * + pow(t, (2. / (dim + 2 - omega))); + U = (2 / (dim + 2 - omega)) * (r2 / t); + rho1 = rho_0 * pow(r2, -omega); + rho2 = ((gamma + 1) / (gamma - 1)) * rho1; + v2 = (2 / (gamma + 1)) * U; + p2 = (2 / (gamma + 1)) * rho1 * U * U; } -void SedovSol::EvalSol(double r, double &rho, double &v, double &P) const { - if (r >= r2) { - // pre-shock state - rho = rho_0 * pow(r, -omega); - v = 0; - P = 0; - return; - } - // post-shock state - if (V2 == Vs) { - // singular - rho = rho2 * pow((r / r2), (dim - 2)); - v = v2 * r / r2; - P = p2 * pow((r / r2), dim); - } else { - // find V(r) - auto x1 = [&](double V) { return a * V; }; - auto x2 = [&](double V) { return b * (c * V - 1); }; - auto x3 = [&](double V) { return d * (1 - e * V); }; - auto x4 = [&](double V) { return b * (1 - c * V / gamma); }; - auto lmbda = [&](double V) { - return pow(x1(V), -alpha0) * pow(x2(V), -alpha2) * pow(x3(V), -alpha1); - }; - auto f = [&](double V) { return x1(V) * lmbda(V); }; - auto g = [&](double V) { - return pow(x1(V), alpha0 * omega) * - pow(x2(V), (alpha3 + alpha2 * omega)) * - pow(x3(V), (alpha4 + alpha1 * omega)) * pow(x4(V), alpha5); - }; - auto h = [&](double V) { - return pow(x1(V), (alpha0 * dim)) * - pow(x3(V), (alpha4 + alpha1 * (omega - 2))) * - pow(x4(V), (1 + alpha5)); - }; - double V; - if (V2 < Vs) { - // standard - V = bisection([&](double V_) { return r2 * lmbda(V_) - r; }, V0, V2); - } else { - // vacuum - V = bisection([&](double V_) { return r2 * lmbda(V_) - r; }, Vv, V2); - double r_vacuum = r2 * lmbda(Vv); - if (r <= r_vacuum) { - // vacuum part - rho = 0; - v = 0; - P = 0; - return; +void SedovSol::EvalSol(double r, double &rho, double &v, double &P) const +{ + if (r >= r2) + { + // pre-shock state + rho = rho_0 * pow(r, -omega); + v = 0; + P = 0; + return; + } + // post-shock state + if (V2 == Vs) + { + // singular + rho = rho2 * pow((r / r2), (dim - 2)); + v = v2 * r / r2; + P = p2 * pow((r / r2), dim); + } + else + { + // find V(r) + auto x1 = [&](double V) { return a * V; }; + auto x2 = [&](double V) { return b * (c * V - 1); }; + auto x3 = [&](double V) { return d * (1 - e * V); }; + auto x4 = [&](double V) { return b * (1 - c * V / gamma); }; + auto lmbda = [&](double V) + { + return pow(x1(V), -alpha0) * pow(x2(V), -alpha2) * pow(x3(V), -alpha1); + }; + auto f = [&](double V) { return x1(V) * lmbda(V); }; + auto g = [&](double V) + { + return pow(x1(V), alpha0 * omega) * + pow(x2(V), (alpha3 + alpha2 * omega)) * + pow(x3(V), (alpha4 + alpha1 * omega)) * pow(x4(V), alpha5); + }; + auto h = [&](double V) + { + return pow(x1(V), (alpha0 * dim)) * + pow(x3(V), (alpha4 + alpha1 * (omega - 2))) * + pow(x4(V), (1 + alpha5)); + }; + double V; + if (V2 < Vs) + { + // standard + V = bisection([&](double V_) { return r2 * lmbda(V_) - r; }, V0, V2); + } + else + { + // vacuum + V = bisection([&](double V_) { return r2 * lmbda(V_) - r; }, Vv, V2); + double r_vacuum = r2 * lmbda(Vv); + if (r <= r_vacuum) + { + // vacuum part + rho = 0; + v = 0; + P = 0; + return; + } } - } - rho = rho2 * g(V); - v = v2 * f(V); - P = p2 * h(V); - } + rho = rho2 * g(V); + v = v2 * f(V); + P = p2 * h(V); + } } diff --git a/sedov_sol.hpp b/sedov_sol.hpp index 29aa3df7..2e558a2a 100644 --- a/sedov_sol.hpp +++ b/sedov_sol.hpp @@ -18,60 +18,61 @@ #define LAGHOS_SEDOV_SOL_HPP /// Taylor-von Neumann-Sedov blast wave solution -struct SedovSol { - /// 1 for plane wave, 2 for cylinder, 3 for sphere - int dim; - /// time to compute the solution at - double t = 0; - /// ideal gas gamma - double gamma; - /// initial density = rho_0 * pow(r, -omega) - double rho_0; - double omega; - /// initial blast energy - double blast_energy; +struct SedovSol +{ + /// 1 for plane wave, 2 for cylinder, 3 for sphere + int dim; + /// time to compute the solution at + double t = 0; + /// ideal gas gamma + double gamma; + /// initial density = rho_0 * pow(r, -omega) + double rho_0; + double omega; + /// initial blast energy + double blast_energy; - /// currently only supports uniform initial density - /// computed quantities used for computing the solution - /// these values don't depend on time - double a; - double b; - double c; - double d; - double e; - - double alpha0; - double alpha1; - double alpha2; - double alpha3; - double alpha4; - double alpha5; - - double V0; - double Vv; - double V2; - double Vs; + /// currently only supports uniform initial density + /// computed quantities used for computing the solution + /// these values don't depend on time + double a; + double b; + double c; + double d; + double e; - double alpha; + double alpha0; + double alpha1; + double alpha2; + double alpha3; + double alpha4; + double alpha5; - /// these values depend on time - /// shock position - double r2; - /// shock speed - double U; - /// pre-shock density - double rho1; - /// post-shock state - double rho2; - double v2; - double p2; + double V0; + double Vv; + double V2; + double Vs; - void SetTime(double t); + double alpha; - void EvalSol(double r, double &rho, double &v, double &P) const; + /// these values depend on time + /// shock position + double r2; + /// shock speed + double U; + /// pre-shock density + double rho1; + /// post-shock state + double rho2; + double v2; + double p2; - SedovSol(int dim, double gamma, double rho_0, double blast_energy, - double omega = 0); + void SetTime(double t); + + void EvalSol(double r, double &rho, double &v, double &P) const; + + SedovSol(int dim, double gamma, double rho_0, double blast_energy, + double omega = 0); }; #endif From 5fd5fd8bb175946d6a7e0d8c46017d3492f32f42 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Fri, 16 Jan 2026 00:33:17 -0800 Subject: [PATCH 20/28] cleanup --- laghos.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index 98c9ae90..312b1ff5 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -1012,9 +1012,7 @@ int main(int argc, char *argv[]) if (mem_usage) { cout << "Maximum memory resident set size: " << mmax << "/" << msum - << "MB , " << dmmax << "/" << dmsum << " MB" << endl; - // cout << "Maximum memory resident set size: " << mmax << "/" << msum - // << " MB" << endl; + << " MB, " << dmmax << "/" << dmsum << " MB" << endl; } } From 113a8473477e44bbca3753d5789455be8305b7a4 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Fri, 16 Jan 2026 12:01:32 -0800 Subject: [PATCH 21/28] updated readme verification command --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 38ddf130..b8d32546 100644 --- a/README.md +++ b/README.md @@ -371,7 +371,7 @@ To make sure the results are correct, we tabulate reference final iterations 1. `mpirun -np 8 ./laghos -p 0 -dim 2 -rs 3 -tf 0.75 -pa` 2. `mpirun -np 8 ./laghos -p 0 -dim 3 -rs 1 -tf 0.75 -pa` 3. `mpirun -np 8 ./laghos -p 1 -dim 2 -rs 3 -tf 0.8 -pa` -4. `mpirun -np 8 ./laghos -p 1 -dim 3 -rs 2 -tf 0.6 -pa` +4. `mpirun -np 8 ./laghos -p 1 -dim 3 -E0 2 -rs 2 -tf 0.6 -pa` 5. `mpirun -np 8 ./laghos -p 2 -dim 1 -rs 5 -tf 0.2 -fa` 6. `mpirun -np 8 ./laghos -p 3 -m data/rectangle01_quad.mesh -rs 2 -tf 3.0 -pa` 7. `mpirun -np 8 ./laghos -p 3 -m data/box01_hex.mesh -rs 1 -tf 5.0 -pa` From fe595316f4d6a14561cd6852830026d60ff34d29 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 29 Jan 2026 09:08:21 -0800 Subject: [PATCH 22/28] makefile build --- README.md | 91 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 90 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index b8d32546..e6c99bf0 100644 --- a/README.md +++ b/README.md @@ -140,10 +140,99 @@ Laghos has the following external dependencies: - Umpire, used for device memory pools in hypre and MFEM. This is only recommended for GPU-accelerated builds. (optional)
https://github.com/LLNL/Umpire.git -- CMake 3.24.0+ +- CMake 3.24.0+ or GNU Make - C and C++17 compiler - MPI +### Makefile + +To build the miniapp, first download *hypre* and METIS from the links above +and put everything on the same level as the `Laghos` directory: +```sh +~> ls +Laghos/ v2.11.2.tar.gz metis-4.0.3.tar.gz +``` + +Build *hypre*: +```sh +~> tar -zxvf v2.11.2.tar.gz +~> cd hypre-2.11.2/src/ +~/hypre-2.11.2/src> ./configure --disable-fortran +~/hypre-2.11.2/src> make -j +~/hypre-2.11.2/src> cd ../.. +~> ln -s hypre-2.11.2 hypre +``` +For large runs (problem size above 2 billion unknowns), add the +`--enable-bigint` option to the above `configure` line. + +Build METIS: +```sh +~> tar -zxvf metis-4.0.3.tar.gz +~> cd metis-4.0.3 +~/metis-4.0.3> make +~/metis-4.0.3> cd .. +~> ln -s metis-4.0.3 metis-4.0 +``` +This build is optional, as MFEM can be build without METIS by specifying +`MFEM_USE_METIS = NO` below. + +Clone and build the parallel version of MFEM: +```sh +~> git clone https://github.com/mfem/mfem.git ./mfem +~> cd mfem/ +~/mfem> git checkout master +~/mfem> make parallel -j +~/mfem> cd .. +``` +The above uses the `master` branch of MFEM. +See the [MFEM building page](http://mfem.org/building/) for additional details. + +(Optional) Clone and build GLVis: +```sh +~> git clone https://github.com/GLVis/glvis.git ./glvis +~> cd glvis/ +~/glvis> make +~/glvis> cd .. +``` +The easiest way to visualize Laghos results is to have GLVis running in a +separate terminal. Then the `-vis` option in Laghos will stream results directly +to the GLVis socket. + +(Optional) Build Caliper +1. Clone and build Adiak: +```sh +~> git clone --recursive https://github.com/LLNL/Adiak.git +~> cd Adiak +~/Adiak> mkdir build && cd build +~/Adiak> cmake -DBUILD_SHARED_LIBS=On -DENABLE_MPI=On \ + -DCMAKE_INSTALL_PREFIX=../../adiak .. +~/Adiak> make && make install +~/Adiak> cd ../.. +``` +2. Clone and build Caliper: +```sh +~> git clone https://github.com/LLNL/Caliper.git +~> cd Caliper +~/Caliper> mkdir build && cd build +~/Caliper> cmake -DWITH_MPI=True -DWITH_ADIAK=True -Dadiak_ROOT=../../adiak/ \ + -DCMAKE_INSTALL_PREFIX=../../caliper .. +~/Caliper> make && make install +~/Caliper> cd ../.. +``` + +Build Laghos +```sh +~> cd Laghos/ +~/Laghos> make -j +``` +This can be followed by `make test` and `make install` to check and install the +build respectively. See `make help` for additional options. + +See also the `make setup` target that can be used to automated the +download and building of hypre, METIS and MFEM. + +### CMake + This installs built dependencies to `INSTALLDIR` using the `CC` C-compiler and `CXX` C++17-compiler. Build METIS (optional): From 44c7a0740a80c8caa0bccfe6b0cac6ef81eda420 Mon Sep 17 00:00:00 2001 From: Vladimir Z Tomov Date: Wed, 18 Feb 2026 17:06:05 -0800 Subject: [PATCH 23/28] Tolerance for the projection of the delta function (Sedov tests). --- laghos.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/laghos.cpp b/laghos.cpp index 312b1ff5..ca7fdf91 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -124,6 +124,7 @@ int main(int argc, char *argv[]) double cfl = 0.5; double cg_tol = 1e-8; double ftz_tol = 0.0; + double delta_tol = 1e-12; int cg_max_iter = 300; int max_tsteps = -1; bool p_assembly = true; @@ -190,6 +191,8 @@ int main(int argc, char *argv[]) "Relative CG tolerance (velocity linear solve)."); args.AddOption(&ftz_tol, "-ftz", "--ftz-tol", "Absolute flush-to-zero tolerance."); + args.AddOption(&delta_tol, "-dtol", "--delta-tol", + "Tolerance for projecting Delta functions."); args.AddOption(&cg_max_iter, "-cgm", "--cg-max-steps", "Maximum number of CG iterations (velocity linear solve)."); args.AddOption(&max_tsteps, "-ms", "--max-steps", @@ -662,7 +665,18 @@ int main(int argc, char *argv[]) // of the symmetric blast. DeltaCoefficient e_coeff(blast_position[0], blast_position[1], blast_position[2], blast_energy / pow(2, dim)); + e_coeff.SetTol(delta_tol); l2_e.ProjectCoefficient(e_coeff); + + int non_finite = l2_e.CheckFinite(); + MPI_Allreduce(MPI_IN_PLACE, &non_finite, 1, MPI_INT, MPI_SUM, pmesh->GetComm()); + if (non_finite > 0) + { + cout << "Delta function coult not be initialized!\n"; + delete ode_solver; + delete pmesh; + return 1; + } } else { From 736dd77f31045c2fac76620e0a6d25bbdd8dbbd9 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 19 Feb 2026 01:26:52 -0800 Subject: [PATCH 24/28] enable serial refinement for metis partitioned mesh --- laghos.cpp | 106 +++++++++++++++++++++++++++++------------------------ 1 file changed, 59 insertions(+), 47 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index 8e258a5f..00f81f66 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -143,9 +143,9 @@ int main(int argc, char *argv[]) args.AddOption(&dim, "-dim", "--dimension", "Dimension of the problem."); args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use."); args.AddOption( - &elem_per_mpi, "-epm", "--elem-per-mpi", - "Number of element per mpi task. Note: this is mutually-exclusive with " - "-nx, -ny, and -nz. Use -epm 0 to use -nx, -ny, and -nz."); + &elem_per_mpi, "-epm", "--elem-per-mpi", + "Number of element per mpi task. Note: this is mutually-exclusive with " + "-nx, -ny, and -nz. Use -epm 0 to use -nx, -ny, and -nz."); args.AddOption(&nx, "-nx", "--xelems", "Elements in x-dimension (do not specify mesh_file). Note: " "this is mutually-exclusive with -nx, -ny, and -nz. Use -epm " @@ -335,50 +335,62 @@ int main(int argc, char *argv[]) } else { - if (elem_per_mpi) { - mesh = PartitionMPI(dim, Mpi::WorldSize(), elem_per_mpi, myid == 0, - rp_levels, mpi_partitioning); - // scale mesh by Sx, Sy, Sz - switch (dim) { - case 1: - mesh.Transform([=](const Vector &x, Vector &y) { y[0] = x[0] * Sx; }); - break; - case 2: - mesh.Transform([=](const Vector &x, Vector &y) { y[0] = x[0] * Sx; - y[1] = x[1] * Sy;}); - break; - case 3: - mesh.Transform([=](const Vector &x, Vector &y) { - y[0] = x[0] * Sx; - y[1] = x[1] * Sy; - y[2] = x[2] * Sz; - }); - break; - } - } else { - if (dim == 1) { - mesh = Mesh::MakeCartesian1D(nx, Sx); - } - if (dim == 2) { - mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true, Sx, - Sy); - } - if (dim == 3) { - mesh = Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, Sx, Sy, - Sz, true); - } - } - - if (dim == 1) { - mesh.GetBdrElement(0)->SetAttribute(1); - mesh.GetBdrElement(1)->SetAttribute(1); - } - if (dim == 2) { - AssignMeshBdrAttrs2D(mesh, 0_r, Sx); - } - if (dim == 3) { - AssignMeshBdrAttrs3D(mesh, 0_r, Sx, 0_r, Sy); - } + if (elem_per_mpi) + { + mesh = PartitionMPI(dim, Mpi::WorldSize(), elem_per_mpi, myid == 0, + rp_levels, mpi_partitioning); + // scale mesh by Sx, Sy, Sz + switch (dim) + { + case 1: + mesh.Transform([=](const Vector &x, Vector &y) { y[0] = x[0] * Sx; }); + mesh.GetBdrElement(0)->SetAttribute(1); + mesh.GetBdrElement(1)->SetAttribute(1); + break; + case 2: + mesh.Transform([=](const Vector &x, Vector &y) + { + y[0] = x[0] * Sx; + y[1] = x[1] * Sy; + }); + AssignMeshBdrAttrs2D(mesh, 0_r, Sx); + break; + case 3: + mesh.Transform([=](const Vector &x, Vector &y) + { + y[0] = x[0] * Sx; + y[1] = x[1] * Sy; + y[2] = x[2] * Sz; + }); + AssignMeshBdrAttrs3D(mesh, 0_r, Sx, 0_r, Sy); + break; + } + } + else + { + if (dim == 1) + { + mesh = Mesh::MakeCartesian1D(nx, Sx); + mesh.GetBdrElement(0)->SetAttribute(1); + mesh.GetBdrElement(1)->SetAttribute(1); + } + if (dim == 2) + { + mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true, Sx, + Sy); + AssignMeshBdrAttrs2D(mesh, 0_r, Sx); + } + if (dim == 3) + { + mesh = Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON, Sx, Sy, + Sz, true); + AssignMeshBdrAttrs3D(mesh, 0_r, Sx, 0_r, Sy); + } + for (int lev = 0; lev < rs_levels; lev++) + { + mesh.UniformRefinement(); + } + } } dim = mesh.Dimension(); From 6e4603ef161dda46ce48051c360522caf0df656d Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 19 Feb 2026 01:48:13 -0800 Subject: [PATCH 25/28] makefile build --- makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/makefile b/makefile index 211e9d73..6f333be2 100644 --- a/makefile +++ b/makefile @@ -121,8 +121,8 @@ CCC = $(strip $(CXX) $(LAGHOS_FLAGS) $(if $(EXTRA_INC_DIR),-I$(EXTRA_INC_DIR))) LAGHOS_LIBS = $(MFEM_LIBS) $(MFEM_EXT_LIBS) $(CALIPER_LIBS) $(ADIAK_LIBS) LIBS = $(strip $(LAGHOS_LIBS) $(LDFLAGS)) -SOURCE_FILES = $(sort $(wildcard *.cpp)) -HEADER_FILES = $(sort $(wildcard *.hpp)) +SOURCE_FILES = $(sort $(wildcard *.cpp) sedov/sedov_sol.cpp) +HEADER_FILES = $(sort $(wildcard *.hpp) $(wildcard sedov/*.hpp)) OBJECT_FILES = $(SOURCE_FILES:.cpp=.o) # Targets From c5ca422f8311b8c6286d90f93281d1c844e95dd4 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Mon, 23 Feb 2026 11:29:31 -0800 Subject: [PATCH 26/28] fixed help text --- laghos.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index 00f81f66..bfba11b8 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -148,15 +148,15 @@ int main(int argc, char *argv[]) "-nx, -ny, and -nz. Use -epm 0 to use -nx, -ny, and -nz."); args.AddOption(&nx, "-nx", "--xelems", "Elements in x-dimension (do not specify mesh_file). Note: " - "this is mutually-exclusive with -nx, -ny, and -nz. Use -epm " + "this is mutually-exclusive with -epm. Use -epm " "0 to use -nx, -ny, and -nz."); args.AddOption(&ny, "-ny", "--yelems", "Elements in y-dimension (do not specify mesh_file). Note: " - "this is mutually-exclusive with -nx, -ny, and -nz. Use -epm " + "this is mutually-exclusive with -epm. Use -epm " "0 to use -nx, -ny, and -nz."); args.AddOption(&nz, "-nz", "--zelems", "Elements in z-dimension (do not specify mesh_file). Note: " - "this is mutually-exclusive with -nx, -ny, and -nz. Use -epm " + "this is mutually-exclusive with -epm. Use -epm " "0 to use -nx, -ny, and -nz."); args.AddOption(&blast_energy, "-E0", "--blast-energy", "Sedov initial blast energy (for problem 1)"); From ef277ea1d347db4fce1204835a8f960bd6311aec Mon Sep 17 00:00:00 2001 From: camierjs Date: Tue, 24 Feb 2026 09:03:43 -0800 Subject: [PATCH 27/28] Fix check Sedov results with new blast initial dim scaling --- laghos.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/laghos.cpp b/laghos.cpp index bfba11b8..e9c0fa15 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -1266,7 +1266,7 @@ static void Checks(const int ti, const double nrm, int &chk) }, { {{5, 1.198510951452527e+03}, {188, 1.199384410059154e+03}}, - {{5, 1.339163718592566e+01}, { 28, 7.521073677397994e+00}}, + {{5, 6.695818592962833e+00}, { 20, 4.267902387082487e+00}}, {{5, 2.041491591302486e+01}, { 59, 3.443180411803796e+01}}, {{5, 1.600000000000000e+01}, { 16, 1.600000000000000e+01}}, {{5, 6.892649884704898e+01}, { 18, 6.893688067534482e+01}}, From 78bf2e7ff36c58b4ab509bb79d77a6d868a18767 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Tue, 24 Feb 2026 10:16:11 -0800 Subject: [PATCH 28/28] review comments --- laghos.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index e9c0fa15..a73d4b9c 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -127,7 +127,7 @@ int main(int argc, char *argv[]) const char *basename = "results/Laghos"; const char *device = "cpu"; bool check = false; - bool check_exact = false; + bool check_exact_sedov = false; bool mem_usage = false; bool fom = false; bool gpu_aware_mpi = false; @@ -215,7 +215,7 @@ int main(int argc, char *argv[]) "Device configuration string, see Device::Configure()."); args.AddOption(&check, "-chk", "--checks", "-no-chk", "--no-checks", "Enable 2D checks."); - args.AddOption(&check_exact, "-err", "--exact-error", "-no-err", + args.AddOption(&check_exact_sedov, "-err", "--exact-error", "-no-err", "--no-exact-error", "Enable comparing the Sedov problem (problem 1) against the " "exact solution."); @@ -239,7 +239,7 @@ int main(int argc, char *argv[]) } if (Mpi::Root()) { args.PrintOptions(cout); } - if (check_exact) + if (check_exact_sedov) { MFEM_VERIFY( problem == 1, @@ -552,7 +552,7 @@ int main(int argc, char *argv[]) MPI_Allreduce(MPI_IN_PLACE, &non_finite, 1, MPI_INT, MPI_SUM, pmesh.GetComm()); if (non_finite > 0) { - cout << "Delta function coult not be initialized!\n"; + cout << "Delta function could not be initialized!\n"; delete ode_solver; return 1; } @@ -935,7 +935,7 @@ int main(int argc, char *argv[]) adiak::fini(); #endif - if (check_exact) + if (check_exact_sedov) { // compare against the exact Sedov solution double gamma = 1.4; @@ -985,7 +985,7 @@ int main(int argc, char *argv[]) r += dr[i] * dr[i]; } r = sqrt(r); - if (r) + if (r > 0) { for (int i = 0; i < dim; ++i) {