From 17a298c0395674e26f19f17664e777abb9eba068 Mon Sep 17 00:00:00 2001 From: Andrew Ho Date: Thu, 30 Jul 2026 17:24:34 -0700 Subject: [PATCH] Enable building Laghos with MFEM_PRECISION=single --- laghos.cpp | 366 +++++++++++++++++++++-------------------- laghos_assembly.cpp | 296 +++++++++++++++++----------------- laghos_assembly.hpp | 4 +- laghos_solver.cpp | 385 +++++++++++++++++++++++--------------------- laghos_solver.hpp | 44 ++--- sedov/sedov.cpp | 2 +- 6 files changed, 565 insertions(+), 532 deletions(-) diff --git a/laghos.cpp b/laghos.cpp index f7eaac5b..43e6bfe0 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -61,30 +61,32 @@ #include #endif +using std::cout; +using std::endl; +using namespace mfem; + #if (defined(HYPRE_USING_UMPIRE) || defined(MFEM_USE_UMPIRE)) && (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)) #define LAGHOS_USE_DEVICE_UMPIRE #define LAGHOS_DEVICE_ALLOCATOR_NAME "LAGHOS_DEVICE_POOL" #include #include -double getDeviceMemoryHighWatermark() { +real_t getDeviceMemoryHighWatermark() +{ auto &rm = umpire::ResourceManager::getInstance(); auto allocator = rm.getAllocator(LAGHOS_DEVICE_ALLOCATOR_NAME); - return ((double) allocator.getHighWatermark()) / (1024 * 1024 * 1024); + return static_cast(allocator.getHighWatermark()) / + (1024.0_r * 1024.0_r * 1024.0_r); } #endif -using std::cout; -using std::endl; -using namespace mfem; - // Choice for the problem setup. static int problem, dim; // Forward declarations. -double e0(const Vector &); -double rho0(const Vector &); -double gamma_func(const Vector &); +real_t e0(const Vector &); +real_t rho0(const Vector &); +real_t gamma_func(const Vector &); void v0(const Vector &, Vector &); static void AssignMeshBdrAttrs2D(Mesh &, real_t, real_t); @@ -92,7 +94,7 @@ 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); +static void Checks(const int ti, const real_t norm, int &checks); static bool HasOption(int argc, char *argv[], const char *short_name, const char *long_name); static bool ValidateElemPerMpiOptions(int elem_per_mpi, int argc, char *argv[], @@ -100,12 +102,13 @@ static bool ValidateElemPerMpiOptions(int elem_per_mpi, int argc, char *argv[], #ifdef LAGHOS_USE_CALIPER static void RecordAdiakMetadata(int dim, const char *mesh_file, int elem_per_mpi, - int nx, int ny, int nz, double blast_energy, - double Sx, double Sy, double Sz, + int nx, int ny, int nz, real_t blast_energy, + real_t Sx, real_t Sy, real_t Sz, int rs_levels, int rp_levels, int problem, int order_v, int order_e, int order_q, - int ode_solver_type, double t_final, double cfl, - double cg_tol, double ftz_tol, double delta_tol, + int ode_solver_type, real_t t_final, + real_t cfl, real_t cg_tol, real_t ftz_tol, + real_t delta_tol, int cg_max_iter, int max_tsteps, bool p_assembly, bool impose_visc, bool visualization, int vis_steps, bool visit, @@ -140,11 +143,11 @@ int main(int argc, char *argv[]) int order_e = 1; int order_q = -1; int ode_solver_type = 4; - double t_final = 0.6; - double cfl = 0.5; - double cg_tol = 1e-8; - double ftz_tol = 0.0; - double delta_tol = 1e-12; + real_t t_final = 0.6_r; + real_t cfl = 0.5_r; + real_t cg_tol = 1e-8_r; + real_t ftz_tol = 0.0_r; + real_t delta_tol = 1e-12_r; int cg_max_iter = 300; int max_tsteps = -1; bool p_assembly = true; @@ -163,8 +166,8 @@ 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; + real_t blast_energy = 1.0_r; + real_t Sx = 1.0_r, Sy = 1.0_r, Sz = 1.0_r; bool enable_nc = true; @@ -519,7 +522,7 @@ int main(int argc, char *argv[]) switch (ode_solver_type) { case 1: ode_solver = new ForwardEulerSolver; break; - case 2: ode_solver = new RK2Solver(0.5); break; + case 2: ode_solver = new RK2Solver(0.5_r); break; case 3: ode_solver = new RK3SSPSolver; break; case 4: ode_solver = new RK4Solver; break; case 6: ode_solver = new RK6Solver; break; @@ -575,7 +578,7 @@ int main(int argc, char *argv[]) v_gf.ProjectCoefficient(v_coeff); for (int i = 0; i < ess_vdofs.Size(); i++) { - v_gf(ess_vdofs[i]) = 0.0; + v_gf(ess_vdofs[i]) = 0.0_r; } // Sync the data location of v_gf with its base, S v_gf.SyncAliasMemory(S); @@ -594,7 +597,7 @@ int main(int argc, char *argv[]) l2_rho0_gf.ProjectCoefficient(rho0_coeff); rho0_gf.ProjectGridFunction(l2_rho0_gf); - double blast_position[] = {0.0, 0.0, 0.0}; + real_t blast_position[] = {0.0_r, 0.0_r, 0.0_r}; if (problem == 1) { // For the Sedov test, we use a delta function at the origin. @@ -661,7 +664,7 @@ int main(int argc, char *argv[]) ParGridFunction rho_gf; if (visualization || visit) { hydro.ComputeDensity(rho_gf); } - const double energy_init = hydro.InternalEnergy(e_gf) + + const real_t energy_init = hydro.InternalEnergy(e_gf) + hydro.KineticEnergy(v_gf); if (visualization) @@ -696,7 +699,7 @@ int main(int argc, char *argv[]) visit_dc.RegisterField("Velocity", &v_gf); visit_dc.RegisterField("Specific Internal Energy", &e_gf); visit_dc.SetCycle(0); - visit_dc.SetTime(0.0); + visit_dc.SetTime(0.0_r); visit_dc.Save(); } @@ -705,15 +708,15 @@ int main(int argc, char *argv[]) // defines the Mult() method that used by the time integrators. ode_solver->Init(hydro); hydro.ResetTimeStepEstimate(); - double t = 0.0, dt = hydro.GetTimeStepEstimate(S), t_old; + real_t t = 0.0_r, dt = hydro.GetTimeStepEstimate(S), t_old; bool last_step = false; 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); + // const real_t internal_energy = hydro.InternalEnergy(e_gf); + // const real_t kinetic_energy = hydro.KineticEnergy(v_gf); // if (mpi.Root()) // { // cout << std::fixed; @@ -760,13 +763,13 @@ int main(int argc, char *argv[]) steps++; // Adaptive time step control. - const double dt_est = hydro.GetTimeStepEstimate(S); + const real_t dt_est = hydro.GetTimeStepEstimate(S); if (dt_est < dt) { // Repeat (solve again) with a decreased time step - decrease of the // time estimate suggests appearance of oscillations. - dt *= 0.85; - if (dt < std::numeric_limits::epsilon()) + dt *= 0.85_r; + if (dt < std::numeric_limits::epsilon()) { MFEM_ABORT("The time step crashed!"); } t = t_old; S = S_old; @@ -775,7 +778,7 @@ int main(int argc, char *argv[]) if (steps < max_tsteps) { last_step = false; } ti--; continue; } - else if (dt_est > 1.25 * dt) { dt *= 1.02; } + else if (dt_est > 1.25_r * dt) { dt *= 1.02_r; } // Ensure the sub-vectors x_gf, v_gf, and e_gf know the location of the // data in S. This operation simply updates the Memory validity flags of @@ -791,8 +794,9 @@ int main(int argc, char *argv[]) if (last_step || (ti % vis_steps) == 0) { - double lnorm = e_gf * e_gf, norm; - MPI_Allreduce(&lnorm, &norm, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm()); + real_t lnorm = e_gf * e_gf, norm; + MPI_Allreduce(&lnorm, &norm, 1, MPITypeMap::mpi_type, + MPI_SUM, pmesh.GetComm()); if (mem_usage) { mem = GetMaxRssMB(); @@ -811,11 +815,11 @@ int main(int argc, char *argv[]) MPI_Reduce(&mem, &mmax, 1, MPI_LONG, MPI_MAX, 0, pmesh.GetComm()); MPI_Reduce(&mem, &msum, 1, MPI_LONG, MPI_SUM, 0, pmesh.GetComm()); } - // const double internal_energy = hydro.InternalEnergy(e_gf); - // const double kinetic_energy = hydro.KineticEnergy(v_gf); + // const real_t internal_energy = hydro.InternalEnergy(e_gf); + // const real_t kinetic_energy = hydro.KineticEnergy(v_gf); if (Mpi::Root()) { - const double sqrt_norm = sqrt(norm); + const real_t sqrt_norm = sqrt(norm); cout << std::fixed; cout << "step " << std::setw(5) << ti @@ -903,15 +907,16 @@ int main(int argc, char *argv[]) // Problems checks if (check) { - double lnorm = e_gf * e_gf, norm; - MPI_Allreduce(&lnorm, &norm, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm()); - const double e_norm = sqrt(norm); + real_t lnorm = e_gf * e_gf, norm; + MPI_Allreduce(&lnorm, &norm, 1, MPITypeMap::mpi_type, + MPI_SUM, pmesh.GetComm()); + const real_t e_norm = sqrt(norm); MFEM_VERIFY(rs_levels == 0 && rp_levels == 0, "check: rs, rp"); MFEM_VERIFY(order_v == 2, "check: order_v"); MFEM_VERIFY(order_e == 1, "check: order_e"); MFEM_VERIFY(ode_solver_type == 4, "check: ode_solver_type"); - MFEM_VERIFY(t_final == 0.6, "check: t_final"); - MFEM_VERIFY(cfl == 0.5, "check: cfl"); + MFEM_VERIFY(t_final == 0.6_r, "check: t_final"); + MFEM_VERIFY(cfl == 0.5_r, "check: cfl"); MFEM_VERIFY(dim == 2 || dim == 3, "check: dimension"); MFEM_VERIFY(std::string(mesh_file) == "data/square01_quad.mesh" || std::string(mesh_file) == "data/cube01_hex.mesh", "check: mesh_file"); @@ -953,7 +958,7 @@ int main(int argc, char *argv[]) MPI_Reduce(&mem, &msum, 1, MPI_LONG, MPI_SUM, 0, pmesh.GetComm()); } - const double energy_final = hydro.InternalEnergy(e_gf) + + const real_t energy_final = hydro.InternalEnergy(e_gf) + hydro.KineticEnergy(v_gf); if (Mpi::Root()) { @@ -971,7 +976,7 @@ int main(int argc, char *argv[]) // For problems 0 and 4 the exact velocity is constant in time. if (problem == 0 || problem == 4) { - const double error_max = v_gf.ComputeMaxError(v_coeff), + const real_t error_max = v_gf.ComputeMaxError(v_coeff), error_l1 = v_gf.ComputeL1Error(v_coeff), error_l2 = v_gf.ComputeL2Error(v_coeff); if (Mpi::Root()) @@ -1091,58 +1096,61 @@ int main(int argc, char *argv[]) return 0; } -double rho0(const Vector &x) +real_t rho0(const Vector &x) { switch (problem) { - case 0: return 1.0; - case 1: return 1.0; - case 2: return (x(0) < 0.5) ? 1.0 : 0.1; - case 3: return (dim == 2) ? (x(0) > 1.0 && x(1) > 1.5) ? 0.125 : 1.0 - : x(0) > 1.0 && ((x(1) < 1.5 && x(2) < 1.5) || - (x(1) > 1.5 && x(2) > 1.5)) ? 0.125 : 1.0; - case 4: return 1.0; + case 0: return 1.0_r; + case 1: return 1.0_r; + case 2: return (x(0) < 0.5_r) ? 1.0_r : 0.1_r; + case 3: return (dim == 2) ? (x(0) > 1.0_r && x(1) > 1.5_r) ? 0.125_r : 1.0_r + : x(0) > 1.0_r && ((x(1) < 1.5_r && x(2) < 1.5_r) || + (x(1) > 1.5_r && x(2) > 1.5_r)) + ? 0.125_r : 1.0_r; + case 4: return 1.0_r; case 5: { - if (x(0) >= 0.5 && x(1) >= 0.5) { return 0.5313; } - if (x(0) < 0.5 && x(1) < 0.5) { return 0.8; } - return 1.0; + if (x(0) >= 0.5_r && x(1) >= 0.5_r) { return 0.5313_r; } + if (x(0) < 0.5_r && x(1) < 0.5_r) { return 0.8_r; } + return 1.0_r; } case 6: { - if (x(0) < 0.5 && x(1) >= 0.5) { return 2.0; } - if (x(0) >= 0.5 && x(1) < 0.5) { return 3.0; } - return 1.0; + if (x(0) < 0.5_r && x(1) >= 0.5_r) { return 2.0_r; } + if (x(0) >= 0.5_r && x(1) < 0.5_r) { return 3.0_r; } + return 1.0_r; } - case 7: return x(1) >= 0.0 ? 2.0 : 1.0; - default: MFEM_ABORT("Bad number given for problem id!"); return 0.0; + case 7: return x(1) >= 0.0_r ? 2.0_r : 1.0_r; + default: MFEM_ABORT("Bad number given for problem id!"); return 0.0_r; } } -double gamma_func(const Vector &x) +real_t gamma_func(const Vector &x) { switch (problem) { - case 0: return 5.0 / 3.0; - case 1: return 1.4; - case 2: return 1.4; + case 0: return 5.0_r / 3.0_r; + case 1: return 1.4_r; + case 2: return 1.4_r; case 3: - if (dim == 1) { return (x(0) > 0.5) ? 1.4 : 1.5; } - else { return (x(0) > 1.0 && x(1) <= 1.5) ? 1.4 : 1.5; } - case 4: return 5.0 / 3.0; - case 5: return 1.4; - case 6: return 1.4; - case 7: return 5.0 / 3.0; - default: MFEM_ABORT("Bad number given for problem id!"); return 0.0; + if (dim == 1) { return (x(0) > 0.5_r) ? 1.4_r : 1.5_r; } + else { return (x(0) > 1.0_r && x(1) <= 1.5_r) ? 1.4_r : 1.5_r; } + case 4: return 5.0_r / 3.0_r; + case 5: return 1.4_r; + case 6: return 1.4_r; + case 7: return 5.0_r / 3.0_r; + default: MFEM_ABORT("Bad number given for problem id!"); return 0.0_r; } } -static double rad(double x, double y) { return sqrt(x*x + y*y); } +static real_t rad(real_t x, real_t y) { return sqrt(x*x + y*y); } void v0(const Vector &x, Vector &v) { - const double atn = dim!=1 ? pow((x(0)*(1.0-x(0))*4*x(1)*(1.0-x(1))*4.0), - 0.4) : 0.0; + const real_t atn = dim != 1 ? + pow((x(0)*(1.0_r-x(0))*4*x(1)*(1.0_r-x(1))*4.0_r), + 0.4_r) : + 0.0_r; switch (problem) { case 0: @@ -1152,125 +1160,136 @@ void v0(const Vector &x, Vector &v) { v(0) *= cos(M_PI*x(2)); v(1) *= cos(M_PI*x(2)); - v(2) = 0.0; + v(2) = 0.0_r; } break; - case 1: v = 0.0; break; - case 2: v = 0.0; break; - case 3: v = 0.0; break; + case 1: v = 0.0_r; break; + case 2: v = 0.0_r; break; + case 3: v = 0.0_r; break; case 4: { - v = 0.0; - const double r = rad(x(0), x(1)); - if (r < 0.2) + v = 0.0_r; + const real_t r = rad(x(0), x(1)); + if (r < 0.2_r) { - v(0) = 5.0 * x(1); - v(1) = -5.0 * x(0); + v(0) = 5.0_r * x(1); + v(1) = -5.0_r * x(0); } - else if (r < 0.4) + else if (r < 0.4_r) { - v(0) = 2.0 * x(1) / r - 5.0 * x(1); - v(1) = -2.0 * x(0) / r + 5.0 * x(0); + v(0) = 2.0_r * x(1) / r - 5.0_r * x(1); + v(1) = -2.0_r * x(0) / r + 5.0_r * x(0); } else { } break; } case 5: { - v = 0.0; - if (x(0) >= 0.5 && x(1) >= 0.5) { v(0)=0.0*atn, v(1)=0.0*atn; return;} - if (x(0) < 0.5 && x(1) >= 0.5) { v(0)=0.7276*atn, v(1)=0.0*atn; return;} - if (x(0) < 0.5 && x(1) < 0.5) { v(0)=0.0*atn, v(1)=0.0*atn; return;} - if (x(0) >= 0.5 && x(1) < 0.5) { v(0)=0.0*atn, v(1)=0.7276*atn; return; } + v = 0.0_r; + if (x(0) >= 0.5_r && x(1) >= 0.5_r) + { v(0)=0.0_r*atn, v(1)=0.0_r*atn; return; } + if (x(0) < 0.5_r && x(1) >= 0.5_r) + { v(0)=0.7276_r*atn, v(1)=0.0_r*atn; return; } + if (x(0) < 0.5_r && x(1) < 0.5_r) + { v(0)=0.0_r*atn, v(1)=0.0_r*atn; return; } + if (x(0) >= 0.5_r && x(1) < 0.5_r) + { v(0)=0.0_r*atn, v(1)=0.7276_r*atn; return; } MFEM_ABORT("Error in problem 5!"); return; } case 6: { - v = 0.0; - if (x(0) >= 0.5 && x(1) >= 0.5) { v(0)=+0.75*atn, v(1)=-0.5*atn; return;} - if (x(0) < 0.5 && x(1) >= 0.5) { v(0)=+0.75*atn, v(1)=+0.5*atn; return;} - if (x(0) < 0.5 && x(1) < 0.5) { v(0)=-0.75*atn, v(1)=+0.5*atn; return;} - if (x(0) >= 0.5 && x(1) < 0.5) { v(0)=-0.75*atn, v(1)=-0.5*atn; return;} + v = 0.0_r; + if (x(0) >= 0.5_r && x(1) >= 0.5_r) + { v(0)=+0.75_r*atn, v(1)=-0.5_r*atn; return; } + if (x(0) < 0.5_r && x(1) >= 0.5_r) + { v(0)=+0.75_r*atn, v(1)=+0.5_r*atn; return; } + if (x(0) < 0.5_r && x(1) < 0.5_r) + { v(0)=-0.75_r*atn, v(1)=+0.5_r*atn; return; } + if (x(0) >= 0.5_r && x(1) < 0.5_r) + { v(0)=-0.75_r*atn, v(1)=-0.5_r*atn; return; } MFEM_ABORT("Error in problem 6!"); return; } case 7: { - v = 0.0; - v(1) = 0.02 * exp(-2*M_PI*x(1)*x(1)) * cos(2*M_PI*x(0)); + v = 0.0_r; + v(1) = 0.02_r * exp(-2.0_r*M_PI*x(1)*x(1)) * cos(2.0_r*M_PI*x(0)); break; } default: MFEM_ABORT("Bad number given for problem id!"); } } -double e0(const Vector &x) +real_t e0(const Vector &x) { switch (problem) { case 0: { - const double denom = 2.0 / 3.0; // (5/3 - 1) * density. - double val; + const real_t denom = 2.0_r / 3.0_r; // (5/3 - 1) * density. + real_t val; if (x.Size() == 2) { - val = 1.0 + (cos(2*M_PI*x(0)) + cos(2*M_PI*x(1))) / 4.0; + val = 1.0_r + (cos(2.0_r*M_PI*x(0)) + cos(2.0_r*M_PI*x(1))) / 4.0_r; } else { - val = 100.0 + ((cos(2*M_PI*x(2)) + 2) * - (cos(2*M_PI*x(0)) + cos(2*M_PI*x(1))) - 2) / 16.0; + val = 100.0_r + ((cos(2.0_r*M_PI*x(2)) + 2.0_r) * + (cos(2.0_r*M_PI*x(0)) + cos(2.0_r*M_PI*x(1))) - + 2.0_r) / 16.0_r; } return val/denom; } - case 1: return 0.0; // This case in initialized in main(). - case 2: return (x(0) < 0.5) ? 1.0 / rho0(x) / (gamma_func(x) - 1.0) - : 0.1 / rho0(x) / (gamma_func(x) - 1.0); - case 3: return (x(0) > 1.0) ? 0.1 / rho0(x) / (gamma_func(x) - 1.0) - : 1.0 / rho0(x) / (gamma_func(x) - 1.0); + case 1: return 0.0_r; // This case in initialized in main(). + case 2: return (x(0) < 0.5_r) ? 1.0_r / rho0(x) / (gamma_func(x) - 1.0_r) + : 0.1_r / rho0(x) / (gamma_func(x) - 1.0_r); + case 3: return (x(0) > 1.0_r) ? 0.1_r / rho0(x) / (gamma_func(x) - 1.0_r) + : 1.0_r / rho0(x) / (gamma_func(x) - 1.0_r); case 4: { - const double r = rad(x(0), x(1)), rsq = x(0) * x(0) + x(1) * x(1); - const double gamma = 5.0 / 3.0; - if (r < 0.2) + const real_t r = rad(x(0), x(1)), + rsq = x(0) * x(0) + x(1) * x(1); + const real_t gamma = 5.0_r / 3.0_r; + if (r < 0.2_r) { - return (5.0 + 25.0 / 2.0 * rsq) / (gamma - 1.0); + return (5.0_r + 25.0_r / 2.0_r * rsq) / (gamma - 1.0_r); } - else if (r < 0.4) + else if (r < 0.4_r) { - const double t1 = 9.0 - 4.0 * log(0.2) + 25.0 / 2.0 * rsq; - const double t2 = 20.0 * r - 4.0 * log(r); - return (t1 - t2) / (gamma - 1.0); + const real_t t1 = 9.0_r - 4.0_r * log(0.2_r) + + 25.0_r / 2.0_r * rsq; + const real_t t2 = 20.0_r * r - 4.0_r * log(r); + return (t1 - t2) / (gamma - 1.0_r); } - else { return (3.0 + 4.0 * log(2.0)) / (gamma - 1.0); } + else { return (3.0_r + 4.0_r * log(2.0_r)) / (gamma - 1.0_r); } } case 5: { - const double irg = 1.0 / rho0(x) / (gamma_func(x) - 1.0); - if (x(0) >= 0.5 && x(1) >= 0.5) { return 0.4 * irg; } - if (x(0) < 0.5 && x(1) >= 0.5) { return 1.0 * irg; } - if (x(0) < 0.5 && x(1) < 0.5) { return 1.0 * irg; } - if (x(0) >= 0.5 && x(1) < 0.5) { return 1.0 * irg; } + const real_t irg = 1.0_r / rho0(x) / (gamma_func(x) - 1.0_r); + if (x(0) >= 0.5_r && x(1) >= 0.5_r) { return 0.4_r * irg; } + if (x(0) < 0.5_r && x(1) >= 0.5_r) { return 1.0_r * irg; } + if (x(0) < 0.5_r && x(1) < 0.5_r) { return 1.0_r * irg; } + if (x(0) >= 0.5_r && x(1) < 0.5_r) { return 1.0_r * irg; } MFEM_ABORT("Error in problem 5!"); - return 0.0; + return 0.0_r; } case 6: { - const double irg = 1.0 / rho0(x) / (gamma_func(x) - 1.0); - if (x(0) >= 0.5 && x(1) >= 0.5) { return 1.0 * irg; } - if (x(0) < 0.5 && x(1) >= 0.5) { return 1.0 * irg; } - if (x(0) < 0.5 && x(1) < 0.5) { return 1.0 * irg; } - if (x(0) >= 0.5 && x(1) < 0.5) { return 1.0 * irg; } + const real_t irg = 1.0_r / rho0(x) / (gamma_func(x) - 1.0_r); + if (x(0) >= 0.5_r && x(1) >= 0.5_r) { return 1.0_r * irg; } + if (x(0) < 0.5_r && x(1) >= 0.5_r) { return 1.0_r * irg; } + if (x(0) < 0.5_r && x(1) < 0.5_r) { return 1.0_r * irg; } + if (x(0) >= 0.5_r && x(1) < 0.5_r) { return 1.0_r * irg; } MFEM_ABORT("Error in problem 6!"); - return 0.0; + return 0.0_r; } case 7: { - const double rho = rho0(x), gamma = gamma_func(x); - return (6.0 - rho * x(1)) / (gamma - 1.0) / rho; + const real_t rho = rho0(x), gamma = gamma_func(x); + return (6.0_r - rho * x(1)) / (gamma - 1.0_r) / rho; } - default: MFEM_ABORT("Bad number given for problem id!"); return 0.0; + default: MFEM_ABORT("Bad number given for problem id!"); return 0.0_r; } } @@ -1287,12 +1306,13 @@ static void display_banner(std::ostream &os) #ifdef LAGHOS_USE_CALIPER static void RecordAdiakMetadata(int dim, const char *mesh_file, int elem_per_mpi, - int nx, int ny, int nz, double blast_energy, - double Sx, double Sy, double Sz, + int nx, int ny, int nz, real_t blast_energy, + real_t Sx, real_t Sy, real_t Sz, int rs_levels, int rp_levels, int problem, int order_v, int order_e, int order_q, - int ode_solver_type, double t_final, double cfl, - double cg_tol, double ftz_tol, double delta_tol, + int ode_solver_type, real_t t_final, real_t cfl, + real_t cg_tol, real_t ftz_tol, + real_t delta_tol, int cg_max_iter, int max_tsteps, bool p_assembly, bool impose_visc, bool visualization, int vis_steps, bool visit, @@ -1309,9 +1329,9 @@ static void RecordAdiakMetadata(int dim, const char *mesh_file, int elem_per_mpi adiak::value("yelems", ny); adiak::value("zelems", nz); adiak::value("blast-energy", blast_energy); - adiak::value("xwidth", (double)Sx); - adiak::value("ywidth", (double)Sy); - adiak::value("zwidth", (double)Sz); + adiak::value("xwidth", Sx); + adiak::value("ywidth", Sy); + adiak::value("zwidth", Sz); adiak::value("refine-serial", rs_levels); adiak::value("refine-parallel", rp_levels); adiak::value("problem", problem); @@ -1414,17 +1434,17 @@ static bool ValidateElemPerMpiOptions(int elem_per_mpi, int argc, char *argv[], return false; } -static void Checks(const int ti, const double nrm, int &chk) +static void Checks(const int ti, const real_t nrm, int &chk) { - const double eps = 1.e-13; + const real_t eps = 1.e-13_r; - auto check = [&](int p, int i, const double res) + auto check = [&](int p, int i, const real_t res) { - auto rerr = [](const double a, const double v, const double eps) + auto rerr = [](const real_t a, const real_t v, const real_t eps) { MFEM_VERIFY(fabs(a) > eps && fabs(v) > eps, "One value is near zero!"); - const double err_a = fabs((a-v)/a); - const double err_v = fabs((a-v)/v); + const real_t err_a = fabs((a-v)/a); + const real_t err_v = fabs((a-v)/v); return fmax(err_a, err_v) < eps; }; if (problem == p && ti == i) @@ -1438,27 +1458,27 @@ static void Checks(const int ti, const double nrm, int &chk) } }; - const double it_norms[2][8][2][2] = // dim, problem, {it,norm} + const real_t it_norms[2][8][2][2] = // dim, problem, {it,norm} { { - {{5, 6.546538624534384e+00}, { 27, 7.588576357792927e+00}}, - {{5, 3.508254945225794e+00}, { 15, 2.756444596823211e+00}}, - {{5, 1.020745795651244e+01}, { 59, 1.721590205901898e+01}}, - {{5, 8.000000000000000e+00}, { 16, 8.000000000000000e+00}}, - {{5, 3.446324942352448e+01}, { 18, 3.446844033767240e+01}}, - {{5, 1.030899557252528e+01}, { 36, 1.057362418574309e+01}}, - {{5, 8.039707010835693e+00}, { 36, 8.316970976817373e+00}}, - {{5, 1.514929259650760e+01}, { 25, 1.514931278155159e+01}}, + {{5, 6.546538624534384e+00_r}, { 27, 7.588576357792927e+00_r}}, + {{5, 3.508254945225794e+00_r}, { 15, 2.756444596823211e+00_r}}, + {{5, 1.020745795651244e+01_r}, { 59, 1.721590205901898e+01_r}}, + {{5, 8.000000000000000e+00_r}, { 16, 8.000000000000000e+00_r}}, + {{5, 3.446324942352448e+01_r}, { 18, 3.446844033767240e+01_r}}, + {{5, 1.030899557252528e+01_r}, { 36, 1.057362418574309e+01_r}}, + {{5, 8.039707010835693e+00_r}, { 36, 8.316970976817373e+00_r}}, + {{5, 1.514929259650760e+01_r}, { 25, 1.514931278155159e+01_r}}, }, { - {{5, 1.198510951452527e+03}, {188, 1.199384410059154e+03}}, - {{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}}, - {{5, 2.061984481890964e+01}, { 36, 2.114519664792607e+01}}, - {{5, 1.607988713996459e+01}, { 36, 1.662736010353023e+01}}, - {{5, 3.029858112572883e+01}, { 24, 3.029858832743707e+01}} + {{5, 1.198510951452527e+03_r}, {188, 1.199384410059154e+03_r}}, + {{5, 6.695818592962833e+00_r}, { 20, 4.267902387082487e+00_r}}, + {{5, 2.041491591302486e+01_r}, { 59, 3.443180411803796e+01_r}}, + {{5, 1.600000000000000e+01_r}, { 16, 1.600000000000000e+01_r}}, + {{5, 6.892649884704898e+01_r}, { 18, 6.893688067534482e+01_r}}, + {{5, 2.061984481890964e+01_r}, { 36, 2.114519664792607e+01_r}}, + {{5, 1.607988713996459e+01_r}, { 36, 1.662736010353023e+01_r}}, + {{5, 3.029858112572883e+01_r}, { 24, 3.029858832743707e+01_r}} } }; @@ -1467,7 +1487,7 @@ static void Checks(const int ti, const double nrm, int &chk) for (int i=0; i<2; i++) { const int it = static_cast(it_norms[dim-2][p][i][0]); - const double norm = it_norms[dim-2][p][i][1]; + const real_t norm = it_norms[dim-2][p][i][1]; check(p, it, norm); } } @@ -1476,12 +1496,12 @@ 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; + constexpr real_t tol = 1e-6_r; const int NBE = mesh.GetNBE(); IntegrationPoint center; - center.x = 0.5; - center.y = 0.5; - center.z = 0.5; + center.x = 0.5_r; + center.y = 0.5_r; + center.z = 0.5_r; for (int b = 0; b < NBE; b++) { Element *bel = mesh.GetBdrElement(b); @@ -1500,12 +1520,12 @@ 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; + constexpr real_t tol = 1e-6_r; const int NBE = mesh.GetNBE(); IntegrationPoint center; - center.x = 0.5; - center.y = 0.5; - center.z = 0.5; + center.x = 0.5_r; + center.y = 0.5_r; + center.z = 0.5_r; for (int b = 0; b < NBE; b++) { Element *bel = mesh.GetBdrElement(b); diff --git a/laghos_assembly.cpp b/laghos_assembly.cpp index fab72aab..e5cd0f87 100644 --- a/laghos_assembly.cpp +++ b/laghos_assembly.cpp @@ -30,7 +30,7 @@ void DensityIntegrator::AssembleRHSElementVect(const FiniteElement &fe, const int nqp = IntRule->GetNPoints(); Vector shape(fe.GetDof()); elvect.SetSize(fe.GetDof()); - elvect = 0.0; + elvect = 0.0_r; for (int q = 0; q < nqp; q++) { fe.CalcShape(IntRule->IntPoint(q), shape); @@ -51,7 +51,7 @@ void ForceIntegrator::AssembleElementMatrix2(const FiniteElement &trial_fe, const int h1dofs_cnt = test_fe.GetDof(); const int l2dofs_cnt = trial_fe.GetDof(); elmat.SetSize(h1dofs_cnt*dim, l2dofs_cnt); - elmat = 0.0; + elmat = 0.0_r; DenseMatrix vshape(h1dofs_cnt, dim), loc_force(h1dofs_cnt, dim); Vector shape(l2dofs_cnt), Vloc_force(loc_force.Data(), h1dofs_cnt*dim); for (int q = 0; q < nqp; q++) @@ -63,11 +63,11 @@ void ForceIntegrator::AssembleElementMatrix2(const FiniteElement &trial_fe, { for (int vd = 0; vd < dim; vd++) // Velocity components. { - loc_force(i, vd) = 0.0; + loc_force(i, vd) = 0.0_r; for (int gd = 0; gd < dim; gd++) // Gradient components. { const int eq = e*nqp + q; - const double stressJinvT = qdata.stressJinvT(vd)(eq, gd); + const real_t stressJinvT = qdata.stressJinvT(vd)(eq, gd); loc_force(i, vd) += stressJinvT * vshape(i,gd); } } @@ -111,13 +111,13 @@ void MassPAOperator::SetEssentialTrueDofs(Array &dofs) void MassPAOperator::EliminateRHS(Vector &b) const { - if (ess_tdofs_count > 0) { b.SetSubVector(ess_tdofs, 0.0); } + if (ess_tdofs_count > 0) { b.SetSubVector(ess_tdofs, 0.0_r); } } void MassPAOperator::Mult(const Vector &x, Vector &y) const { mass->Mult(x, y); - if (ess_tdofs_count > 0) { y.SetSubVector(ess_tdofs, 0.0); } + if (ess_tdofs_count > 0) { y.SetSubVector(ess_tdofs, 0.0_r); } } ForcePAOperator::ForcePAOperator(const QuadratureData &qdata, @@ -144,41 +144,41 @@ ForcePAOperator::ForcePAOperator(const QuadratureData &qdata, template static void ForceMult2D(const int NE, - const Array &B_, - const Array &Bt_, - const Array &Gt_, + const Array &B_, + const Array &Bt_, + const Array &Gt_, const DenseTensor &sJit_, const Vector &x, Vector &y) { auto b = Reshape(B_.Read(), Q1D, L1D); auto bt = Reshape(Bt_.Read(), D1D, Q1D); auto gt = Reshape(Gt_.Read(), D1D, Q1D); - const double *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*NE*DIM*DIM); + const real_t *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*NE*DIM*DIM); auto sJit = Reshape(StressJinvT, Q1D, Q1D, NE, DIM, DIM); auto energy = Reshape(x.Read(), L1D, L1D, NE); - const double eps1 = std::numeric_limits::epsilon(); - const double eps2 = eps1*eps1; + const real_t eps1 = std::numeric_limits::epsilon(); + const real_t eps2 = eps1*eps1; auto velocity = Reshape(y.Write(), D1D, D1D, DIM, NE); MFEM_FORALL_2D(e, NE, Q1D, Q1D, 1, { const int z = MFEM_THREAD_ID(z); - MFEM_SHARED double B[Q1D][L1D]; - MFEM_SHARED double Bt[D1D][Q1D]; - MFEM_SHARED double Gt[D1D][Q1D]; + MFEM_SHARED real_t B[Q1D][L1D]; + MFEM_SHARED real_t Bt[D1D][Q1D]; + MFEM_SHARED real_t Gt[D1D][Q1D]; - MFEM_SHARED double Ez[NBZ][L1D][L1D]; - double (*E)[L1D] = (double (*)[L1D])(Ez + z); + MFEM_SHARED real_t Ez[NBZ][L1D][L1D]; + real_t (*E)[L1D] = (real_t (*)[L1D])(Ez + z); - MFEM_SHARED double LQz[2][NBZ][D1D][Q1D]; - double (*LQ0)[Q1D] = (double (*)[Q1D])(LQz[0] + z); - double (*LQ1)[Q1D] = (double (*)[Q1D])(LQz[1] + z); + MFEM_SHARED real_t LQz[2][NBZ][D1D][Q1D]; + real_t (*LQ0)[Q1D] = (real_t (*)[Q1D])(LQz[0] + z); + real_t (*LQ1)[Q1D] = (real_t (*)[Q1D])(LQz[1] + z); - MFEM_SHARED double QQz[3][NBZ][Q1D][Q1D]; - double (*QQ)[Q1D] = (double (*)[Q1D])(QQz[0] + z); - double (*QQ0)[Q1D] = (double (*)[Q1D])(QQz[1] + z); - double (*QQ1)[Q1D] = (double (*)[Q1D])(QQz[2] + z); + MFEM_SHARED real_t QQz[3][NBZ][Q1D][Q1D]; + real_t (*QQ)[Q1D] = (real_t (*)[Q1D])(QQz[0] + z); + real_t (*QQ0)[Q1D] = (real_t (*)[Q1D])(QQz[1] + z); + real_t (*QQ1)[Q1D] = (real_t (*)[Q1D])(QQz[2] + z); if (z == 0) { @@ -207,7 +207,7 @@ void ForceMult2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; + real_t u = 0.0_r; for (int lx = 0; lx < L1D; ++lx) { u += B[qx][lx] * E[lx][ly]; @@ -220,7 +220,7 @@ void ForceMult2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; + real_t u = 0.0_r; for (int ly = 0; ly < L1D; ++ly) { u += B[qy][ly] * LQ0[ly][qx]; @@ -236,8 +236,8 @@ void ForceMult2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - const double esx = QQ[qy][qx] * sJit(qx,qy,e,0,c); - const double esy = QQ[qy][qx] * sJit(qx,qy,e,1,c); + const real_t esx = QQ[qy][qx] * sJit(qx,qy,e,0,c); + const real_t esy = QQ[qy][qx] * sJit(qx,qy,e,1,c); QQ0[qy][qx] = esx; QQ1[qy][qx] = esy; } @@ -247,8 +247,8 @@ void ForceMult2D(const int NE, { MFEM_FOREACH_THREAD(dx,x,D1D) { - double u = 0.0; - double v = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; for (int qx = 0; qx < Q1D; ++qx) { u += Gt[dx][qx] * QQ0[qy][qx]; @@ -263,8 +263,8 @@ void ForceMult2D(const int NE, { MFEM_FOREACH_THREAD(dx,x,D1D) { - double u = 0.0; - double v = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; for (int qy = 0; qy < Q1D; ++qy) { u += LQ0[dx][qy] * Bt[dy][qy]; @@ -281,10 +281,10 @@ void ForceMult2D(const int NE, { MFEM_FOREACH_THREAD(dx,x,D1D) { - const double v = velocity(dx,dy,c,e); + const real_t v = velocity(dx,dy,c,e); if (fabs(v) < eps2) { - velocity(dx,dy,c,e) = 0.0; + velocity(dx,dy,c,e) = 0.0_r; } } } @@ -295,47 +295,47 @@ void ForceMult2D(const int NE, template static void ForceMult3D(const int NE, - const Array &B_, - const Array &Bt_, - const Array &Gt_, + const Array &B_, + const Array &Bt_, + const Array &Gt_, const DenseTensor &sJit_, const Vector &x, Vector &y) { auto b = Reshape(B_.Read(), Q1D, L1D); auto bt = Reshape(Bt_.Read(), D1D, Q1D); auto gt = Reshape(Gt_.Read(), D1D, Q1D); - const double *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*Q1D*NE*DIM*DIM); + const real_t *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*Q1D*NE*DIM*DIM); auto sJit = Reshape(StressJinvT, Q1D, Q1D, Q1D, NE, DIM, DIM); auto energy = Reshape(x.Read(), L1D, L1D, L1D, NE); - const double eps1 = std::numeric_limits::epsilon(); - const double eps2 = eps1*eps1; + const real_t eps1 = std::numeric_limits::epsilon(); + const real_t eps2 = eps1*eps1; auto velocity = Reshape(y.Write(), D1D, D1D, D1D, DIM, NE); MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D, { const int z = MFEM_THREAD_ID(z); - MFEM_SHARED double B[Q1D][L1D]; - MFEM_SHARED double Bt[D1D][Q1D]; - MFEM_SHARED double Gt[D1D][Q1D]; + MFEM_SHARED real_t B[Q1D][L1D]; + MFEM_SHARED real_t Bt[D1D][Q1D]; + MFEM_SHARED real_t Gt[D1D][Q1D]; - MFEM_SHARED double E[L1D][L1D][L1D]; + MFEM_SHARED real_t E[L1D][L1D][L1D]; - MFEM_SHARED double sm0[3][Q1D*Q1D*Q1D]; - MFEM_SHARED double sm1[3][Q1D*Q1D*Q1D]; + MFEM_SHARED real_t sm0[3][Q1D*Q1D*Q1D]; + MFEM_SHARED real_t sm1[3][Q1D*Q1D*Q1D]; - double (*MMQ0)[D1D][Q1D] = (double (*)[D1D][Q1D]) (sm0+0); - double (*MMQ1)[D1D][Q1D] = (double (*)[D1D][Q1D]) (sm0+1); - double (*MMQ2)[D1D][Q1D] = (double (*)[D1D][Q1D]) (sm0+2); + real_t (*MMQ0)[D1D][Q1D] = (real_t (*)[D1D][Q1D]) (sm0+0); + real_t (*MMQ1)[D1D][Q1D] = (real_t (*)[D1D][Q1D]) (sm0+1); + real_t (*MMQ2)[D1D][Q1D] = (real_t (*)[D1D][Q1D]) (sm0+2); - double (*MQQ0)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm1+0); - double (*MQQ1)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm1+1); - double (*MQQ2)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm1+2); + real_t (*MQQ0)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm1+0); + real_t (*MQQ1)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm1+1); + real_t (*MQQ2)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm1+2); - MFEM_SHARED double QQQ[Q1D][Q1D][Q1D]; - double (*QQQ0)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm0+0); - double (*QQQ1)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm0+1); - double (*QQQ2)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm0+2); + MFEM_SHARED real_t QQQ[Q1D][Q1D][Q1D]; + real_t (*QQQ0)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm0+0); + real_t (*QQQ1)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm0+1); + real_t (*QQQ2)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm0+2); if (z == 0) { @@ -367,7 +367,7 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; + real_t u = 0.0_r; for (int lx = 0; lx < L1D; ++lx) { u += B[qx][lx] * E[lx][ly][lz]; @@ -383,7 +383,7 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; + real_t u = 0.0_r; for (int ly = 0; ly < L1D; ++ly) { u += B[qy][ly] * MMQ0[lz][ly][qx]; @@ -399,7 +399,7 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; + real_t u = 0.0_r; for (int lz = 0; lz < L1D; ++lz) { u += B[qz][lz] * MQQ0[lz][qy][qx]; @@ -417,9 +417,9 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - const double esx = QQQ[qz][qy][qx] * sJit(qx,qy,qz,e,0,c); - const double esy = QQQ[qz][qy][qx] * sJit(qx,qy,qz,e,1,c); - const double esz = QQQ[qz][qy][qx] * sJit(qx,qy,qz,e,2,c); + const real_t esx = QQQ[qz][qy][qx] * sJit(qx,qy,qz,e,0,c); + const real_t esy = QQQ[qz][qy][qx] * sJit(qx,qy,qz,e,1,c); + const real_t esz = QQQ[qz][qy][qx] * sJit(qx,qy,qz,e,2,c); QQQ0[qz][qy][qx] = esx; QQQ1[qz][qy][qx] = esy; QQQ2[qz][qy][qx] = esz; @@ -433,9 +433,9 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(hx,x,D1D) { - double u = 0.0; - double v = 0.0; - double w = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; + real_t w = 0.0_r; for (int qx = 0; qx < Q1D; ++qx) { u += Gt[hx][qx] * QQQ0[qz][qy][qx]; @@ -455,9 +455,9 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(hx,x,D1D) { - double u = 0.0; - double v = 0.0; - double w = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; + real_t w = 0.0_r; for (int qy = 0; qy < Q1D; ++qy) { u += MQQ0[hx][qy][qz] * Bt[hy][qy]; @@ -477,9 +477,9 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(hx,x,D1D) { - double u = 0.0; - double v = 0.0; - double w = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; + real_t w = 0.0_r; for (int qz = 0; qz < Q1D; ++qz) { u += MMQ0[hx][hy][qz] * Bt[hz][qz]; @@ -500,10 +500,10 @@ void ForceMult3D(const int NE, { MFEM_FOREACH_THREAD(hx,x,D1D) { - const double v = velocity(hx,hy,hz,c,e); + const real_t v = velocity(hx,hy,hz,c,e); if (fabs(v) < eps2) { - velocity(hx,hy,hz,c,e) = 0.0; + velocity(hx,hy,hz,c,e) = 0.0_r; } } } @@ -514,17 +514,17 @@ void ForceMult3D(const int NE, } typedef void (*fForceMult)(const int E, - const Array &B, - const Array &Bt, - const Array &Gt, + const Array &B, + const Array &Bt, + const Array &Gt, const DenseTensor &stressJinvT, const Vector &X, Vector &Y); static void ForceMult(const int DIM, const int D1D, const int Q1D, const int L1D, const int H1D, const int NE, - const Array &B, - const Array &Bt, - const Array &Gt, + const Array &B, + const Array &Bt, + const Array &Gt, const DenseTensor &stressJinvT, const Vector &e, Vector &v) @@ -566,16 +566,16 @@ void ForcePAOperator::Mult(const Vector &x, Vector &y) const template static void ForceMultTranspose2D(const int NE, - const Array &Bt_, - const Array &B_, - const Array &G_, + const Array &Bt_, + const Array &B_, + const Array &G_, const DenseTensor &sJit_, const Vector &x, Vector &y) { auto b = Reshape(B_.Read(), Q1D, D1D); auto g = Reshape(G_.Read(), Q1D, D1D); auto bt = Reshape(Bt_.Read(), L1D, Q1D); - const double *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*NE*DIM*DIM); + const real_t *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*NE*DIM*DIM); auto sJit = Reshape(StressJinvT, Q1D, Q1D, NE, DIM, DIM); auto velocity = Reshape(x.Read(), D1D, D1D, DIM, NE); auto energy = Reshape(y.Write(), L1D, L1D, NE); @@ -584,24 +584,24 @@ void ForceMultTranspose2D(const int NE, { const int z = MFEM_THREAD_ID(z); - MFEM_SHARED double Bt[L1D][Q1D]; - MFEM_SHARED double B[Q1D][D1D]; - MFEM_SHARED double G[Q1D][D1D]; + MFEM_SHARED real_t Bt[L1D][Q1D]; + MFEM_SHARED real_t B[Q1D][D1D]; + MFEM_SHARED real_t G[Q1D][D1D]; - MFEM_SHARED double Vz[NBZ][D1D*D1D]; - double (*V)[D1D] = (double (*)[D1D])(Vz + z); + MFEM_SHARED real_t Vz[NBZ][D1D*D1D]; + real_t (*V)[D1D] = (real_t (*)[D1D])(Vz + z); - MFEM_SHARED double DQz[DIM][NBZ][D1D*Q1D]; - double (*DQ0)[Q1D] = (double (*)[Q1D])(DQz[0] + z); - double (*DQ1)[Q1D] = (double (*)[Q1D])(DQz[1] + z); + MFEM_SHARED real_t DQz[DIM][NBZ][D1D*Q1D]; + real_t (*DQ0)[Q1D] = (real_t (*)[Q1D])(DQz[0] + z); + real_t (*DQ1)[Q1D] = (real_t (*)[Q1D])(DQz[1] + z); - MFEM_SHARED double QQz[3][NBZ][Q1D*Q1D]; - double (*QQ)[Q1D] = (double (*)[Q1D])(QQz[0] + z); - double (*QQ0)[Q1D] = (double (*)[Q1D])(QQz[1] + z); - double (*QQ1)[Q1D] = (double (*)[Q1D])(QQz[2] + z); + MFEM_SHARED real_t QQz[3][NBZ][Q1D*Q1D]; + real_t (*QQ)[Q1D] = (real_t (*)[Q1D])(QQz[0] + z); + real_t (*QQ0)[Q1D] = (real_t (*)[Q1D])(QQz[1] + z); + real_t (*QQ1)[Q1D] = (real_t (*)[Q1D])(QQz[2] + z); - MFEM_SHARED double QLz[NBZ][Q1D*L1D]; - double (*QL)[L1D] = (double (*)[L1D]) (QLz + z); + MFEM_SHARED real_t QLz[NBZ][Q1D*L1D]; + real_t (*QL)[L1D] = (real_t (*)[L1D]) (QLz + z); if (z == 0) { @@ -621,7 +621,7 @@ void ForceMultTranspose2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - QQ[qy][qx] = 0.0; + QQ[qy][qx] = 0.0_r; } } MFEM_SYNC_THREAD; @@ -641,11 +641,11 @@ void ForceMultTranspose2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; - double v = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; for (int dx = 0; dx < D1D; ++dx) { - const double input = V[dx][dy]; + const real_t input = V[dx][dy]; u += B[qx][dx] * input; v += G[qx][dx] * input; } @@ -658,8 +658,8 @@ void ForceMultTranspose2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; - double v = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; for (int dy = 0; dy < D1D; ++dy) { u += DQ1[dy][qx] * B[qy][dy]; @@ -674,8 +674,8 @@ void ForceMultTranspose2D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - const double esx = QQ0[qy][qx] * sJit(qx,qy,e,0,c); - const double esy = QQ1[qy][qx] * sJit(qx,qy,e,1,c); + const real_t esx = QQ0[qy][qx] * sJit(qx,qy,e,0,c); + const real_t esy = QQ1[qy][qx] * sJit(qx,qy,e,1,c); QQ[qy][qx] += esx + esy; } } @@ -687,7 +687,7 @@ void ForceMultTranspose2D(const int NE, { MFEM_FOREACH_THREAD(lx,x,L1D) { - double u = 0.0; + real_t u = 0.0_r; for (int qx = 0; qx < Q1D; ++qx) { u += QQ[qy][qx] * Bt[lx][qx]; @@ -700,7 +700,7 @@ void ForceMultTranspose2D(const int NE, { MFEM_FOREACH_THREAD(lx,x,L1D) { - double u = 0.0; + real_t u = 0.0_r; for (int qy = 0; qy < Q1D; ++qy) { u += QL[qy][lx] * Bt[ly][qy]; @@ -714,9 +714,9 @@ void ForceMultTranspose2D(const int NE, template static void ForceMultTranspose3D(const int NE, - const Array &Bt_, - const Array &B_, - const Array &G_, + const Array &Bt_, + const Array &B_, + const Array &G_, const DenseTensor &sJit_, const Vector &v_, Vector &e_) @@ -724,7 +724,7 @@ void ForceMultTranspose3D(const int NE, auto b = Reshape(B_.Read(), Q1D, D1D); auto g = Reshape(G_.Read(), Q1D, D1D); auto bt = Reshape(Bt_.Read(), L1D, Q1D); - const double *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*Q1D*NE*DIM*DIM); + const real_t *StressJinvT = Read(sJit_.GetMemory(), Q1D*Q1D*Q1D*NE*DIM*DIM); auto sJit = Reshape(StressJinvT, Q1D, Q1D, Q1D, NE, DIM, DIM); auto velocity = Reshape(v_.Read(), D1D, D1D, D1D, DIM, NE); auto energy = Reshape(e_.Write(), L1D, L1D, L1D, NE); @@ -733,25 +733,25 @@ void ForceMultTranspose3D(const int NE, { const int z = MFEM_THREAD_ID(z); - MFEM_SHARED double Bt[L1D][Q1D]; - MFEM_SHARED double B[Q1D][D1D]; - MFEM_SHARED double G[Q1D][D1D]; + MFEM_SHARED real_t Bt[L1D][Q1D]; + MFEM_SHARED real_t B[Q1D][D1D]; + MFEM_SHARED real_t G[Q1D][D1D]; - MFEM_SHARED double sm0[3][Q1D*Q1D*Q1D]; - MFEM_SHARED double sm1[3][Q1D*Q1D*Q1D]; - double (*V)[D1D][D1D] = (double (*)[D1D][D1D]) (sm0+0); - double (*MMQ0)[D1D][Q1D] = (double (*)[D1D][Q1D]) (sm0+1); - double (*MMQ1)[D1D][Q1D] = (double (*)[D1D][Q1D]) (sm0+2); + MFEM_SHARED real_t sm0[3][Q1D*Q1D*Q1D]; + MFEM_SHARED real_t sm1[3][Q1D*Q1D*Q1D]; + real_t (*V)[D1D][D1D] = (real_t (*)[D1D][D1D]) (sm0+0); + real_t (*MMQ0)[D1D][Q1D] = (real_t (*)[D1D][Q1D]) (sm0+1); + real_t (*MMQ1)[D1D][Q1D] = (real_t (*)[D1D][Q1D]) (sm0+2); - double (*MQQ0)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm1+0); - double (*MQQ1)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm1+1); - double (*MQQ2)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm1+2); + real_t (*MQQ0)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm1+0); + real_t (*MQQ1)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm1+1); + real_t (*MQQ2)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm1+2); - double (*QQQ0)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm0+0); - double (*QQQ1)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm0+1); - double (*QQQ2)[Q1D][Q1D] = (double (*)[Q1D][Q1D]) (sm0+2); + real_t (*QQQ0)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm0+0); + real_t (*QQQ1)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm0+1); + real_t (*QQQ2)[Q1D][Q1D] = (real_t (*)[Q1D][Q1D]) (sm0+2); - MFEM_SHARED double QQQ[Q1D][Q1D][Q1D]; + MFEM_SHARED real_t QQQ[Q1D][Q1D][Q1D]; if (z == 0) { @@ -773,7 +773,7 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - QQQ[qz][qy][qx] = 0.0; + QQQ[qz][qy][qx] = 0.0_r; } } } @@ -798,11 +798,11 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; - double v = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; for (int dx = 0; dx < D1D; ++dx) { - const double input = V[dx][dy][dz]; + const real_t input = V[dx][dy][dz]; u += G[qx][dx] * input; v += B[qx][dx] * input; } @@ -818,9 +818,9 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; - double v = 0.0; - double w = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; + real_t w = 0.0_r; for (int dy = 0; dy < D1D; ++dy) { u += MMQ0[dz][dy][qx] * B[qy][dy]; @@ -840,9 +840,9 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - double u = 0.0; - double v = 0.0; - double w = 0.0; + real_t u = 0.0_r; + real_t v = 0.0_r; + real_t w = 0.0_r; for (int dz = 0; dz < D1D; ++dz) { u += MQQ0[dz][qy][qx] * B[qz][dz]; @@ -862,9 +862,9 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(qx,x,Q1D) { - const double esx = QQQ0[qz][qy][qx] * sJit(qx,qy,qz,e,0,c); - const double esy = QQQ1[qz][qy][qx] * sJit(qx,qy,qz,e,1,c); - const double esz = QQQ2[qz][qy][qx] * sJit(qx,qy,qz,e,2,c); + const real_t esx = QQQ0[qz][qy][qx] * sJit(qx,qy,qz,e,0,c); + const real_t esy = QQQ1[qz][qy][qx] * sJit(qx,qy,qz,e,1,c); + const real_t esz = QQQ2[qz][qy][qx] * sJit(qx,qy,qz,e,2,c); QQQ[qz][qy][qx] += esx + esy + esz; } } @@ -878,7 +878,7 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(lx,x,L1D) { - double u = 0.0; + real_t u = 0.0_r; for (int qx = 0; qx < Q1D; ++qx) { u += QQQ[qz][qy][qx] * Bt[lx][qx]; @@ -894,7 +894,7 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(lx,x,L1D) { - double u = 0.0; + real_t u = 0.0_r; for (int qy = 0; qy < Q1D; ++qy) { u += MQQ0[qz][qy][lx] * Bt[ly][qy]; @@ -910,7 +910,7 @@ void ForceMultTranspose3D(const int NE, { MFEM_FOREACH_THREAD(lx,x,L1D) { - double u = 0.0; + real_t u = 0.0_r; for (int qz = 0; qz < Q1D; ++qz) { u += MMQ0[qz][ly][lx] * Bt[lz][qz]; @@ -924,17 +924,17 @@ void ForceMultTranspose3D(const int NE, } typedef void (*fForceMultTranspose)(const int NE, - const Array &Bt, - const Array &B, - const Array &G, + const Array &Bt, + const Array &B, + const Array &G, const DenseTensor &sJit, const Vector &X, Vector &Y); static void ForceMultTranspose(const int DIM, const int D1D, const int Q1D, const int L1D, const int NE, - const Array &L2Bt, - const Array &H1B, - const Array &H1G, + const Array &L2Bt, + const Array &H1B, + const Array &H1G, const DenseTensor &stressJinvT, const Vector &v, Vector &e) diff --git a/laghos_assembly.hpp b/laghos_assembly.hpp index 0699352e..c0e87df8 100644 --- a/laghos_assembly.hpp +++ b/laghos_assembly.hpp @@ -49,11 +49,11 @@ struct QuadratureData // Initial length scale. This represents a notion of local mesh size. // We assume that all initial zones have similar size. - double h0; + real_t h0; // Estimate of the minimum time step over all quadrature points. This is // recomputed at every time step to achieve adaptive time stepping. - double dt_est; + real_t dt_est; QuadratureData(int dim, int NE, int quads_per_el) : Jac0inv(dim, dim, NE * quads_per_el), diff --git a/laghos_solver.cpp b/laghos_solver.cpp index a6f2ded0..9f6997a9 100644 --- a/laghos_solver.cpp +++ b/laghos_solver.cpp @@ -99,7 +99,7 @@ static void Rho0DetJ0Vol(const int dim, const int NE, ParFiniteElementSpace &L2, const ParGridFunction &rho0, QuadratureData &qdata, - double &volume); + real_t &volume); LagrangianHydroOperator::LagrangianHydroOperator(const int size, ParFiniteElementSpace &h1, @@ -109,13 +109,13 @@ LagrangianHydroOperator::LagrangianHydroOperator(const int size, ParGridFunction &rho0_gf, ParGridFunction &gamma_gf, const int source, - const double cfl, + const real_t cfl, const bool visc, const bool vort, const bool p_assembly, - const double cgt, + const real_t cgt, const int cgiter, - double ftz, + real_t ftz, const int oq) : TimeDependentOperator(size), H1(h1), L2(l2), H1c(H1.GetParMesh(), H1.FEColl(), 1), @@ -144,7 +144,7 @@ LagrangianHydroOperator::LagrangianHydroOperator(const int size, Me_inv(l2dofs_cnt, l2dofs_cnt, NE), ir(IntRules.Get(pmesh->GetElementBaseGeometry(0), (oq > 0) ? oq : 3 * H1.GetOrder(0) + L2.GetOrder(0) - 1)), - Q1D(int(floor(0.7 + pow(ir.GetNPoints(), 1.0 / dim)))), + Q1D(int(floor(0.7_r + pow(ir.GetNPoints(), 1.0_r / dim)))), qdata(dim, NE, ir.GetNPoints()), qdata_is_current(false), forcemat_is_assembled(false), @@ -168,7 +168,7 @@ LagrangianHydroOperator::LagrangianHydroOperator(const int size, block_offsets[2] = block_offsets[1] + H1Vsize; block_offsets[3] = block_offsets[2] + L2Vsize; one.UseDevice(true); - one = 1.0; + one = 1.0_r; if (p_assembly) { @@ -223,7 +223,7 @@ LagrangianHydroOperator::LagrangianHydroOperator(const int size, // Values of rho0DetJ0 and Jac0inv at all quadrature points. // Initial local mesh size (assumes all mesh elements are the same). HYPRE_BigInt Ne, ne = NE; - double Volume, vol = 0.0; + real_t Volume, vol = 0.0_r; if (dim > 1 && p_assembly) { Rho0DetJ0Vol(dim, NE, ir, pmesh, L2, rho0_gf, qdata, vol); @@ -242,24 +242,26 @@ LagrangianHydroOperator::LagrangianHydroOperator(const int size, Tr.SetIntPoint(&ip); DenseMatrixInverse Jinv(Tr.Jacobian()); Jinv.GetInverseMatrix(qdata.Jac0inv(e*NQ + q)); - const double rho0DetJ0 = Tr.Weight() * rho_vals(q); + const real_t rho0DetJ0 = Tr.Weight() * rho_vals(q); qdata.rho0DetJ0w(e*NQ + q) = rho0DetJ0 * ir.IntPoint(q).weight; } } for (int e = 0; e < NE; e++) { vol += pmesh->GetElementVolume(e); } } - MPI_Allreduce(&vol, &Volume, 1, MPI_DOUBLE, MPI_SUM, pmesh->GetComm()); + MPI_Allreduce(&vol, &Volume, 1, MPITypeMap::mpi_type, MPI_SUM, + pmesh->GetComm()); MPI_Allreduce(&ne, &Ne, 1, HYPRE_MPI_BIG_INT, MPI_SUM, pmesh->GetComm()); switch (pmesh->GetElementBaseGeometry(0)) { case Geometry::SEGMENT: qdata.h0 = Volume / Ne; break; case Geometry::SQUARE: qdata.h0 = sqrt(Volume / Ne); break; - case Geometry::TRIANGLE: qdata.h0 = sqrt(2.0 * Volume / Ne); break; - case Geometry::CUBE: qdata.h0 = pow(Volume / Ne, 1./3.); break; - case Geometry::TETRAHEDRON: qdata.h0 = pow(6.0 * Volume / Ne, 1./3.); break; + case Geometry::TRIANGLE: qdata.h0 = sqrt(2.0_r * Volume / Ne); break; + case Geometry::CUBE: qdata.h0 = pow(Volume / Ne, 1.0_r/3.0_r); break; + case Geometry::TETRAHEDRON: + qdata.h0 = pow(6.0_r * Volume / Ne, 1.0_r/3.0_r); break; default: MFEM_ABORT("Unknown zone type!"); } - qdata.h0 /= (double) H1.GetOrder(0); + qdata.h0 /= static_cast(H1.GetOrder(0)); if (p_assembly) { @@ -271,14 +273,14 @@ LagrangianHydroOperator::LagrangianHydroOperator(const int size, CG_VMass.SetOperator(*VMassPA); CG_VMass.SetRelTol(cg_rel_tol); - CG_VMass.SetAbsTol(0.0); + CG_VMass.SetAbsTol(0.0_r); CG_VMass.SetMaxIter(cg_max_iter); CG_VMass.SetPrintLevel(-1); CG_EMass.SetOperator(*EMassPA); CG_EMass.iterative_mode = false; CG_EMass.SetRelTol(cg_rel_tol); - CG_EMass.SetAbsTol(0.0); + CG_EMass.SetAbsTol(0.0_r); CG_EMass.SetMaxIter(cg_max_iter); CG_EMass.SetPrintLevel(-1); } @@ -335,7 +337,7 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S, // (Position, Velocity, Specific Internal Energy). ParGridFunction dv; dv.MakeRef(&H1, dS_dt, H1Vsize); - dv = 0.0; + dv = 0.0_r; ParGridFunction accel_src_gf; if (source_type == 2) @@ -424,7 +426,7 @@ void LagrangianHydroOperator::SolveVelocity(const Vector &S, cg.SetPreconditioner(prec); cg.SetOperator(A); cg.SetRelTol(cg_rel_tol); - cg.SetAbsTol(0.0); + cg.SetAbsTol(0.0_r); cg.SetMaxIter(cg_max_iter); cg.SetPrintLevel(-1); LAGHOS_DEVICE_SYNC; @@ -449,7 +451,7 @@ void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v, // (Position, Velocity, Specific Internal Energy). ParGridFunction de; de.MakeRef(&L2, dS_dt, H1Vsize*2); - de = 0.0; + de = 0.0_r; // Solve for energy, assemble the energy source if such exists. LinearForm *e_source = nullptr; @@ -524,19 +526,20 @@ void LagrangianHydroOperator::UpdateMesh(const Vector &S) const H1.GetParMesh()->NewNodes(x_gf, false); } -double LagrangianHydroOperator::GetTimeStepEstimate(const Vector &S) const +real_t LagrangianHydroOperator::GetTimeStepEstimate(const Vector &S) const { UpdateMesh(S); UpdateQuadratureData(S); - double glob_dt_est; + real_t glob_dt_est; const MPI_Comm comm = H1.GetParMesh()->GetComm(); - MPI_Allreduce(&qdata.dt_est, &glob_dt_est, 1, MPI_DOUBLE, MPI_MIN, comm); + MPI_Allreduce(&qdata.dt_est, &glob_dt_est, 1, MPITypeMap::mpi_type, + MPI_MIN, comm); return glob_dt_est; } void LagrangianHydroOperator::ResetTimeStepEstimate() const { - qdata.dt_est = std::numeric_limits::infinity(); + qdata.dt_est = std::numeric_limits::infinity(); } void LagrangianHydroOperator::ComputeDensity(ParGridFunction &rho) const @@ -562,9 +565,9 @@ void LagrangianHydroOperator::ComputeDensity(ParGridFunction &rho) const } } -double ComputeVolumeIntegral(const ParFiniteElementSpace &pfes, +real_t ComputeVolumeIntegral(const ParFiniteElementSpace &pfes, const int DIM, const int NE, const int NQ, - const int Q1D, const int VDIM, const double norm, + const int Q1D, const int VDIM, const real_t norm, const Vector& mass, const Vector& f) { MFEM_VERIFY(pfes.GetNE() > 0, "Empty local mesh should have been handled!"); @@ -584,7 +587,7 @@ double ComputeVolumeIntegral(const ParFiniteElementSpace &pfes, { for (int q = 0; q < NQ; ++q) { - double vmag = 0; + real_t vmag = 0.0_r; for (int k = 0; k < VDIM; k++) { vmag += pow(f_vals(k,q,e), norm); @@ -602,7 +605,7 @@ double ComputeVolumeIntegral(const ParFiniteElementSpace &pfes, MFEM_FOREACH_THREAD(qx,x,QX) { const int q = qx + qy * QX; - double vmag = 0.0; + real_t vmag = 0.0_r; for (int k = 0; k < VDIM; k++) { vmag += pow(f_vals(k, q, e), norm); @@ -623,7 +626,7 @@ double ComputeVolumeIntegral(const ParFiniteElementSpace &pfes, MFEM_FOREACH_THREAD(qx,x,QX) { const int q = qx + (qy + qz * QY) * QX; - double vmag = 0; + real_t vmag = 0.0_r; for (int k = 0; k < VDIM; k++) { vmag += pow(f_vals(k, q, e), norm); @@ -637,9 +640,9 @@ double ComputeVolumeIntegral(const ParFiniteElementSpace &pfes, return integrand * mass; } -double LagrangianHydroOperator::InternalEnergy(const ParGridFunction &gf) const +real_t LagrangianHydroOperator::InternalEnergy(const ParGridFunction &gf) const { - double glob_ie = 0.0, internal_energy = 0.0; + real_t glob_ie = 0.0_r, internal_energy = 0.0_r; if (L2.GetNE() > 0) // UsesTensorBasis does not handle empty local mesh { @@ -657,18 +660,20 @@ double LagrangianHydroOperator::InternalEnergy(const ParGridFunction &gf) const L2r->Mult(gf, e_vec); L2qi->Values(e_vec, q_val); internal_energy = - ComputeVolumeIntegral(L2, dim, NE, NQ, Q1D, 1, 1.0, qdata.rho0DetJ0w, q_val); + ComputeVolumeIntegral(L2, dim, NE, NQ, Q1D, 1, 1.0_r, + qdata.rho0DetJ0w, q_val); } - MPI_Allreduce(&internal_energy, &glob_ie, 1, MPI_DOUBLE, MPI_SUM, + MPI_Allreduce(&internal_energy, &glob_ie, 1, MPITypeMap::mpi_type, + MPI_SUM, L2.GetParMesh()->GetComm()); return glob_ie; } -double LagrangianHydroOperator::KineticEnergy(const ParGridFunction &v) const +real_t LagrangianHydroOperator::KineticEnergy(const ParGridFunction &v) const { - double glob_ke = 0.0, kinetic_energy = 0.0; + real_t glob_ke = 0.0_r, kinetic_energy = 0.0_r; if (H1.GetNE() > 0) // UsesTensorBasis does not handle empty local mesh { @@ -687,13 +692,15 @@ double LagrangianHydroOperator::KineticEnergy(const ParGridFunction &v) const h1_interpolator->Values(e_vec, q_val); // Get the IE, initial weighted mass kinetic_energy = - ComputeVolumeIntegral(H1, dim, NE, NQ, Q1D, dim, 2.0, qdata.rho0DetJ0w, q_val); + ComputeVolumeIntegral(H1, dim, NE, NQ, Q1D, dim, 2.0_r, + qdata.rho0DetJ0w, q_val); } - MPI_Allreduce(&kinetic_energy, &glob_ke, 1, MPI_DOUBLE, MPI_SUM, + MPI_Allreduce(&kinetic_energy, &glob_ke, 1, MPITypeMap::mpi_type, + MPI_SUM, H1.GetParMesh()->GetComm()); - return 0.5*glob_ke; + return 0.5_r*glob_ke; } void LagrangianHydroOperator::PrintTimingData(bool IamRoot, int steps, @@ -796,12 +803,12 @@ void LagrangianHydroOperator::PrintTimingData(bool IamRoot, int steps, } // Smooth transition between 0 and 1 for x in [-eps, eps]. -MFEM_HOST_DEVICE inline double smooth_step_01(double x, double eps) +MFEM_HOST_DEVICE inline real_t smooth_step_01(real_t x, real_t eps) { - const double y = (x + eps) / (2.0 * eps); - if (y < 0.0) { return 0.0; } - if (y > 1.0) { return 1.0; } - return (3.0 - 2.0 * y) * y * y; + const real_t y = (x + eps) / (2.0_r * eps); + if (y < 0.0_r) { return 0.0_r; } + if (y > 1.0_r) { return 1.0_r; } + return (3.0_r - 2.0_r * y) * y * y; } void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const @@ -833,11 +840,11 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const int nzones_batch = 3; const int nbatches = NE / nzones_batch + 1; // +1 for the remainder. int nqp_batch = nqp * nzones_batch; - double *gamma_b = new double[nqp_batch], - *rho_b = new double[nqp_batch], - *e_b = new double[nqp_batch], - *p_b = new double[nqp_batch], - *cs_b = new double[nqp_batch]; + real_t *gamma_b = new real_t[nqp_batch], + *rho_b = new real_t[nqp_batch], + *e_b = new real_t[nqp_batch], + *p_b = new real_t[nqp_batch], + *cs_b = new real_t[nqp_batch]; // Jacobians of reference->physical transformations for all quadrature points // in the batch. DenseTensor *Jpr_b = new DenseTensor[nzones_batch]; @@ -852,7 +859,7 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const nqp_batch = nqp * nzones_batch; } - double min_detJ = std::numeric_limits::infinity(); + real_t min_detJ = std::numeric_limits::infinity(); for (int z = 0; z < nzones_batch; z++) { ElementTransformation *T = H1.GetElementTransformation(z_id); @@ -863,13 +870,13 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const const IntegrationPoint &ip = ir.IntPoint(q); T->SetIntPoint(&ip); Jpr_b[z](q) = T->Jacobian(); - const double detJ = Jpr_b[z](q).Det(); + const real_t detJ = Jpr_b[z](q).Det(); min_detJ = fmin(min_detJ, detJ); const int idx = z * nqp + q; // Assuming piecewise constant gamma that moves with the mesh. gamma_b[idx] = gamma_gf(z_id); rho_b[idx] = qdata.rho0DetJ0w(z_id*nqp + q) / detJ / ip.weight; - e_b[idx] = fmax(0.0, e_vals(q)); + e_b[idx] = fmax(0.0_r, e_vals(q)); } ++z_id; } @@ -889,11 +896,11 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const // not to store the Jacobians for all batched quadrature points. const DenseMatrix &Jpr = Jpr_b[z](q); CalcInverse(Jpr, Jinv); - const double detJ = Jpr.Det(), rho = rho_b[z*nqp + q], + const real_t detJ = Jpr.Det(), rho = rho_b[z*nqp + q], p = p_b[z*nqp + q], sound_speed = cs_b[z*nqp + q]; - stress = 0.0; + stress = 0.0_r; for (int d = 0; d < dim; d++) { stress(d, d) = -p; } - double visc_coeff = 0.0; + real_t visc_coeff = 0.0_r; if (use_viscosity) { // Compression-based length scale at the point. The first @@ -902,20 +909,21 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const // relative change of the initial length scale. v.GetVectorGradient(*T, sgrad_v); - double vorticity_coeff = 1.0; + real_t vorticity_coeff = 1.0_r; if (use_vorticity) { - const double grad_norm = sgrad_v.FNorm(); - const double div_v = fabs(sgrad_v.Trace()); - vorticity_coeff = (grad_norm > 0.0) ? div_v / grad_norm : 1.0; + const real_t grad_norm = sgrad_v.FNorm(); + const real_t div_v = fabs(sgrad_v.Trace()); + vorticity_coeff = + (grad_norm > 0.0_r) ? div_v / grad_norm : 1.0_r; } sgrad_v.Symmetrize(); - double eig_val_data[3], eig_vec_data[9]; + real_t eig_val_data[3], eig_vec_data[9]; if (dim==1) { eig_val_data[0] = sgrad_v(0, 0); - eig_vec_data[0] = 1.; + eig_vec_data[0] = 1.0_r; } else { sgrad_v.CalcEigenvalues(eig_val_data, eig_vec_data); } Vector compr_dir(eig_vec_data, dim); @@ -923,38 +931,38 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const mfem::Mult(Jpr, qdata.Jac0inv(z_id*nqp + q), Jpi); Vector ph_dir(dim); Jpi.Mult(compr_dir, ph_dir); // Change of the initial mesh size in the compression direction. - const double h = qdata.h0 * ph_dir.Norml2() / + const real_t h = qdata.h0 * ph_dir.Norml2() / compr_dir.Norml2(); // Measure of maximal compression. - const double mu = eig_val_data[0]; - visc_coeff = 2.0 * rho * h * h * fabs(mu); + const real_t mu = eig_val_data[0]; + visc_coeff = 2.0_r * rho * h * h * fabs(mu); // The following represents a "smooth" version of the statement // "if (mu < 0) visc_coeff += 0.5 rho h sound_speed". Note that // eps must be scaled appropriately if a different unit system is // being used. - const double eps = 1e-12; - visc_coeff += 0.5 * rho * h * sound_speed * vorticity_coeff * - (1.0 - smooth_step_01(mu - 2.0 * eps, eps)); + const real_t eps = 1e-12_r; + visc_coeff += 0.5_r * rho * h * sound_speed * vorticity_coeff * + (1.0_r - smooth_step_01(mu - 2.0_r * eps, eps)); stress.Add(visc_coeff, sgrad_v); } // Time step estimate at the point. Here the more relevant length // scale is related to the actual mesh deformation; we use the min // singular value of the ref->physical Jacobian. In addition, the // time step estimate should be aware of the presence of shocks. - const double h_min = - Jpr.CalcSingularvalue(dim-1) / (double) H1.GetOrder(0); - const double inv_dt = sound_speed / h_min + - 2.5 * visc_coeff / rho / h_min / h_min; - if (min_detJ < 0.0) + const real_t h_min = Jpr.CalcSingularvalue(dim-1) / + static_cast(H1.GetOrder(0)); + const real_t inv_dt = sound_speed / h_min + + 2.5_r * visc_coeff / rho / h_min / h_min; + if (min_detJ < 0.0_r) { // This will force repetition of the step with smaller dt. - qdata.dt_est = 0.0; + qdata.dt_est = 0.0_r; } else { - if (inv_dt>0.0) + if (inv_dt > 0.0_r) { - qdata.dt_est = fmin(qdata.dt_est, cfl*(1.0/inv_dt)); + qdata.dt_est = fmin(qdata.dt_est, cfl*(1.0_r/inv_dt)); } } // Quadrature data for partial assembly of the force operator. @@ -987,21 +995,21 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const /// Trace of a square matrix template MFEM_HOST_DEVICE inline -double Trace(const T * __restrict__ data) +real_t Trace(const T * __restrict__ data) { - double t = 0.0; + real_t t = 0.0_r; for (int i = 0; i < W; i++) { t += data[i+i*H]; } return t; } template MFEM_HOST_DEVICE static inline -void SFNorm(double &scale_factor, double &scaled_fnorm2, +void SFNorm(real_t &scale_factor, real_t &scaled_fnorm2, const T * __restrict__ data) { int i; constexpr int hw = H * W; - T max_norm = 0.0, entry, fnorm2; + T max_norm = 0.0_r, entry, fnorm2; for (i = 0; i < hw; i++) { @@ -1012,13 +1020,13 @@ void SFNorm(double &scale_factor, double &scaled_fnorm2, } } - if (max_norm == 0.0) + if (max_norm == 0.0_r) { - scale_factor = scaled_fnorm2 = 0.0; + scale_factor = scaled_fnorm2 = 0.0_r; return; } - fnorm2 = 0.0; + fnorm2 = 0.0_r; for (i = 0; i < hw; i++) { entry = data[i] / max_norm; @@ -1032,9 +1040,9 @@ void SFNorm(double &scale_factor, double &scaled_fnorm2, /// Compute the Frobenius norm of the matrix template MFEM_HOST_DEVICE inline -double FNorm(const T * __restrict__ data) +real_t FNorm(const T * __restrict__ data) { - double s, n2; + real_t s, n2; SFNorm(s, n2, data); return s*sqrt(n2); } @@ -1044,69 +1052,69 @@ void QUpdateBody(const int NE, const int e, const int NQ, const int q, const bool use_viscosity, const bool use_vorticity, - const double h0, - const double h1order, - const double cfl, - const double infinity, - double* __restrict__ Jinv, - double* __restrict__ stress, - double* __restrict__ sgrad_v, - double* __restrict__ eig_val_data, - double* __restrict__ eig_vec_data, - double* __restrict__ compr_dir, - double* __restrict__ Jpi, - double* __restrict__ ph_dir, - double* __restrict__ stressJiT, - const double* __restrict__ d_gamma, - const double* __restrict__ d_weights, - const double* __restrict__ d_Jacobians, - const double* __restrict__ d_rho0DetJ0w, - const double* __restrict__ d_e_quads, - const double* __restrict__ d_grad_v_ext, - const double* __restrict__ d_Jac0inv, - double *d_dt_est, - double *d_stressJinvT) + const real_t h0, + const real_t h1order, + const real_t cfl, + const real_t infinity, + real_t* __restrict__ Jinv, + real_t* __restrict__ stress, + real_t* __restrict__ sgrad_v, + real_t* __restrict__ eig_val_data, + real_t* __restrict__ eig_vec_data, + real_t* __restrict__ compr_dir, + real_t* __restrict__ Jpi, + real_t* __restrict__ ph_dir, + real_t* __restrict__ stressJiT, + const real_t* __restrict__ d_gamma, + const real_t* __restrict__ d_weights, + const real_t* __restrict__ d_Jacobians, + const real_t* __restrict__ d_rho0DetJ0w, + const real_t* __restrict__ d_e_quads, + const real_t* __restrict__ d_grad_v_ext, + const real_t* __restrict__ d_Jac0inv, + real_t *d_dt_est, + real_t *d_stressJinvT) { constexpr int DIM2 = DIM*DIM; - double min_detJ = infinity; + real_t min_detJ = infinity; const int eq = e * NQ + q; - const double gamma = d_gamma[e]; - const double weight = d_weights[q]; - const double inv_weight = 1. / weight; - const double *J = d_Jacobians + DIM2*(NQ*e + q); - const double detJ = kernels::Det(J); + const real_t gamma = d_gamma[e]; + const real_t weight = d_weights[q]; + const real_t inv_weight = 1.0_r / weight; + const real_t *J = d_Jacobians + DIM2*(NQ*e + q); + const real_t detJ = kernels::Det(J); min_detJ = fmin(min_detJ, detJ); kernels::CalcInverse(J, Jinv); - const double R = inv_weight * d_rho0DetJ0w[eq] / detJ; - const double E = fmax(0.0, d_e_quads[eq]); - const double P = (gamma - 1.0) * R * E; - const double S = sqrt(gamma * (gamma - 1.0) * E); - for (int k = 0; k < DIM2; k++) { stress[k] = 0.0; } + const real_t R = inv_weight * d_rho0DetJ0w[eq] / detJ; + const real_t E = fmax(0.0_r, d_e_quads[eq]); + const real_t P = (gamma - 1.0_r) * R * E; + const real_t S = sqrt(gamma * (gamma - 1.0_r) * E); + for (int k = 0; k < DIM2; k++) { stress[k] = 0.0_r; } for (int d = 0; d < DIM; d++) { stress[d*DIM+d] = -P; } - double visc_coeff = 0.0; + real_t visc_coeff = 0.0_r; if (use_viscosity) { // Compression-based length scale at the point. The first // eigenvector of the symmetric velocity gradient gives the // direction of maximal compression. This is used to define the // relative change of the initial length scale. - const double *dV = d_grad_v_ext + DIM2*(NQ*e + q); + const real_t *dV = d_grad_v_ext + DIM2*(NQ*e + q); kernels::Mult(DIM, DIM, DIM, dV, Jinv, sgrad_v); - double vorticity_coeff = 1.0; + real_t vorticity_coeff = 1.0_r; if (use_vorticity) { - const double grad_norm = FNorm(sgrad_v); - const double div_v = fabs(Trace(sgrad_v)); - vorticity_coeff = (grad_norm > 0.0) ? div_v / grad_norm : 1.0; + const real_t grad_norm = FNorm(sgrad_v); + const real_t div_v = fabs(Trace(sgrad_v)); + vorticity_coeff = (grad_norm > 0.0_r) ? div_v / grad_norm : 1.0_r; } kernels::Symmetrize(DIM, sgrad_v); if (DIM == 1) { eig_val_data[0] = sgrad_v[0]; - eig_vec_data[0] = 1.; + eig_vec_data[0] = 1.0_r; } else { @@ -1117,40 +1125,40 @@ void QUpdateBody(const int NE, const int e, kernels::Mult(DIM, DIM, DIM, J, d_Jac0inv + eq*DIM*DIM, Jpi); kernels::Mult(DIM, DIM, Jpi, compr_dir, ph_dir); // Change of the initial mesh size in the compression direction. - const double ph_dir_nl2 = kernels::Norml2(DIM, ph_dir); - const double compr_dir_nl2 = kernels::Norml2(DIM, compr_dir); - const double H = h0 * ph_dir_nl2 / compr_dir_nl2; + const real_t ph_dir_nl2 = kernels::Norml2(DIM, ph_dir); + const real_t compr_dir_nl2 = kernels::Norml2(DIM, compr_dir); + const real_t H = h0 * ph_dir_nl2 / compr_dir_nl2; // Measure of maximal compression. - const double mu = eig_val_data[0]; - visc_coeff = 2.0 * R * H * H * fabs(mu); + const real_t mu = eig_val_data[0]; + visc_coeff = 2.0_r * R * H * H * fabs(mu); // The following represents a "smooth" version of the statement // "if (mu < 0) visc_coeff += 0.5 rho h sound_speed". Note that // eps must be scaled appropriately if a different unit system is // being used. - const double eps = 1e-12; - visc_coeff += 0.5 * R * H * S * vorticity_coeff * - (1.0 - smooth_step_01(mu-2.0*eps, eps)); + const real_t eps = 1e-12_r; + visc_coeff += 0.5_r * R * H * S * vorticity_coeff * + (1.0_r - smooth_step_01(mu-2.0_r*eps, eps)); kernels::Add(DIM, DIM, visc_coeff, stress, sgrad_v, stress); } // Time step estimate at the point. Here the more relevant length // scale is related to the actual mesh deformation; we use the min // singular value of the ref->physical Jacobian. In addition, the // time step estimate should be aware of the presence of shocks. - const double sv = kernels::CalcSingularvalue(J, DIM - 1); - const double h_min = sv / h1order; - const double ih_min = 1. / h_min; - const double irho_ih_min_sq = ih_min * ih_min / R ; - const double idt = S * ih_min + 2.5 * visc_coeff * irho_ih_min_sq; - if (min_detJ < 0.0) + const real_t sv = kernels::CalcSingularvalue(J, DIM - 1); + const real_t h_min = sv / h1order; + const real_t ih_min = 1.0_r / h_min; + const real_t irho_ih_min_sq = ih_min * ih_min / R ; + const real_t idt = S * ih_min + 2.5_r * visc_coeff * irho_ih_min_sq; + if (min_detJ < 0.0_r) { // This will force repetition of the step with smaller dt. - d_dt_est[eq] = 0.0; + d_dt_est[eq] = 0.0_r; } else { - if (idt > 0.0) + if (idt > 0.0_r) { - const double cfl_inv_dt = cfl / idt; + const real_t cfl_inv_dt = cfl / idt; d_dt_est[eq] = fmin(d_dt_est[eq], cfl_inv_dt); } } @@ -1173,7 +1181,7 @@ static void Rho0DetJ0Vol(const int dim, const int NE, ParFiniteElementSpace &L2, const ParGridFunction &rho0, QuadratureData &qdata, - double &volume) + real_t &volume) { const int NQ = ir.GetNPoints(); const int Q1D = IntRules.Get(Geometry::SEGMENT,ir.GetOrder()).GetNPoints(); @@ -1189,7 +1197,7 @@ static void Rho0DetJ0Vol(const int dim, const int NE, const auto J = Reshape(geom->J.Read(), NQ, dim, dim, NE); const auto detJ = Reshape(geom->detJ.Read(), NQ, NE); auto V = Reshape(qdata.rho0DetJ0w.Write(), NQ, NE); - Memory &Jinv_m = qdata.Jac0inv.GetMemory(); + Memory &Jinv_m = qdata.Jac0inv.GetMemory(); const MemoryClass mc = Device::GetMemoryClass(); const int Ji_total_size = qdata.Jac0inv.TotalSize(); auto invJ = Reshape(Jinv_m.Write(mc, Ji_total_size), dim, dim, NQ, NE); @@ -1206,19 +1214,19 @@ static void Rho0DetJ0Vol(const int dim, const int NE, MFEM_FOREACH_THREAD(qx,x,Q1D) { const int q = qx + qy * Q1D; - const double J11 = J(q,0,0,e); - const double J12 = J(q,1,0,e); - const double J21 = J(q,0,1,e); - const double J22 = J(q,1,1,e); - const double det = detJ(q,e); + const real_t J11 = J(q,0,0,e); + const real_t J12 = J(q,1,0,e); + const real_t J21 = J(q,0,1,e); + const real_t J22 = J(q,1,1,e); + const real_t det = detJ(q,e); V(q,e) = W[q] * R(q,e) * det; - const double r_idetJ = 1.0 / det; + const real_t r_idetJ = 1.0_r / det; invJ(0,0,q,e) = J22 * r_idetJ; invJ(1,0,q,e) = -J12 * r_idetJ; invJ(0,1,q,e) = -J21 * r_idetJ; invJ(1,1,q,e) = J11 * r_idetJ; A(q,e) = W[q] * det; - O(q,e) = 1.0; + O(q,e) = 1.0_r; } } }); @@ -1234,12 +1242,15 @@ static void Rho0DetJ0Vol(const int dim, const int NE, MFEM_FOREACH_THREAD(qx,x,Q1D) { const int q = qx + (qy + qz * Q1D) * Q1D; - const double J11 = J(q,0,0,e), J12 = J(q,0,1,e), J13 = J(q,0,2,e); - const double J21 = J(q,1,0,e), J22 = J(q,1,1,e), J23 = J(q,1,2,e); - const double J31 = J(q,2,0,e), J32 = J(q,2,1,e), J33 = J(q,2,2,e); - const double det = detJ(q,e); + const real_t J11 = J(q,0,0,e), J12 = J(q,0,1,e), + J13 = J(q,0,2,e); + const real_t J21 = J(q,1,0,e), J22 = J(q,1,1,e), + J23 = J(q,1,2,e); + const real_t J31 = J(q,2,0,e), J32 = J(q,2,1,e), + J33 = J(q,2,2,e); + const real_t det = detJ(q,e); V(q,e) = W[q] * R(q,e) * det; - const double r_idetJ = 1.0 / det; + const real_t r_idetJ = 1.0_r / det; invJ(0,0,q,e) = r_idetJ * ((J22 * J33)-(J23 * J32)); invJ(1,0,q,e) = r_idetJ * ((J32 * J13)-(J33 * J12)); invJ(2,0,q,e) = r_idetJ * ((J12 * J23)-(J13 * J22)); @@ -1250,7 +1261,7 @@ static void Rho0DetJ0Vol(const int dim, const int NE, invJ(1,2,q,e) = r_idetJ * ((J31 * J12)-(J32 * J11)); invJ(2,2,q,e) = r_idetJ * ((J11 * J22)-(J12 * J21)); A(q,e) = W[q] * det; - O(q,e) = 1.0; + O(q,e) = 1.0_r; } } } @@ -1264,12 +1275,12 @@ template static inline void QKernel(const int NE, const int NQ, const bool use_viscosity, const bool use_vorticity, - const double h0, - const double h1order, - const double cfl, - const double infinity, + const real_t h0, + const real_t h1order, + const real_t cfl, + const real_t infinity, const ParGridFunction &gamma_gf, - const Array &weights, + const Array &weights, const Vector &Jacobians, const Vector &rho0DetJ0w, const Vector &e_quads, @@ -1292,15 +1303,15 @@ void QKernel(const int NE, const int NQ, { MFEM_FORALL_2D(e, NE, Q1D, Q1D, 1, { - double Jinv[DIM2]; - double stress[DIM2]; - double sgrad_v[DIM2]; - double eig_val_data[3]; - double eig_vec_data[9]; - double compr_dir[DIM]; - double Jpi[DIM2]; - double ph_dir[DIM]; - double stressJiT[DIM2]; + real_t Jinv[DIM2]; + real_t stress[DIM2]; + real_t sgrad_v[DIM2]; + real_t eig_val_data[3]; + real_t eig_vec_data[9]; + real_t compr_dir[DIM]; + real_t Jpi[DIM2]; + real_t ph_dir[DIM]; + real_t stressJiT[DIM2]; MFEM_FOREACH_THREAD(qx,x,Q1D) { MFEM_FOREACH_THREAD(qy,y,Q1D) @@ -1321,15 +1332,15 @@ void QKernel(const int NE, const int NQ, { MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D, { - double Jinv[DIM2]; - double stress[DIM2]; - double sgrad_v[DIM2]; - double eig_val_data[3]; - double eig_vec_data[9]; - double compr_dir[DIM]; - double Jpi[DIM2]; - double ph_dir[DIM]; - double stressJiT[DIM2]; + real_t Jinv[DIM2]; + real_t stress[DIM2]; + real_t sgrad_v[DIM2]; + real_t eig_val_data[3]; + real_t eig_vec_data[9]; + real_t compr_dir[DIM]; + real_t Jpi[DIM2]; + real_t ph_dir[DIM]; + real_t stressJiT[DIM2]; MFEM_FOREACH_THREAD(qx,x,Q1D) { MFEM_FOREACH_THREAD(qy,y,Q1D) @@ -1358,8 +1369,8 @@ void QUpdate::UpdateQuadratureData(const Vector &S, QuadratureData &qdata) LAGHOS_CALI_MARK_BEGIN("QUpdate-UpdateQuadratureData"); Vector* S_p = const_cast(&S); const int H1_size = H1.GetVSize(); - const double h1order = (double) H1.GetOrder(0); - const double infinity = std::numeric_limits::infinity(); + const real_t h1order = static_cast(H1.GetOrder(0)); + const real_t infinity = std::numeric_limits::infinity(); ParGridFunction x, v, e; x.MakeRef(&H1,*S_p, 0); H1R->Mult(x, e_vec); @@ -1376,10 +1387,10 @@ void QUpdate::UpdateQuadratureData(const Vector &S, QuadratureData &qdata) typedef void (*fQKernel)(const int NE, const int NQ, const bool use_viscosity, const bool use_vorticity, - const double h0, const double h1order, - const double cfl, const double infinity, + const real_t h0, const real_t h1order, + const real_t cfl, const real_t infinity, const ParGridFunction &gamma_gf, - const Array &weights, + const Array &weights, const Vector &Jacobians, const Vector &rho0DetJ0w, const Vector &e_quads, const Vector &grad_v_ext, const DenseTensor &Jac0inv, @@ -1413,7 +1424,7 @@ void QUpdate::UpdateQuadratureData(const Vector &S, QuadratureData &qdata) void LagrangianHydroOperator::AssembleForceMatrix() const { if (forcemat_is_assembled || p_assembly) { return; } - Force = 0.0; + Force = 0.0_r; LAGHOS_DEVICE_SYNC; timer.sw_force.Start(); LAGHOS_CALI_MARK_BEGIN("LagrangianHydroOperator-AssembleForceMatrix"); @@ -1440,11 +1451,11 @@ void RK2AvgSolver::Init(TimeDependentOperator &tdop) V.SetSize(block_offsets[1], mem_type); V.UseDevice(true); dS_dt.Update(block_offsets, mem_type); - dS_dt = 0.0; + dS_dt = 0.0_r; S0.Update(block_offsets, mem_type); } -void RK2AvgSolver::Step(Vector &S, double &t, double &dt) +void RK2AvgSolver::Step(Vector &S, real_t &t, real_t &dt) { // The monolithic BlockVector stores the unknown fields as follows: // (Position, Velocity, Specific Internal Energy). @@ -1464,18 +1475,18 @@ void RK2AvgSolver::Step(Vector &S, double &t, double &dt) hydro_oper->UpdateMesh(S); hydro_oper->SolveVelocity(S, dS_dt); // V = v0 + 0.5 * dt * dv_dt; - add(v0, 0.5 * dt, dv_dt, V); + add(v0, 0.5_r * dt, dv_dt, V); hydro_oper->SolveEnergy(S, V, dS_dt); dx_dt = V; // -- 2. // S = S0 + 0.5 * dt * dS_dt; - add(S0, 0.5 * dt, dS_dt, S); + add(S0, 0.5_r * dt, dS_dt, S); hydro_oper->ResetQuadratureData(); hydro_oper->UpdateMesh(S); hydro_oper->SolveVelocity(S, dS_dt); // V = v0 + 0.5 * dt * dv_dt; - add(v0, 0.5 * dt, dv_dt, V); + add(v0, 0.5_r * dt, dv_dt, V); hydro_oper->SolveEnergy(S, V, dS_dt); dx_dt = V; diff --git a/laghos_solver.hpp b/laghos_solver.hpp index f1c5e42f..8b7b2d57 100644 --- a/laghos_solver.hpp +++ b/laghos_solver.hpp @@ -60,7 +60,7 @@ class QUpdate private: const int dim, vdim, NQ, NE, Q1D; const bool use_viscosity, use_vorticity; - const double cfl; + const real_t cfl; TimingData *timer; const IntegrationRule &ir; ParFiniteElementSpace &H1, &L2; @@ -71,7 +71,7 @@ class QUpdate public: QUpdate(const int d, const int ne, const int q1d, const bool visc, const bool vort, - const double cfl, TimingData *t, + const real_t cfl, TimingData *t, const ParGridFunction &gamma_gf, const IntegrationRule &ir, ParFiniteElementSpace &h1, ParFiniteElementSpace &l2): @@ -112,11 +112,11 @@ class LagrangianHydroOperator : public TimeDependentOperator mutable ParGridFunction x_gf; const Array &ess_tdofs; const int dim, NE, l2dofs_cnt, h1dofs_cnt, source_type; - const double cfl; + const real_t cfl; const bool use_viscosity, use_vorticity, p_assembly; - const double cg_rel_tol; + const real_t cg_rel_tol; const int cg_max_iter; - const double ftz_tol; + const real_t ftz_tol; const ParGridFunction &gamma_gf; // Velocity mass matrix and local inverses of the energy mass matrices. These // are constant in time, due to the pointwise mass conservation property. @@ -148,14 +148,14 @@ class LagrangianHydroOperator : public TimeDependentOperator mutable ParGridFunction rhs_c_gf, dvc_gf; mutable Array c_tdofs[3]; - virtual void ComputeMaterialProperties(int nvalues, const double gamma[], - const double rho[], const double e[], - double p[], double cs[]) const + virtual void ComputeMaterialProperties(int nvalues, const real_t gamma[], + const real_t rho[], const real_t e[], + real_t p[], real_t cs[]) const { for (int v = 0; v < nvalues; v++) { - p[v] = (gamma[v] - 1.0) * rho[v] * e[v]; - cs[v] = sqrt(gamma[v] * (gamma[v]-1.0) * e[v]); + p[v] = (gamma[v] - 1.0_r) * rho[v] * e[v]; + cs[v] = sqrt(gamma[v] * (gamma[v]-1.0_r) * e[v]); } } @@ -171,9 +171,10 @@ class LagrangianHydroOperator : public TimeDependentOperator ParGridFunction &rho0_gf, ParGridFunction &gamma_gf, const int source, - const double cfl, + const real_t cfl, const bool visc, const bool vort, const bool pa, - const double cgt, const int cgiter, double ftz_tol, + const real_t cgt, const int cgiter, + real_t ftz_tol, const int order_q); ~LagrangianHydroOperator(); @@ -188,15 +189,15 @@ class LagrangianHydroOperator : public TimeDependentOperator void UpdateMesh(const Vector &S) const; // Calls UpdateQuadratureData to compute the new qdata.dt_estimate. - double GetTimeStepEstimate(const Vector &S) const; + real_t GetTimeStepEstimate(const Vector &S) const; void ResetTimeStepEstimate() const; void ResetQuadratureData() const { qdata_is_current = false; } // The density values, which are stored only at some quadrature points, // are projected as a ParGridFunction. void ComputeDensity(ParGridFunction &rho) const; - double InternalEnergy(const ParGridFunction &e) const; - double KineticEnergy(const ParGridFunction &v) const; + real_t InternalEnergy(const ParGridFunction &e) const; + real_t KineticEnergy(const ParGridFunction &v) const; int GetH1VSize() const { return H1.GetVSize(); } const Array &GetBlockOffsets() const { return block_offsets; } @@ -208,12 +209,13 @@ class LagrangianHydroOperator : public TimeDependentOperator class TaylorCoefficient : public Coefficient { public: - virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip) + virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) { Vector x(2); T.Transform(ip, x); - return 3.0 / 8.0 * M_PI * ( cos(3.0*M_PI*x(0)) * cos(M_PI*x(1)) - - cos(M_PI*x(0)) * cos(3.0*M_PI*x(1)) ); + return 3.0_r / 8.0_r * M_PI * + ( cos(3.0_r*M_PI*x(0)) * cos(M_PI*x(1)) - + cos(M_PI*x(0)) * cos(3.0_r*M_PI*x(1)) ); } }; @@ -226,7 +228,7 @@ class RTCoefficient : public VectorCoefficient virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip) { - V = 0.0; V(1) = -1.0; + V = 0.0_r; V(1) = -1.0_r; } }; @@ -239,7 +241,7 @@ class HydroODESolver : public ODESolver public: HydroODESolver() : hydro_oper(NULL) { } virtual void Init(TimeDependentOperator&); - virtual void Step(Vector&, double&, double&) + virtual void Step(Vector&, real_t&, real_t&) { MFEM_ABORT("Time stepping is undefined."); } }; @@ -251,7 +253,7 @@ class RK2AvgSolver : public HydroODESolver public: RK2AvgSolver() { } virtual void Init(TimeDependentOperator &_f); - virtual void Step(Vector &S, double &t, double &dt); + virtual void Step(Vector &S, real_t &t, real_t &dt); }; } // namespace mfem diff --git a/sedov/sedov.cpp b/sedov/sedov.cpp index 7fe7e5f1..80233a55 100644 --- a/sedov/sedov.cpp +++ b/sedov/sedov.cpp @@ -43,7 +43,7 @@ int main(int argc, char *argv[]) { int ny = 2; int nz = 2; int order_q = 4; - double t_final = 0.6; + real_t t_final = 0.6; const char *basename = "results/Sedov"; real_t Sx = 1, Sy = 1, Sz = 1;