diff --git a/README.md b/README.md index bf994bc9..094758ae 100644 --- a/README.md +++ b/README.md @@ -165,7 +165,7 @@ Clone and build the parallel version of MFEM: ``` The above uses the `laghos-v1.0` tag of MFEM, which is guaranteed to work with Laghos v1.0. Alternatively, one can use the latest versions of the MFEM and -Laghos `master` branches (provided there are no conflicts. See the [MFEM +Laghos `master` branches (provided there are no conflicts). See the [MFEM building page](http://mfem.org/building/) for additional details. (Optional) Clone and build GLVis: @@ -255,31 +255,33 @@ partial assembly option (`-pa`). Some sample runs in 2D and 3D respectively are: ```sh -mpirun -np 8 laghos -p 1 -m data/square01_quad.mesh -rs 3 -tf 0.8 -no-vis -pa -mpirun -np 8 laghos -p 1 -m data/cube01_hex.mesh -rs 2 -tf 0.6 -no-vis -pa +mpirun -np 8 laghos -p 1 -m data/square01_quad.mesh -rs 3 -tf 0.8 -pa +mpirun -np 8 laghos -p 1 -m data/cube01_hex.mesh -rs 2 -tf 0.6 -vis -pa ``` -The latter produces the following density plot (when run with the `-vis` instead -of the `-no-vis` option) +The latter produces the following density plot (notice the `-vis` option) ![Sedov blast image](data/sedov.png) -#### Taylor-Green vortex +#### Taylor-Green and Gresho vortices -Laghos includes also a smooth test problem, that exposes all the principal +Laghos includes also smooth test problems that expose all the principal computational kernels of the problem except for the artificial viscosity evaluation. Some sample runs in 2D and 3D respectively are: ```sh -mpirun -np 8 laghos -p 0 -m data/square01_quad.mesh -rs 3 -tf 0.5 -no-vis -pa -mpirun -np 8 laghos -p 0 -m data/cube01_hex.mesh -rs 1 -cfl 0.1 -tf 0.25 -no-vis -pa +mpirun -np 8 laghos -p 0 -m data/square01_quad.mesh -rs 3 -tf 0.5 -pa +mpirun -np 8 laghos -p 0 -m data/cube01_hex.mesh -rs 1 -cfl 0.1 -tf 0.25 -vis -pa +mpirun -np 8 laghos -p 4 -m data/square_gresho.mesh -rs 3 -ok 3 -ot 2 -tf 0.62 -s 7 -vis -pa ``` -The latter produces the following velocity magnitude plot (when run with the -`-vis` instead of the `-no-vis` option) +The latter produce the following velocity magnitude plots (notice the `-vis` option) -![Taylor-Green image](data/tg.png) + +
+ +
#### Triple-point problem @@ -288,27 +290,26 @@ thus examining the complex computational abilities of Laghos. Some sample runs in 2D and 3D respectively are: ```sh -mpirun -np 8 laghos -p 3 -m data/rectangle01_quad.mesh -rs 2 -tf 2.5 -cfl 0.025 -no-vis -pa -mpirun -np 8 laghos -p 3 -m data/box01_hex.mesh -rs 1 -tf 2.5 -cfl 0.05 -no-vis -pa +mpirun -np 8 laghos -p 3 -m data/rectangle01_quad.mesh -rs 2 -tf 2.5 -cfl 0.025 -pa +mpirun -np 8 laghos -p 3 -m data/box01_hex.mesh -rs 1 -tf 2.5 -cfl 0.05 -vis -pa ``` -The latter produces the following specific internal energy plot (when run with -the `-vis` instead of the `-no-vis` option) +The latter produces the following specific internal energy plot (notice the `-vis` option) ![Triple-point image](data/tp.png) ## Verification of Results To make sure the results are correct, we tabulate reference final iterations -(`step`), time steps (`dt`) and energies (`|e|`) for the nine runs listed above: +(`step`), time steps (`dt`) and energies (`|e|`) for the runs listed below: -1. `mpirun -np 8 laghos -p 0 -m data/square01_quad.mesh -rs 3 -tf 0.75 -no-vis -pa` -2. `mpirun -np 8 laghos -p 0 -m data/cube01_hex.mesh -rs 1 -tf 0.75 -no-vis -pa` -3. `mpirun -np 8 laghos -p 1 -m data/square01_quad.mesh -rs 3 -tf 0.8 -no-vis -pa` -4. `mpirun -np 8 laghos -p 1 -m data/cube01_hex.mesh -rs 2 -tf 0.6 -no-vis -pa` -5. `mpirun -np 8 laghos -p 2 -m data/segment01.mesh -rs 5 -tf 0.2 -no-vis -fa` -6. `mpirun -np 8 laghos -p 3 -m data/rectangle01_quad.mesh -rs 2 -tf 3.0 -no-vis -pa` -7. `mpirun -np 8 laghos -p 3 -m data/box01_hex.mesh -rs 1 -tf 3.0 -no-vis -pa` +1. `mpirun -np 8 laghos -p 0 -m data/square01_quad.mesh -rs 3 -tf 0.75 -pa` +2. `mpirun -np 8 laghos -p 0 -m data/cube01_hex.mesh -rs 1 -tf 0.75 -pa` +3. `mpirun -np 8 laghos -p 1 -m data/square01_quad.mesh -rs 3 -tf 0.8 -pa` +4. `mpirun -np 8 laghos -p 1 -m data/cube01_hex.mesh -rs 2 -tf 0.6 -pa` +5. `mpirun -np 8 laghos -p 2 -m data/segment01.mesh -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 4 -m data/square_gresho.mesh -rs 3 -ok 3 -ot 2 -tf 0.62831853 -s 7 -pa` | `run` | `step` | `dt` | `e` | | ----- | ------ | ---- | --- | @@ -318,7 +319,7 @@ To make sure the results are correct, we tabulate reference final iterations | 4. | 561 | 0.000360 | 134.0937837919 | | 5. | 414 | 0.000339 | 32.0120759615 | | 6. | 5310 | 0.000264 | 141.8348694390 | -| 7. | 937 | 0.002285 | 144.0012514765 | +| 7. | 776 | 0.000045 | 409.8243172608 | An implementation is considered valid if the final energy values are all within @@ -348,10 +349,8 @@ A sample run on the [Vulcan](https://computation.llnl.gov/computers/vulcan) BG/Q machine at LLNL is: ``` -srun -n 393216 laghos -pa -p 1 -tf 0.6 -no-vis - -pt 322 -m data/cube_12_hex.mesh - --cg-tol 0 --cg-max-iter 50 --max-steps 2 - -ok 3 -ot 2 -rs 5 -rp 3 +srun -n 393216 laghos -pa -p 1 -tf 0.6 -pt 322 -m data/cube_12_hex.mesh \ + --cg-tol 0 --cg-max-iter 50 --max-steps 2 -ok 3 -ot 2 -rs 5 -rp 3 ``` This is Q3-Q2 3D computation on 393,216 MPI ranks (24,576 nodes) that produces rates of approximately 168497, 74221, and 16696 megadofs, and a total FOM of @@ -373,6 +372,10 @@ the following versions of Laghos have been developed [OCCA](http://libocca.org/). - A [RAJA](https://software.llnl.gov/RAJA/)-based version in the [raja-dev](https://github.com/CEED/Laghos/tree/raja-dev) branch. +- An MFEM/engines-based version in the + [engines-kernels](https://github.com/CEED/Laghos/tree/engines-kernels) branches. +- Version with adaptive mesh refinement in the + [amr-dev](https://github.com/CEED/Laghos/tree/amr-dev) branch. ## Contact diff --git a/data/gresho.png b/data/gresho.png new file mode 100644 index 00000000..ba72748b Binary files /dev/null and b/data/gresho.png differ diff --git a/data/square_gresho.mesh b/data/square_gresho.mesh new file mode 100644 index 00000000..f47f2723 --- /dev/null +++ b/data/square_gresho.mesh @@ -0,0 +1,61 @@ +MFEM mesh v1.0 + +# +# MFEM Geometry Types (see mesh/geom.hpp): +# +# POINT = 0 +# SEGMENT = 1 +# TRIANGLE = 2 +# SQUARE = 3 +# TETRAHEDRON = 4 +# CUBE = 5 +# + +dimension +2 + +elements +4 +1 3 0 1 4 3 +1 3 1 2 5 4 +1 3 3 4 7 6 +1 3 4 5 8 7 + +boundary +8 +2 1 0 1 +2 1 1 2 +2 1 7 6 +2 1 8 7 +1 1 3 0 +1 1 6 3 +1 1 2 5 +1 1 5 8 + +vertices +9 + +nodes +FiniteElementSpace +FiniteElementCollection: Linear +VDim: 2 +Ordering: 0 + +-0.5 +0 +0.5 +-0.5 +0 +0.5 +-0.5 +0 +0.5 +-0.5 +-0.5 +-0.5 +0 +0 +0 +0.5 +0.5 +0.5 diff --git a/laghos.cpp b/laghos.cpp index 53ea382f..146a1c25 100644 --- a/laghos.cpp +++ b/laghos.cpp @@ -33,27 +33,18 @@ // methods for Lagrangian hydrodynamics", SIAM Journal on Scientific // Computing, (34) 2012, pp. B606–B641, https://doi.org/10.1137/120864672. // -// Sample runs: -// mpirun -np 8 laghos -p 0 -m data/square01_quad.mesh -rs 3 -tf 0.75 -// mpirun -np 8 laghos -p 0 -m data/square01_tri.mesh -rs 1 -tf 0.75 -// mpirun -np 8 laghos -p 0 -m data/cube01_hex.mesh -rs 1 -tf 2.0 -// mpirun -np 8 laghos -p 1 -m data/square01_quad.mesh -rs 3 -tf 0.8 -// mpirun -np 8 laghos -p 1 -m data/square01_quad.mesh -rs 0 -tf 0.8 -ok 7 -ot 6 -// mpirun -np 8 laghos -p 1 -m data/cube01_hex.mesh -rs 2 -tf 0.6 -// mpirun -np 8 laghos -p 2 -m data/segment01.mesh -rs 5 -tf 0.2 -// mpirun -np 8 laghos -p 3 -m data/rectangle01_quad.mesh -rs 2 -tf 3.0 -// mpirun -np 8 laghos -p 3 -m data/box01_hex.mesh -rs 1 -tf 3.0 -// // Test problems: // p = 0 --> Taylor-Green vortex (smooth problem). // p = 1 --> Sedov blast. // p = 2 --> 1D Sod shock tube. // p = 3 --> Triple point. - +// p = 4 --> Gresho vortex (smooth problem). +// +// Sample runs: see README.md, section 'Verification of Results'. +// #include "laghos_solver.hpp" -#include -#include +#include "laghos_timeinteg.hpp" #include using namespace std; @@ -281,6 +272,7 @@ int main(int argc, char *argv[]) case 3: ode_solver = new RK3SSPSolver; break; case 4: ode_solver = new RK4Solver; break; case 6: ode_solver = new RK6Solver; break; + case 7: ode_solver = new RK2AvgSolver; break; default: if (myid == 0) { @@ -326,8 +318,7 @@ int main(int argc, char *argv[]) v_gf.MakeRef(&H1FESpace, S, true_offset[1]); e_gf.MakeRef(&L2FESpace, S, true_offset[2]); - // Initialize x_gf using the starting mesh coordinates. This also links the - // mesh positions to the values in x_gf. + // Initialize x_gf using the starting mesh coordinates. pmesh->SetNodalGridFunction(&x_gf); // Initialize the velocity. @@ -371,7 +362,7 @@ int main(int argc, char *argv[]) GridFunctionCoefficient *mat_gf_coeff = new GridFunctionCoefficient(&mat_gf); // Additional details, depending on the problem. - int source = 0; bool visc; + int source = 0; bool visc = true; switch (problem) { case 0: if (pmesh->Dimension() == 2) { source = 1; } @@ -379,6 +370,7 @@ int main(int argc, char *argv[]) case 1: visc = true; break; case 2: visc = true; break; case 3: visc = true; break; + case 4: visc = false; break; default: MFEM_ABORT("Wrong problem specification!"); } @@ -393,6 +385,9 @@ int main(int argc, char *argv[]) ParGridFunction rho_gf; if (visualization || visit) { oper.ComputeDensity(rho_gf); } + const double energy_init = oper.InternalEnergy(e_gf) + + oper.KineticEnergy(v_gf); + if (visualization) { // Make sure all MPI ranks have sent their 'v' solution before initiating @@ -407,8 +402,12 @@ int main(int argc, char *argv[]) const int Ww = 350, Wh = 350; // window size int offx = Ww+10; // window offsets - VisualizeField(vis_rho, vishost, visport, rho_gf, - "Density", Wx, Wy, Ww, Wh); + if (problem != 0 && problem != 4) + { + VisualizeField(vis_rho, vishost, visport, rho_gf, + "Density", Wx, Wy, Ww, Wh); + } + Wx += offx; VisualizeField(vis_v, vishost, visport, v_gf, "Velocity", Wx, Wy, Ww, Wh); @@ -473,7 +472,9 @@ int main(int argc, char *argv[]) } else if (dt_est > 1.25 * dt) { dt *= 1.02; } - // Make sure that the mesh corresponds to the new solution state. + // Make sure that the mesh corresponds to the new solution state. This is + // needed, because some time integrators use different S-type vectors + // and the oper object might have redirected the mesh positions to those. pmesh->NewNodes(x_gf, false); if (last_step || (ti % vis_steps) == 0) @@ -502,8 +503,12 @@ int main(int argc, char *argv[]) int Ww = 350, Wh = 350; // window size int offx = Ww+10; // window offsets - VisualizeField(vis_rho, vishost, visport, rho_gf, - "Density", Wx, Wy, Ww, Wh); + if (problem != 0 && problem != 4) + { + VisualizeField(vis_rho, vishost, visport, rho_gf, + "Density", Wx, Wy, Ww, Wh); + } + Wx += offx; VisualizeField(vis_v, vishost, visport, v_gf, "Velocity", Wx, Wy, Ww, Wh); @@ -564,6 +569,30 @@ int main(int argc, char *argv[]) } oper.PrintTimingData(mpi.Root(), steps); + const double energy_final = oper.InternalEnergy(e_gf) + + oper.KineticEnergy(v_gf); + if (mpi.Root()) + { + cout << endl; + cout << "Energy diff: " << scientific << setprecision(2) + << fabs(energy_init - energy_final) << endl; + } + + // Print the error. + // 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), + error_l1 = v_gf.ComputeL1Error(v_coeff), + error_l2 = v_gf.ComputeL2Error(v_coeff); + if (mpi.Root()) + { + cout << "L_inf error: " << error_max << endl + << "L_1 error: " << error_l1 << endl + << "L_2 error: " << error_l2 << endl; + } + } + if (visualization) { vis_v.close(); @@ -590,10 +619,9 @@ double rho0(const Vector &x) { case 0: return 1.0; case 1: return 1.0; - case 2: if (x(0) < 0.5) { return 1.0; } - else { return 0.1; } - case 3: if (x(0) > 1.0 && x(1) <= 1.5) { return 1.0; } - else { return 0.125; } + case 2: return (x(0) < 0.5) ? 1.0 : 0.1; + case 3: return (x(0) > 1.0 && x(1) <= 1.5) ? 1.0 : 0.125; + case 4: return 1.0; default: MFEM_ABORT("Bad number given for problem id!"); return 0.0; } } @@ -605,12 +633,17 @@ double gamma(const Vector &x) case 0: return 5./3.; case 1: return 1.4; case 2: return 1.4; - case 3: if (x(0) > 1.0 && x(1) <= 1.5) { return 1.4; } - else { return 1.5; } + case 3: return (x(0) > 1.0 && x(1) <= 1.5) ? 1.4 : 1.5; + case 4: return 5.0 / 3.0; default: MFEM_ABORT("Bad number given for problem id!"); return 0.0; } } +double rad(double x, double y) +{ + return sqrt(x*x + y*y); +} + void v0(const Vector &x, Vector &v) { switch (problem) @@ -628,6 +661,22 @@ void v0(const Vector &x, Vector &v) case 1: v = 0.0; break; case 2: v = 0.0; break; case 3: v = 0.0; break; + case 4: + { + const double r = rad(x(0), x(1)); + if (r < 0.2) + { + v(0) = 5.0 * x(1); + v(1) = -5.0 * x(0); + } + else if (r < 0.4) + { + v(0) = 2.0 * x(1) / r - 5.0 * x(1); + v(1) = -2.0 * x(0) / r + 5.0 * x(0); + } + else { v = 0.0; } + break; + } default: MFEM_ABORT("Bad number given for problem id!"); } } @@ -652,10 +701,26 @@ double e0(const Vector &x) return val/denom; } case 1: return 0.0; // This case in initialized in main(). - case 2: if (x(0) < 0.5) { return 1.0 / rho0(x) / (gamma(x) - 1.0); } - else { return 0.1 / rho0(x) / (gamma(x) - 1.0); } - case 3: if (x(0) > 1.0) { return 0.1 / rho0(x) / (gamma(x) - 1.0); } - else { return 1.0 / rho0(x) / (gamma(x) - 1.0); } + case 2: return (x(0) < 0.5) ? 1.0 / rho0(x) / (gamma(x) - 1.0) + : 0.1 / rho0(x) / (gamma(x) - 1.0); + case 3: return (x(0) > 1.0) ? 0.1 / rho0(x) / (gamma(x) - 1.0) + : 1.0 / rho0(x) / (gamma(x) - 1.0); + 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) + { + return (5.0 + 25.0 / 2.0 * rsq) / (gamma - 1.0); + } + else if (r < 0.4) + { + 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); + } + else { return (3.0 + 4.0 * log(2.0)) / (gamma - 1.0); } + } default: MFEM_ABORT("Bad number given for problem id!"); return 0.0; } } diff --git a/laghos_assembly.hpp b/laghos_assembly.hpp index f9eee0d0..78c40805 100644 --- a/laghos_assembly.hpp +++ b/laghos_assembly.hpp @@ -19,12 +19,8 @@ #include "mfem.hpp" - #ifdef MFEM_USE_MPI -#include -#include - namespace mfem { diff --git a/laghos_solver.cpp b/laghos_solver.cpp index e0551ef5..0b638582 100644 --- a/laghos_solver.cpp +++ b/laghos_solver.cpp @@ -94,11 +94,12 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size, source_type(source_type_), cfl(cfl_), use_viscosity(visc), p_assembly(pa), cg_rel_tol(cgt), cg_max_iter(cgiter), material_pcf(material_), - Mv(&h1_fes), Me_inv(l2dofs_cnt, l2dofs_cnt, nzones), + Mv(&h1_fes), Mv_spmat_copy(), + Me(l2dofs_cnt, l2dofs_cnt, nzones), Me_inv(l2dofs_cnt, l2dofs_cnt, nzones), integ_rule(IntRules.Get(h1_fes.GetMesh()->GetElementBaseGeometry(), 3*h1_fes.GetOrder(0) + l2_fes.GetOrder(0) - 1)), quad_data(dim, nzones, integ_rule.GetNPoints()), - quad_data_is_current(false), + quad_data_is_current(false), forcemat_is_assembled(false), Force(&l2_fes, &h1_fes), ForcePA(&quad_data, h1_fes, l2_fes), VMassPA(&quad_data, H1FESpace), VMassPA_prec(H1FESpace), locEMassPA(&quad_data, l2_fes), @@ -107,13 +108,12 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size, GridFunctionCoefficient rho_coeff(&rho0); // Standard local assembly and inversion for energy mass matrices. - DenseMatrix Me(l2dofs_cnt); - DenseMatrixInverse inv(&Me); MassIntegrator mi(rho_coeff, &integ_rule); for (int i = 0; i < nzones; i++) { + DenseMatrixInverse inv(&Me(i)); mi.AssembleElementMatrix(*l2_fes.GetFE(i), - *l2_fes.GetElementTransformation(i), Me); + *l2_fes.GetElementTransformation(i), Me(i)); inv.Factor(); inv.GetInverseMatrix(Me_inv(i)); } @@ -122,6 +122,7 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size, VectorMassIntegrator *vmi = new VectorMassIntegrator(rho_coeff, &integ_rule); Mv.AddDomainIntegrator(vmi); Mv.Assemble(); + Mv_spmat_copy = Mv.SpMat(); // Values of rho0DetJ0 and Jac0inv at all quadrature points. const int nqp = integ_rule.GetNPoints(); @@ -198,47 +199,44 @@ LagrangianHydroOperator::LagrangianHydroOperator(int size, void LagrangianHydroOperator::Mult(const Vector &S, Vector &dS_dt) const { - dS_dt = 0.0; - // Make sure that the mesh positions correspond to the ones in S. This is // needed only because some mfem time integrators don't update the solution // vector at every intermediate stage (hence they don't change the mesh). - Vector* sptr = (Vector*) &S; - ParGridFunction x; - x.MakeRef(&H1FESpace, *sptr, 0); - H1FESpace.GetParMesh()->NewNodes(x, false); - - UpdateQuadratureData(S); + UpdateMesh(S); // The monolithic BlockVector stores the unknown fields as follows: - // - Position - // - Velocity - // - Specific Internal Energy - - const int VsizeL2 = L2FESpace.GetVSize(); + // (Position, Velocity, Specific Internal Energy). + Vector* sptr = (Vector*) &S; + ParGridFunction v; const int VsizeH1 = H1FESpace.GetVSize(); - - ParGridFunction v, e; v.MakeRef(&H1FESpace, *sptr, VsizeH1); - e.MakeRef(&L2FESpace, *sptr, VsizeH1*2); - - ParGridFunction dx, dv, de; - dx.MakeRef(&H1FESpace, dS_dt, 0); - dv.MakeRef(&H1FESpace, dS_dt, VsizeH1); - de.MakeRef(&L2FESpace, dS_dt, VsizeH1*2); // Set dx_dt = v (explicit). + ParGridFunction dx; + dx.MakeRef(&H1FESpace, dS_dt, 0); dx = v; - if (!p_assembly) - { - Force = 0.0; - timer.sw_force.Start(); - Force.Assemble(); - timer.sw_force.Stop(); - } + SolveVelocity(S, dS_dt); + SolveEnergy(S, v, dS_dt); + + quad_data_is_current = false; +} + +void LagrangianHydroOperator::SolveVelocity(const Vector &S, + Vector &dS_dt) const +{ + UpdateQuadratureData(S); + AssembleForceMatrix(); + + const int VsizeL2 = L2FESpace.GetVSize(); + const int VsizeH1 = H1FESpace.GetVSize(); + + // The monolithic BlockVector stores the unknown fields as follows: + // (Position, Velocity, Specific Internal Energy). + ParGridFunction dv; + dv.MakeRef(&H1FESpace, dS_dt, VsizeH1); + dv = 0.0; - // Solve for velocity. Vector one(VsizeL2), rhs(VsizeH1), B, X; one = 1.0; if (p_assembly) { @@ -285,6 +283,22 @@ void LagrangianHydroOperator::Mult(const Vector &S, Vector &dS_dt) const timer.H1cg_iter += cg.GetNumIterations(); Mv.RecoverFEMSolution(X, rhs, dv); } +} + +void LagrangianHydroOperator::SolveEnergy(const Vector &S, const Vector &v, + Vector &dS_dt) const +{ + UpdateQuadratureData(S); + AssembleForceMatrix(); + + const int VsizeL2 = L2FESpace.GetVSize(); + const int VsizeH1 = H1FESpace.GetVSize(); + + // The monolithic BlockVector stores the unknown fields as follows: + // (Position, Velocity, Specific Internal Energy). + ParGridFunction de; + de.MakeRef(&L2FESpace, dS_dt, VsizeH1*2); + de = 0.0; // Solve for energy, assemble the energy source if such exists. LinearForm *e_source = NULL; @@ -335,16 +349,18 @@ void LagrangianHydroOperator::Mult(const Vector &S, Vector &dS_dt) const } } delete e_source; +} - quad_data_is_current = false; +void LagrangianHydroOperator::UpdateMesh(const Vector &S) const +{ + Vector* sptr = (Vector*) &S; + x_gf.MakeRef(&H1FESpace, *sptr, 0); + H1FESpace.GetParMesh()->NewNodes(x_gf, false); } double LagrangianHydroOperator::GetTimeStepEstimate(const Vector &S) const { - Vector* sptr = (Vector*) &S; - ParGridFunction x; - x.MakeRef(&H1FESpace, *sptr, 0); - H1FESpace.GetParMesh()->NewNodes(x, false); + UpdateMesh(S); UpdateQuadratureData(S); double glob_dt_est; @@ -358,7 +374,7 @@ void LagrangianHydroOperator::ResetTimeStepEstimate() const quad_data.dt_est = numeric_limits::infinity(); } -void LagrangianHydroOperator::ComputeDensity(ParGridFunction &rho) +void LagrangianHydroOperator::ComputeDensity(ParGridFunction &rho) const { rho.SetSpace(&L2FESpace); @@ -382,7 +398,37 @@ void LagrangianHydroOperator::ComputeDensity(ParGridFunction &rho) } } -void LagrangianHydroOperator::PrintTimingData(bool IamRoot, int steps) +double LagrangianHydroOperator::InternalEnergy(const ParGridFunction &e) const +{ + Vector one(l2dofs_cnt), loc_e(l2dofs_cnt); + one = 1.0; + Array l2dofs; + + double loc_ie = 0.0; + for (int z = 0; z < nzones; z++) + { + L2FESpace.GetElementDofs(z, l2dofs); + e.GetSubVector(l2dofs, loc_e); + loc_ie += Me(z).InnerProduct(loc_e, one); + } + + double glob_ie; + MPI_Allreduce(&loc_ie, &glob_ie, 1, MPI_DOUBLE, MPI_SUM, + H1FESpace.GetParMesh()->GetComm()); + return glob_ie; +} + +double LagrangianHydroOperator::KineticEnergy(const ParGridFunction &v) const +{ + double loc_ke = 0.5 * Mv_spmat_copy.InnerProduct(v, v); + + double glob_ke; + MPI_Allreduce(&loc_ke, &glob_ke, 1, MPI_DOUBLE, MPI_SUM, + H1FESpace.GetParMesh()->GetComm()); + return glob_ke; +} + +void LagrangianHydroOperator::PrintTimingData(bool IamRoot, int steps) const { double my_rt[5], rt_max[5]; my_rt[0] = timer.sw_cgH1.RealTime(); @@ -577,6 +623,12 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const visc_coeff = 2.0 * rho * h * h * fabs(mu); if (mu < 0.0) { visc_coeff += 0.5 * rho * h * sound_speed; } stress.Add(visc_coeff, sgrad_v); + + // Note that the (mu < 0.0) check introduces discontinuous + // behavior. This can lead to bigger differences in results when + // there are round-offs around the zero for the min eigenvalue. + // We've observed differences between Linux and Mac for some 3D + // calculations. } // Time step estimate at the point. Here the more relevant length @@ -619,11 +671,24 @@ void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const delete [] cs_b; delete [] Jpr_b; quad_data_is_current = true; + forcemat_is_assembled = false; timer.sw_qdata.Stop(); timer.quad_tstep += nzones; } +void LagrangianHydroOperator::AssembleForceMatrix() const +{ + if (forcemat_is_assembled || p_assembly) { return; } + + Force = 0.0; + timer.sw_force.Start(); + Force.Assemble(); + timer.sw_force.Stop(); + + forcemat_is_assembled = true; +} + } // namespace hydrodynamics } // namespace mfem diff --git a/laghos_solver.hpp b/laghos_solver.hpp index a349c437..8b0f9019 100644 --- a/laghos_solver.hpp +++ b/laghos_solver.hpp @@ -22,10 +22,6 @@ #ifdef MFEM_USE_MPI -#include -#include -#include - namespace mfem { @@ -70,6 +66,9 @@ class LagrangianHydroOperator : public TimeDependentOperator ParFiniteElementSpace &H1FESpace; ParFiniteElementSpace &L2FESpace; + // Reference to the current mesh configuration. + mutable ParGridFunction x_gf; + Array &ess_tdofs; const int dim, nzones, l2dofs_cnt, h1dofs_cnt, source_type; @@ -82,7 +81,8 @@ class LagrangianHydroOperator : public TimeDependentOperator // Velocity mass matrix and local inverses of the energy mass matrices. These // are constant in time, due to the pointwise mass conservation property. mutable ParBilinearForm Mv; - DenseTensor Me_inv; + SparseMatrix Mv_spmat_copy; + DenseTensor Me, Me_inv; // Integration rule for all assemblies. const IntegrationRule &integ_rule; @@ -90,7 +90,7 @@ class LagrangianHydroOperator : public TimeDependentOperator // Data associated with each quadrature point in the mesh. These values are // recomputed at each time step. mutable QuadratureData quad_data; - mutable bool quad_data_is_current; + mutable bool quad_data_is_current, forcemat_is_assembled; // Force matrix that combines the kinematic and thermodynamic spaces. It is // assembled in each time step and then it is used to compute the final @@ -123,6 +123,7 @@ class LagrangianHydroOperator : public TimeDependentOperator } void UpdateQuadratureData(const Vector &S) const; + void AssembleForceMatrix() const; public: LagrangianHydroOperator(int size, ParFiniteElementSpace &h1_fes, @@ -135,6 +136,10 @@ class LagrangianHydroOperator : public TimeDependentOperator // Solve for dx_dt, dv_dt and de_dt. virtual void Mult(const Vector &S, Vector &dS_dt) const; + void SolveVelocity(const Vector &S, Vector &dS_dt) const; + void SolveEnergy(const Vector &S, const Vector &v, Vector &dS_dt) const; + void UpdateMesh(const Vector &S) const; + // Calls UpdateQuadratureData to compute the new quad_data.dt_estimate. double GetTimeStepEstimate(const Vector &S) const; void ResetTimeStepEstimate() const; @@ -142,9 +147,14 @@ class LagrangianHydroOperator : public TimeDependentOperator // The density values, which are stored only at some quadrature points, are // projected as a ParGridFunction. - void ComputeDensity(ParGridFunction &rho); + void ComputeDensity(ParGridFunction &rho) const; + + double InternalEnergy(const ParGridFunction &e) const; + double KineticEnergy(const ParGridFunction &v) const; + + void PrintTimingData(bool IamRoot, int steps) const; - void PrintTimingData(bool IamRoot, int steps); + int GetH1VSize() const { return H1FESpace.GetVSize(); } ~LagrangianHydroOperator(); }; diff --git a/laghos_timeinteg.cpp b/laghos_timeinteg.cpp new file mode 100644 index 00000000..04706e29 --- /dev/null +++ b/laghos_timeinteg.cpp @@ -0,0 +1,84 @@ +// 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 "laghos_timeinteg.hpp" +#include "laghos_solver.hpp" + +using namespace std; + +namespace mfem +{ + +namespace hydrodynamics +{ + +void HydroODESolver::Init(TimeDependentOperator &_f) +{ + ODESolver::Init(_f); + + hydro_oper = dynamic_cast(f); + MFEM_VERIFY(hydro_oper, "HydroSolvers expect LagrangianHydroOperator."); +} + +void RK2AvgSolver::Step(Vector &S, double &t, double &dt) +{ + const int Vsize = hydro_oper->GetH1VSize(); + Vector V(Vsize), dS_dt(S.Size()), S0(S); + + // The monolithic BlockVector stores the unknown fields as follows: + // (Position, Velocity, Specific Internal Energy). + Vector dv_dt, v0, dx_dt; + v0.SetDataAndSize(S0.GetData() + Vsize, Vsize); + dv_dt.SetDataAndSize(dS_dt.GetData() + Vsize, Vsize); + dx_dt.SetDataAndSize(dS_dt.GetData(), Vsize); + + // In each sub-step: + // - Update the global state Vector S. + // - Compute dv_dt using S. + // - Update V using dv_dt. + // - Compute de_dt and dx_dt using S and V. + + // -- 1. + // S is S0. + 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); + 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); + 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); + hydro_oper->SolveEnergy(S, V, dS_dt); + dx_dt = V; + + // -- 3. + // S = S0 + dt * dS_dt. + add(S0, dt, dS_dt, S); + hydro_oper->ResetQuadratureData(); + + t += dt; +} + +} // namespace hydrodynamics + +} // namespace mfem diff --git a/laghos_timeinteg.hpp b/laghos_timeinteg.hpp new file mode 100644 index 00000000..bb9cbc7e --- /dev/null +++ b/laghos_timeinteg.hpp @@ -0,0 +1,56 @@ +// 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 MFEM_LAGHOS_TIMEINTEG +#define MFEM_LAGHOS_TIMEINTEG + +#include "mfem.hpp" + +namespace mfem +{ + +namespace hydrodynamics +{ + +class LagrangianHydroOperator; + +class HydroODESolver : public ODESolver +{ +protected: + LagrangianHydroOperator *hydro_oper; + +public: + HydroODESolver() : hydro_oper(NULL) { } + + virtual void Init(TimeDependentOperator &_f); + + virtual void Step(Vector &S, double &t, double &dt) + { MFEM_ABORT("Time stepping is undefined."); } +}; + +class RK2AvgSolver : public HydroODESolver +{ +public: + RK2AvgSolver() { } + + virtual void Step(Vector &S, double &t, double &dt); +}; + +} // namespace hydrodynamics + +} // namespace mfem + +#endif // MFEM_LAGHOS_TIMEINTEG diff --git a/makefile b/makefile index 6595e964..fe40b03e 100644 --- a/makefile +++ b/makefile @@ -105,10 +105,10 @@ LIBS = $(strip $(LAGHOS_LIBS) $(LDFLAGS)) CCC = $(strip $(CXX) $(LAGHOS_FLAGS)) Ccc = $(strip $(CC) $(CFLAGS) $(GL_OPTS)) -SOURCE_FILES = laghos.cpp laghos_solver.cpp laghos_assembly.cpp +SOURCE_FILES = laghos.cpp laghos_solver.cpp laghos_assembly.cpp laghos_timeinteg.cpp OBJECT_FILES1 = $(SOURCE_FILES:.cpp=.o) OBJECT_FILES = $(OBJECT_FILES1:.c=.o) -HEADER_FILES = laghos_solver.hpp laghos_assembly.hpp +HEADER_FILES = laghos_solver.hpp laghos_assembly.hpp laghos_timeinteg.hpp # Targets