1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3 // reserved. See files LICENSE and NOTICE for details.
4 //
5 // This file is part of CEED, a collection of benchmarks, miniapps, software
6 // libraries and APIs for efficient high-order finite element and spectral
7 // element discretizations for exascale applications. For more information and
8 // source code availability see http://github.com/ceed.
9 //
10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11 // a collaborative effort of two U.S. Department of Energy organizations (Office
12 // of Science and the National Nuclear Security Administration) responsible for
13 // the planning and preparation of a capable exascale ecosystem, including
14 // software, applications, hardware, advanced system engineering and early
15 // testbed platforms, in support of the nation's exascale computing imperative.
16
17 #include "laghos_solver_s.hpp"
18
19 using namespace std;
20
21 namespace mfem
22 {
23
24 namespace hydrodynamics
25 {
26
VisualizeField(socketstream & sock,const char * vishost,int visport,GridFunction & gf,const char * title,int x,int y,int w,int h,bool vec)27 void VisualizeField(socketstream &sock, const char *vishost, int visport,
28 GridFunction &gf, const char *title,
29 int x, int y, int w, int h, bool vec)
30 {
31 Mesh &mesh = *gf.FESpace()->GetMesh();
32 bool newly_opened = false;
33 int connection_failed;
34
35 do
36 {
37 if (!sock.is_open() || !sock)
38 {
39 sock.open(vishost, visport);
40 sock.precision(8);
41 newly_opened = true;
42 }
43 sock << "solution\n";
44
45 mesh.Print(sock);
46 gf.Save(sock);
47
48 if (newly_opened)
49 {
50 sock << "window_title '" << title << "'\n"
51 << "window_geometry "
52 << x << " " << y << " " << w << " " << h << "\n"
53 << "keys maaAc";
54 if ( vec ) { sock << "vvv"; }
55 sock << endl;
56 }
57
58 connection_failed = !sock && !newly_opened;
59 }
60 while (connection_failed);
61 }
62
LagrangianHydroOperator(int size,FiniteElementSpace & h1_fes,FiniteElementSpace & l2_fes,Array<int> & essential_tdofs,GridFunction & rho0,int source_type_,double cfl_,Coefficient * material_,bool visc,bool pa)63 LagrangianHydroOperator::LagrangianHydroOperator(int size,
64 FiniteElementSpace &h1_fes,
65 FiniteElementSpace &l2_fes,
66 Array<int> &essential_tdofs,
67 GridFunction &rho0,
68 int source_type_, double cfl_,
69 Coefficient *material_,
70 bool visc, bool pa)
71 : TimeDependentOperator(size),
72 H1FESpace(h1_fes), L2FESpace(l2_fes),
73 H1compFESpace(h1_fes.GetMesh(), h1_fes.FEColl(), 1),
74 ess_tdofs(essential_tdofs),
75 dim(h1_fes.GetMesh()->Dimension()),
76 nzones(h1_fes.GetMesh()->GetNE()),
77 l2dofs_cnt(l2_fes.GetFE(0)->GetDof()),
78 h1dofs_cnt(h1_fes.GetFE(0)->GetDof()),
79 source_type(source_type_), cfl(cfl_),
80 use_viscosity(visc), p_assembly(pa),
81 material_pcf(material_),
82 Mv(&h1_fes), Me_inv(l2dofs_cnt, l2dofs_cnt, nzones),
83 integ_rule(IntRules.Get(h1_fes.GetMesh()->GetElementBaseGeometry(0),
84 3*h1_fes.GetOrder(0) + l2_fes.GetOrder(0) - 1)),
85 quad_data(dim, nzones, integ_rule.GetNPoints()),
86 quad_data_is_current(false),
87 Force(&l2_fes, &h1_fes), ForcePA(&quad_data, h1_fes, l2_fes),
88 VMassPA(&quad_data, H1compFESpace), locEMassPA(&quad_data, l2_fes),
89 locCG()
90 {
91 GridFunctionCoefficient rho_coeff(&rho0);
92
93 // Standard local assembly and inversion for energy mass matrices.
94 DenseMatrix Me(l2dofs_cnt);
95 DenseMatrixInverse inv(&Me);
96 MassIntegrator mi(rho_coeff, &integ_rule);
97 for (int i = 0; i < nzones; i++)
98 {
99 mi.AssembleElementMatrix(*l2_fes.GetFE(i),
100 *l2_fes.GetElementTransformation(i), Me);
101 inv.Factor();
102 inv.GetInverseMatrix(Me_inv(i));
103 }
104
105 // Standard assembly for the velocity mass matrix.
106 VectorMassIntegrator *vmi = new VectorMassIntegrator(rho_coeff, &integ_rule);
107 Mv.AddDomainIntegrator(vmi);
108 Mv.Assemble();
109
110 // Values of rho0DetJ0 and Jac0inv at all quadrature points.
111 const int nqp = integ_rule.GetNPoints();
112 Vector rho_vals(nqp);
113 for (int i = 0; i < nzones; i++)
114 {
115 rho0.GetValues(i, integ_rule, rho_vals);
116 ElementTransformation *T = h1_fes.GetElementTransformation(i);
117 for (int q = 0; q < nqp; q++)
118 {
119 const IntegrationPoint &ip = integ_rule.IntPoint(q);
120 T->SetIntPoint(&ip);
121
122 DenseMatrixInverse Jinv(T->Jacobian());
123 Jinv.GetInverseMatrix(quad_data.Jac0inv(i*nqp + q));
124
125 const double rho0DetJ0 = T->Weight() * rho_vals(q);
126 quad_data.rho0DetJ0w(i*nqp + q) = rho0DetJ0 *
127 integ_rule.IntPoint(q).weight;
128 }
129 }
130
131 // Initial local mesh size (assumes all mesh elements are of the same type).
132 double area = 0.0;
133 Mesh *m = H1FESpace.GetMesh();
134 for (int i = 0; i < nzones; i++) { area += m->GetElementVolume(i); }
135 switch (m->GetElementBaseGeometry(0))
136 {
137 case Geometry::SEGMENT:
138 quad_data.h0 = area / nzones; break;
139 case Geometry::SQUARE:
140 quad_data.h0 = sqrt(area / nzones); break;
141 case Geometry::TRIANGLE:
142 quad_data.h0 = sqrt(2.0 * area / nzones); break;
143 case Geometry::CUBE:
144 quad_data.h0 = pow(area / nzones, 1.0/3.0); break;
145 case Geometry::TETRAHEDRON:
146 quad_data.h0 = pow(6.0 * area / nzones, 1.0/3.0); break;
147 default: MFEM_ABORT("Unknown zone type!");
148 }
149 quad_data.h0 /= (double) H1FESpace.GetOrder(0);
150
151 ForceIntegrator *fi = new ForceIntegrator(quad_data);
152 fi->SetIntRule(&integ_rule);
153 Force.AddDomainIntegrator(fi);
154 // Make a dummy assembly to figure out the sparsity.
155 Force.Assemble(0);
156 Force.Finalize(0);
157
158 if (p_assembly)
159 {
160 tensors1D = new Tensors1D(H1FESpace.GetFE(0)->GetOrder(),
161 L2FESpace.GetFE(0)->GetOrder(),
162 int(floor(0.7 + pow(nqp, 1.0 / dim))));
163 evaluator = new FastEvaluator(H1FESpace);
164 }
165
166 locCG.SetOperator(locEMassPA);
167 locCG.iterative_mode = false;
168 locCG.SetRelTol(1e-8);
169 locCG.SetAbsTol(1e-8 * numeric_limits<double>::epsilon());
170 locCG.SetMaxIter(200);
171 locCG.SetPrintLevel(0);
172 }
173
Mult(const Vector & S,Vector & dS_dt) const174 void LagrangianHydroOperator::Mult(const Vector &S, Vector &dS_dt) const
175 {
176 dS_dt = 0.0;
177
178 // Make sure that the mesh positions correspond to the ones in S. This is
179 // needed only because some mfem time integrators don't update the solution
180 // vector at every intermediate stage (hence they don't change the mesh).
181 Vector* sptr = (Vector*) &S;
182 GridFunction x;
183 x.MakeRef(&H1FESpace, *sptr, 0);
184 H1FESpace.GetMesh()->NewNodes(x, false);
185
186 UpdateQuadratureData(S);
187
188 // The monolithic BlockVector stores the unknown fields as follows:
189 // - Position
190 // - Velocity
191 // - Specific Internal Energy
192
193 const int Vsize_l2 = L2FESpace.GetVSize();
194 const int Vsize_h1 = H1FESpace.GetVSize();
195
196 GridFunction v, e;
197 v.MakeRef(&H1FESpace, *sptr, Vsize_h1);
198 e.MakeRef(&L2FESpace, *sptr, Vsize_h1*2);
199
200 GridFunction dx, dv, de;
201 dx.MakeRef(&H1FESpace, dS_dt, 0);
202 dv.MakeRef(&H1FESpace, dS_dt, Vsize_h1);
203 de.MakeRef(&L2FESpace, dS_dt, Vsize_h1*2);
204
205 // Set dx_dt = v (explicit).
206 dx = v;
207
208 if (!p_assembly)
209 {
210 Force = 0.0;
211 Force.Assemble();
212 }
213
214 // Solve for velocity.
215 Vector one(Vsize_l2), rhs(Vsize_h1); one = 1.0;
216 if (p_assembly)
217 {
218 ForcePA.Mult(one, rhs); rhs.Neg();
219
220 // Partial assembly solve for each velocity component.
221 const int size = H1compFESpace.GetVSize();
222 for (int c = 0; c < dim; c++)
223 {
224 Vector rhs_c(rhs.GetData() + c*size, size),
225 dv_c(dv.GetData() + c*size, size);
226
227 Array<int> vdofs_marker, ess_vdofs;
228 Array<int> ess_bdr(H1FESpace.GetMesh()->bdr_attributes.Max());
229 // Attributes 1/2/3 correspond to fixed-x/y/z boundaries, i.e.,
230 // we must enforce v_x/y/z = 0 for the velocity components.
231 ess_bdr = 0; ess_bdr[c] = 1;
232 // Essential vdofs as if there's only one component.
233 H1compFESpace.GetEssentialVDofs(ess_bdr, vdofs_marker);
234 FiniteElementSpace::MarkerToList(vdofs_marker, ess_vdofs);
235
236 dv_c = 0.0;
237 VMassPA.EliminateRHS(ess_vdofs, rhs_c);
238
239 CGSolver cg;
240 cg.SetOperator(VMassPA);
241 cg.SetRelTol(1e-8);
242 cg.SetAbsTol(0.0);
243 cg.SetMaxIter(300);
244 cg.SetPrintLevel(0);
245 cg.Mult(rhs_c, dv_c);
246 }
247 }
248 else
249 {
250 Force.Mult(one, rhs); rhs.Neg();
251 SparseMatrix A;
252 Vector B, X;
253 dv = 0.0;
254 Mv.FormLinearSystem(ess_tdofs, dv, rhs, A, X, B);
255 CGSolver cg;
256 cg.SetOperator(A);
257 cg.SetRelTol(1e-8);
258 cg.SetAbsTol(0.0);
259 cg.SetMaxIter(200);
260 cg.SetPrintLevel(0);
261 cg.Mult(B, X);
262 Mv.RecoverFEMSolution(X, rhs, dv);
263 }
264
265 // Solve for energy, assemble the energy source if such exists.
266 LinearForm *e_source = NULL;
267 if (source_type == 1) // 2D Taylor-Green.
268 {
269 e_source = new LinearForm(&L2FESpace);
270 TaylorCoefficient coeff;
271 DomainLFIntegrator *d = new DomainLFIntegrator(coeff, &integ_rule);
272 e_source->AddDomainIntegrator(d);
273 e_source->Assemble();
274 }
275 Array<int> l2dofs;
276 Vector e_rhs(Vsize_l2), loc_rhs(l2dofs_cnt), loc_de(l2dofs_cnt);
277 if (p_assembly)
278 {
279 ForcePA.MultTranspose(v, e_rhs);
280 if (e_source) { e_rhs += *e_source; }
281 for (int z = 0; z < nzones; z++)
282 {
283 L2FESpace.GetElementDofs(z, l2dofs);
284 e_rhs.GetSubVector(l2dofs, loc_rhs);
285 locEMassPA.SetZoneId(z);
286 locCG.Mult(loc_rhs, loc_de);
287 de.SetSubVector(l2dofs, loc_de);
288 }
289 }
290 else
291 {
292 Force.MultTranspose(v, e_rhs);
293 if (e_source) { e_rhs += *e_source; }
294 for (int z = 0; z < nzones; z++)
295 {
296 L2FESpace.GetElementDofs(z, l2dofs);
297 e_rhs.GetSubVector(l2dofs, loc_rhs);
298 Me_inv(z).Mult(loc_rhs, loc_de);
299 de.SetSubVector(l2dofs, loc_de);
300 }
301 }
302 delete e_source;
303
304 quad_data_is_current = false;
305 }
306
GetTimeStepEstimate(const Vector & S) const307 double LagrangianHydroOperator::GetTimeStepEstimate(const Vector &S) const
308 {
309 Vector* sptr = (Vector*) &S;
310 GridFunction x;
311 x.MakeRef(&H1FESpace, *sptr, 0);
312 H1FESpace.GetMesh()->NewNodes(x, false);
313 UpdateQuadratureData(S);
314
315 return quad_data.dt_est;
316 }
317
ResetTimeStepEstimate() const318 void LagrangianHydroOperator::ResetTimeStepEstimate() const
319 {
320 quad_data.dt_est = numeric_limits<double>::infinity();
321 }
322
ComputeDensity(GridFunction & rho)323 void LagrangianHydroOperator::ComputeDensity(GridFunction &rho)
324 {
325 rho.SetSpace(&L2FESpace);
326
327 DenseMatrix Mrho(l2dofs_cnt);
328 Vector rhs(l2dofs_cnt), rho_z(l2dofs_cnt);
329 Array<int> dofs(l2dofs_cnt);
330 DenseMatrixInverse inv(&Mrho);
331 MassIntegrator mi(&integ_rule);
332 DensityIntegrator di(quad_data);
333 di.SetIntRule(&integ_rule);
334 for (int i = 0; i < nzones; i++)
335 {
336 di.AssembleRHSElementVect(*L2FESpace.GetFE(i),
337 *L2FESpace.GetElementTransformation(i), rhs);
338 mi.AssembleElementMatrix(*L2FESpace.GetFE(i),
339 *L2FESpace.GetElementTransformation(i), Mrho);
340 inv.Factor();
341 inv.Mult(rhs, rho_z);
342 L2FESpace.GetElementDofs(i, dofs);
343 rho.SetSubVector(dofs, rho_z);
344 }
345 }
346
~LagrangianHydroOperator()347 LagrangianHydroOperator::~LagrangianHydroOperator()
348 {
349 delete tensors1D;
350 }
351
UpdateQuadratureData(const Vector & S) const352 void LagrangianHydroOperator::UpdateQuadratureData(const Vector &S) const
353 {
354 if (quad_data_is_current) { return; }
355
356 const int nqp = integ_rule.GetNPoints();
357
358 GridFunction x, v, e;
359 Vector* sptr = (Vector*) &S;
360 x.MakeRef(&H1FESpace, *sptr, 0);
361 v.MakeRef(&H1FESpace, *sptr, H1FESpace.GetVSize());
362 e.MakeRef(&L2FESpace, *sptr, 2*H1FESpace.GetVSize());
363 Vector e_vals, e_loc(l2dofs_cnt), vector_vals(h1dofs_cnt * dim);
364 DenseMatrix Jpi(dim), sgrad_v(dim), Jinv(dim), stress(dim), stressJiT(dim),
365 vecvalMat(vector_vals.GetData(), h1dofs_cnt, dim);
366 DenseTensor grad_v_ref(dim, dim, nqp);
367 Array<int> L2dofs, H1dofs;
368
369 // Batched computations are needed, because hydrodynamic codes usually
370 // involve expensive computations of material properties. Although this
371 // miniapp uses simple EOS equations, we still want to represent the batched
372 // cycle structure.
373 int nzones_batch = 3;
374 const int nbatches = nzones / nzones_batch + 1; // +1 for the remainder.
375 int nqp_batch = nqp * nzones_batch;
376 double *gamma_b = new double[nqp_batch],
377 *rho_b = new double[nqp_batch],
378 *e_b = new double[nqp_batch],
379 *p_b = new double[nqp_batch],
380 *cs_b = new double[nqp_batch];
381 // Jacobians of reference->physical transformations for all quadrature points
382 // in the batch.
383 DenseTensor *Jpr_b = new DenseTensor[nqp_batch];
384 for (int b = 0; b < nbatches; b++)
385 {
386 int z_id = b * nzones_batch; // Global index over zones.
387 // The last batch might not be full.
388 if (z_id == nzones) { break; }
389 else if (z_id + nzones_batch > nzones)
390 {
391 nzones_batch = nzones - z_id;
392 nqp_batch = nqp * nzones_batch;
393 }
394
395 double min_detJ = numeric_limits<double>::infinity();
396 for (int z = 0; z < nzones_batch; z++)
397 {
398 ElementTransformation *T = H1FESpace.GetElementTransformation(z_id);
399 Jpr_b[z].SetSize(dim, dim, nqp);
400
401 if (p_assembly)
402 {
403 // Energy values at quadrature point.
404 L2FESpace.GetElementDofs(z_id, L2dofs);
405 e.GetSubVector(L2dofs, e_loc);
406 evaluator->GetL2Values(e_loc, e_vals);
407
408 // All reference->physical Jacobians at the quadrature points.
409 H1FESpace.GetElementVDofs(z_id, H1dofs);
410 x.GetSubVector(H1dofs, vector_vals);
411 evaluator->GetVectorGrad(vecvalMat, Jpr_b[z]);
412 }
413 else { e.GetValues(z_id, integ_rule, e_vals); }
414 for (int q = 0; q < nqp; q++)
415 {
416 const IntegrationPoint &ip = integ_rule.IntPoint(q);
417 T->SetIntPoint(&ip);
418 if (!p_assembly) { Jpr_b[z](q) = T->Jacobian(); }
419 const double detJ = Jpr_b[z](q).Det();
420 min_detJ = min(min_detJ, detJ);
421
422 const int idx = z * nqp + q;
423 if (material_pcf == NULL) { gamma_b[idx] = 5./3.; } // Ideal gas.
424 else { gamma_b[idx] = material_pcf->Eval(*T, ip); }
425 rho_b[idx] = quad_data.rho0DetJ0w(z_id*nqp + q) / detJ / ip.weight;
426 e_b[idx] = max(0.0, e_vals(q));
427 }
428 ++z_id;
429 }
430
431 // Batched computation of material properties.
432 ComputeMaterialProperties(nqp_batch, gamma_b, rho_b, e_b, p_b, cs_b);
433
434 z_id -= nzones_batch;
435 for (int z = 0; z < nzones_batch; z++)
436 {
437 ElementTransformation *T = H1FESpace.GetElementTransformation(z_id);
438 if (p_assembly)
439 {
440 // All reference->physical Jacobians at the quadrature points.
441 H1FESpace.GetElementVDofs(z_id, H1dofs);
442 v.GetSubVector(H1dofs, vector_vals);
443 evaluator->GetVectorGrad(vecvalMat, grad_v_ref);
444 }
445 for (int q = 0; q < nqp; q++)
446 {
447 const IntegrationPoint &ip = integ_rule.IntPoint(q);
448 T->SetIntPoint(&ip);
449 // Note that the Jacobian was already computed above. We've chosen
450 // not to store the Jacobians for all batched quadrature points.
451 const DenseMatrix &Jpr = Jpr_b[z](q);
452 CalcInverse(Jpr, Jinv);
453 const double detJ = Jpr.Det(), rho = rho_b[z*nqp + q],
454 p = p_b[z*nqp + q], sound_speed = cs_b[z*nqp + q];
455
456 stress = 0.0;
457 for (int d = 0; d < dim; d++) { stress(d, d) = -p; }
458
459 double visc_coeff = 0.0;
460 if (use_viscosity)
461 {
462 // Compression-based length scale at the point. The first
463 // eigenvector of the symmetric velocity gradient gives the
464 // direction of maximal compression. This is used to define the
465 // relative change of the initial length scale.
466 if (p_assembly)
467 {
468 mfem::Mult(grad_v_ref(q), Jinv, sgrad_v);
469 }
470 else
471 {
472 v.GetVectorGradient(*T, sgrad_v);
473 }
474 sgrad_v.Symmetrize();
475 double eig_val_data[3], eig_vec_data[9];
476 if (dim==1)
477 {
478 eig_val_data[0] = sgrad_v(0, 0);
479 eig_vec_data[0] = 1.;
480 }
481 else { sgrad_v.CalcEigenvalues(eig_val_data, eig_vec_data); }
482 Vector compr_dir(eig_vec_data, dim);
483 // Computes the initial->physical transformation Jacobian.
484 mfem::Mult(Jpr, quad_data.Jac0inv(z_id*nqp + q), Jpi);
485 Vector ph_dir(dim); Jpi.Mult(compr_dir, ph_dir);
486 // Change of the initial mesh size in the compression direction.
487 const double h = quad_data.h0 * ph_dir.Norml2() /
488 compr_dir.Norml2();
489
490 // Measure of maximal compression.
491 const double mu = eig_val_data[0];
492 visc_coeff = 2.0 * rho * h * h * fabs(mu);
493 if (mu < 0.0) { visc_coeff += 0.5 * rho * h * sound_speed; }
494 stress.Add(visc_coeff, sgrad_v);
495 }
496
497 // Time step estimate at the point. Here the more relevant length
498 // scale is related to the actual mesh deformation; we use the min
499 // singular value of the ref->physical Jacobian. In addition, the
500 // time step estimate should be aware of the presence of shocks.
501 const double h_min =
502 Jpr.CalcSingularvalue(dim-1) / (double) H1FESpace.GetOrder(0);
503 const double inv_dt = sound_speed / h_min +
504 2.5 * visc_coeff / rho / h_min / h_min;
505 if (min_detJ < 0.0)
506 {
507 // This will force repetition of the step with smaller dt.
508 quad_data.dt_est = 0.0;
509 }
510 else
511 {
512 quad_data.dt_est = min(quad_data.dt_est, cfl * (1.0 / inv_dt) );
513 }
514
515 // Quadrature data for partial assembly of the force operator.
516 MultABt(stress, Jinv, stressJiT);
517 stressJiT *= integ_rule.IntPoint(q).weight * detJ;
518 for (int vd = 0 ; vd < dim; vd++)
519 {
520 for (int gd = 0; gd < dim; gd++)
521 {
522 quad_data.stressJinvT(vd)(z_id*nqp + q, gd) =
523 stressJiT(vd, gd);
524 }
525 }
526 }
527 ++z_id;
528 }
529 }
530 delete [] gamma_b;
531 delete [] rho_b;
532 delete [] e_b;
533 delete [] p_b;
534 delete [] cs_b;
535 delete [] Jpr_b;
536 quad_data_is_current = true;
537 }
538
539 } // namespace hydrodynamics
540
541 } // namespace mfem
542