1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 // 3D Taylor-Green vortex benchmark example at Re=1600
13 // Unsteady flow of a decaying vortex is computed and compared against a known,
14 // analytical solution.
15
16 #include "navier_solver.hpp"
17 #include <fstream>
18
19 using namespace mfem;
20 using namespace navier;
21
22 struct s_NavierContext
23 {
24 int element_subdivisions = 1;
25 int order = 4;
26 double kinvis = 1.0 / 1600.0;
27 double t_final = 10 * 1e-3;
28 double dt = 1e-3;
29 bool pa = true;
30 bool ni = false;
31 bool visualization = false;
32 bool checkres = false;
33 } ctx;
34
vel_tgv(const Vector & x,double t,Vector & u)35 void vel_tgv(const Vector &x, double t, Vector &u)
36 {
37 double xi = x(0);
38 double yi = x(1);
39 double zi = x(2);
40
41 u(0) = sin(xi) * cos(yi) * cos(zi);
42 u(1) = -cos(xi) * sin(yi) * cos(zi);
43 u(2) = 0.0;
44 }
45
46 class QuantitiesOfInterest
47 {
48 public:
QuantitiesOfInterest(ParMesh * pmesh)49 QuantitiesOfInterest(ParMesh *pmesh)
50 {
51 H1_FECollection h1fec(1);
52 ParFiniteElementSpace h1fes(pmesh, &h1fec);
53
54 onecoeff.constant = 1.0;
55 mass_lf = new ParLinearForm(&h1fes);
56 mass_lf->AddDomainIntegrator(new DomainLFIntegrator(onecoeff));
57 mass_lf->Assemble();
58
59 ParGridFunction one_gf(&h1fes);
60 one_gf.ProjectCoefficient(onecoeff);
61
62 volume = mass_lf->operator()(one_gf);
63 };
64
ComputeKineticEnergy(ParGridFunction & v)65 double ComputeKineticEnergy(ParGridFunction &v)
66 {
67 Vector velx, vely, velz;
68 double integ = 0.0;
69 const FiniteElement *fe;
70 ElementTransformation *T;
71 FiniteElementSpace *fes = v.FESpace();
72
73 for (int i = 0; i < fes->GetNE(); i++)
74 {
75 fe = fes->GetFE(i);
76 double intorder = 2 * fe->GetOrder();
77 const IntegrationRule *ir = &(
78 IntRules.Get(fe->GetGeomType(), intorder));
79
80 v.GetValues(i, *ir, velx, 1);
81 v.GetValues(i, *ir, vely, 2);
82 v.GetValues(i, *ir, velz, 3);
83
84 T = fes->GetElementTransformation(i);
85 for (int j = 0; j < ir->GetNPoints(); j++)
86 {
87 const IntegrationPoint &ip = ir->IntPoint(j);
88 T->SetIntPoint(&ip);
89
90 double vel2 = velx(j) * velx(j) + vely(j) * vely(j)
91 + velz(j) * velz(j);
92
93 integ += ip.weight * T->Weight() * vel2;
94 }
95 }
96
97 double global_integral = 0.0;
98 MPI_Allreduce(&integ,
99 &global_integral,
100 1,
101 MPI_DOUBLE,
102 MPI_SUM,
103 MPI_COMM_WORLD);
104
105 return 0.5 * global_integral / volume;
106 };
107
~QuantitiesOfInterest()108 ~QuantitiesOfInterest() { delete mass_lf; };
109
110 private:
111 ConstantCoefficient onecoeff;
112 ParLinearForm *mass_lf;
113 double volume;
114 };
115
116 template<typename T>
sq(T x)117 T sq(T x)
118 {
119 return x * x;
120 }
121
122 // Computes Q = 0.5*(tr(\nabla u)^2 - tr(\nabla u \cdot \nabla u))
ComputeQCriterion(ParGridFunction & u,ParGridFunction & q)123 void ComputeQCriterion(ParGridFunction &u, ParGridFunction &q)
124 {
125 FiniteElementSpace *v_fes = u.FESpace();
126 FiniteElementSpace *fes = q.FESpace();
127
128 // AccumulateAndCountZones
129 Array<int> zones_per_vdof;
130 zones_per_vdof.SetSize(fes->GetVSize());
131 zones_per_vdof = 0;
132
133 q = 0.0;
134
135 // Local interpolation
136 int elndofs;
137 Array<int> v_dofs, dofs;
138 Vector vals;
139 Vector loc_data;
140 int vdim = v_fes->GetVDim();
141 DenseMatrix grad_hat;
142 DenseMatrix dshape;
143 DenseMatrix grad;
144
145 for (int e = 0; e < fes->GetNE(); ++e)
146 {
147 fes->GetElementVDofs(e, dofs);
148 v_fes->GetElementVDofs(e, v_dofs);
149 u.GetSubVector(v_dofs, loc_data);
150 vals.SetSize(dofs.Size());
151 ElementTransformation *tr = fes->GetElementTransformation(e);
152 const FiniteElement *el = fes->GetFE(e);
153 elndofs = el->GetDof();
154 int dim = el->GetDim();
155 dshape.SetSize(elndofs, dim);
156
157 for (int dof = 0; dof < elndofs; ++dof)
158 {
159 // Project
160 const IntegrationPoint &ip = el->GetNodes().IntPoint(dof);
161 tr->SetIntPoint(&ip);
162
163 // Eval
164 // GetVectorGradientHat
165 el->CalcDShape(tr->GetIntPoint(), dshape);
166 grad_hat.SetSize(vdim, dim);
167 DenseMatrix loc_data_mat(loc_data.GetData(), elndofs, vdim);
168 MultAtB(loc_data_mat, dshape, grad_hat);
169
170 const DenseMatrix &Jinv = tr->InverseJacobian();
171 grad.SetSize(grad_hat.Height(), Jinv.Width());
172 Mult(grad_hat, Jinv, grad);
173
174 double q_val = 0.5 * (sq(grad(0, 0)) + sq(grad(1, 1)) + sq(grad(2, 2)))
175 + grad(0, 1) * grad(1, 0) + grad(0, 2) * grad(2, 0)
176 + grad(1, 2) * grad(2, 1);
177
178 vals(dof) = q_val;
179 }
180
181 // Accumulate values in all dofs, count the zones.
182 for (int j = 0; j < dofs.Size(); j++)
183 {
184 int ldof = dofs[j];
185 q(ldof) += vals[j];
186 zones_per_vdof[ldof]++;
187 }
188 }
189
190 // Communication
191
192 // Count the zones globally.
193 GroupCommunicator &gcomm = q.ParFESpace()->GroupComm();
194 gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
195 gcomm.Bcast(zones_per_vdof);
196
197 // Accumulate for all vdofs.
198 gcomm.Reduce<double>(q.GetData(), GroupCommunicator::Sum);
199 gcomm.Bcast<double>(q.GetData());
200
201 // Compute means
202 for (int i = 0; i < q.Size(); i++)
203 {
204 const int nz = zones_per_vdof[i];
205 if (nz)
206 {
207 q(i) /= nz;
208 }
209 }
210 }
211
main(int argc,char * argv[])212 int main(int argc, char *argv[])
213 {
214 MPI_Session mpi(argc, argv);
215
216 OptionsParser args(argc, argv);
217 args.AddOption(&ctx.element_subdivisions,
218 "-es",
219 "--element-subdivisions",
220 "Number of 1d uniform subdivisions for each element.");
221 args.AddOption(&ctx.order,
222 "-o",
223 "--order",
224 "Order (degree) of the finite elements.");
225 args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
226 args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
227 args.AddOption(&ctx.pa,
228 "-pa",
229 "--enable-pa",
230 "-no-pa",
231 "--disable-pa",
232 "Enable partial assembly.");
233 args.AddOption(&ctx.ni,
234 "-ni",
235 "--enable-ni",
236 "-no-ni",
237 "--disable-ni",
238 "Enable numerical integration rules.");
239 args.AddOption(&ctx.visualization,
240 "-vis",
241 "--visualization",
242 "-no-vis",
243 "--no-visualization",
244 "Enable or disable GLVis visualization.");
245 args.AddOption(
246 &ctx.checkres,
247 "-cr",
248 "--checkresult",
249 "-no-cr",
250 "--no-checkresult",
251 "Enable or disable checking of the result. Returns -1 on failure.");
252 args.Parse();
253 if (!args.Good())
254 {
255 if (mpi.Root())
256 {
257 args.PrintUsage(mfem::out);
258 }
259 return 1;
260 }
261 if (mpi.Root())
262 {
263 args.PrintOptions(mfem::out);
264 }
265
266 Mesh orig_mesh("../../data/periodic-cube.mesh");
267 Mesh mesh = Mesh::MakeRefined(orig_mesh, ctx.element_subdivisions,
268 BasisType::ClosedUniform);
269 orig_mesh.Clear();
270
271 mesh.EnsureNodes();
272 GridFunction *nodes = mesh.GetNodes();
273 *nodes *= M_PI;
274
275 int nel = mesh.GetNE();
276 if (mpi.Root())
277 {
278 mfem::out << "Number of elements: " << nel << std::endl;
279 }
280
281 auto *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
282 mesh.Clear();
283
284 // Create the flow solver.
285 NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
286 flowsolver.EnablePA(ctx.pa);
287 flowsolver.EnableNI(ctx.ni);
288
289 // Set the initial condition.
290 ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
291 VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_tgv);
292 u_ic->ProjectCoefficient(u_excoeff);
293
294 double t = 0.0;
295 double dt = ctx.dt;
296 double t_final = ctx.t_final;
297 bool last_step = false;
298
299 flowsolver.Setup(dt);
300
301 ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
302 ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
303
304 ParGridFunction w_gf(*u_gf);
305 ParGridFunction q_gf(*p_gf);
306 flowsolver.ComputeCurl3D(*u_gf, w_gf);
307 ComputeQCriterion(*u_gf, q_gf);
308
309 QuantitiesOfInterest kin_energy(pmesh);
310
311 ParaViewDataCollection pvdc("shear_output", pmesh);
312 pvdc.SetDataFormat(VTKFormat::BINARY32);
313 pvdc.SetHighOrderOutput(true);
314 pvdc.SetLevelsOfDetail(ctx.order);
315 pvdc.SetCycle(0);
316 pvdc.SetTime(t);
317 pvdc.RegisterField("velocity", u_gf);
318 pvdc.RegisterField("pressure", p_gf);
319 pvdc.RegisterField("vorticity", &w_gf);
320 pvdc.RegisterField("qcriterion", &q_gf);
321 pvdc.Save();
322
323 double u_inf_loc = u_gf->Normlinf();
324 double p_inf_loc = p_gf->Normlinf();
325 double u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
326 double p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
327 double ke = kin_energy.ComputeKineticEnergy(*u_gf);
328
329 std::string fname = "tgv_out_p_" + std::to_string(ctx.order) + ".txt";
330 FILE *f = NULL;
331
332 if (mpi.Root())
333 {
334 int nel1d = std::round(pow(nel, 1.0 / 3.0));
335 int ngridpts = p_gf->ParFESpace()->GlobalVSize();
336 printf("%11s %11s %11s %11s %11s\n", "Time", "dt", "u_inf", "p_inf", "ke");
337 printf("%.5E %.5E %.5E %.5E %.5E\n", t, dt, u_inf, p_inf, ke);
338
339 f = fopen(fname.c_str(), "w");
340 fprintf(f, "3D Taylor Green Vortex\n");
341 fprintf(f, "order = %d\n", ctx.order);
342 fprintf(f, "grid = %d x %d x %d\n", nel1d, nel1d, nel1d);
343 fprintf(f, "dofs per component = %d\n", ngridpts);
344 fprintf(f, "=================================================\n");
345 fprintf(f, " time kinetic energy\n");
346 fprintf(f, "%20.16e %20.16e\n", t, ke);
347 fflush(f);
348 fflush(stdout);
349 }
350
351 for (int step = 0; !last_step; ++step)
352 {
353 if (t + dt >= t_final - dt / 2)
354 {
355 last_step = true;
356 }
357
358 flowsolver.Step(t, dt, step);
359
360 if ((step + 1) % 100 == 0 || last_step)
361 {
362 flowsolver.ComputeCurl3D(*u_gf, w_gf);
363 ComputeQCriterion(*u_gf, q_gf);
364 pvdc.SetCycle(step);
365 pvdc.SetTime(t);
366 pvdc.Save();
367 }
368
369 double u_inf_loc = u_gf->Normlinf();
370 double p_inf_loc = p_gf->Normlinf();
371 double u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
372 double p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
373 double ke = kin_energy.ComputeKineticEnergy(*u_gf);
374 if (mpi.Root())
375 {
376 printf("%.5E %.5E %.5E %.5E %.5E\n", t, dt, u_inf, p_inf, ke);
377 fprintf(f, "%20.16e %20.16e\n", t, ke);
378 fflush(f);
379 fflush(stdout);
380 }
381 }
382
383 flowsolver.PrintTimingData();
384
385 // Test if the result for the test run is as expected.
386 if (ctx.checkres)
387 {
388 double tol = 1e-5;
389 double ke_expected = 1.25e-1;
390 if (fabs(ke - ke_expected) > tol)
391 {
392 if (mpi.Root())
393 {
394 mfem::out << "Result has a larger error than expected."
395 << std::endl;
396 }
397 return -1;
398 }
399 }
400
401 delete pmesh;
402
403 return 0;
404 }
405