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