MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex9p.cpp
Go to the documentation of this file.
1 // MFEM Example 9 - Parallel Version
2 // Nonlinear Constrained Optimization Modification
3 //
4 // Compile with: make ex9p
5 //
6 // Sample runs:
7 // mpirun -np 4 ex9p -m ../../data/periodic-segment.mesh -rs 3 -p 0 -o 2 -dt 0.002 -opt 1
8 // mpirun -np 4 ex9p -m ../../data/periodic-segment.mesh -rs 3 -p 0 -o 2 -dt 0.002 -opt 2
9 //
10 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 0 -rs 2 -dt 0.01 -tf 10 -opt 1
11 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 0 -rs 2 -dt 0.01 -tf 10 -opt 2
12 //
13 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 1 -rs 2 -dt 0.005 -tf 9 -opt 1
14 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 1 -rs 2 -dt 0.005 -tf 9 -opt 2
15 //
16 // mpirun -np 4 ex9p -m ../../data/amr-quad.mesh -p 1 -rs 1 -dt 0.002 -tf 9 -opt 1
17 // mpirun -np 4 ex9p -m ../../data/amr-quad.mesh -p 1 -rs 1 -dt 0.002 -tf 9 -opt 2
18 //
19 // mpirun -np 4 ex9p -m ../../data/disc-nurbs.mesh -p 1 -rs 2 -dt 0.005 -tf 9 -opt 1
20 // mpirun -np 4 ex9p -m ../../data/disc-nurbs.mesh -p 1 -rs 2 -dt 0.005 -tf 9 -opt 2
21 //
22 // mpirun -np 4 ex9p -m ../../data/disc-nurbs.mesh -p 2 -rs 2 -dt 0.01 -tf 9 -opt 1
23 // mpirun -np 4 ex9p -m ../../data/disc-nurbs.mesh -p 2 -rs 2 -dt 0.01 -tf 9 -opt 2
24 //
25 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 3 -rs 3 -dt 0.0025 -tf 9 -opt 1
26 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 3 -rs 3 -dt 0.0025 -tf 9 -opt 2
27 //
28 // mpirun -np 4 ex9p -m ../../data/periodic-cube.mesh -p 0 -rs 2 -o 2 -dt 0.02 -tf 8 -opt 1
29 // mpirun -np 4 ex9p -m ../../data/periodic-cube.mesh -p 0 -rs 2 -o 2 -dt 0.02 -tf 8 -opt 2
30 
31 // Description: This example modifies the standard MFEM ex9 by adding nonlinear
32 // constrained optimization capabilities through the SLBQP and
33 // HIOP solvers. It demonstrates how a user can define a custom
34 // class OptimizationProblem that includes linear/nonlinear
35 // equality/inequality constraints. This optimization is applied
36 // as post-processing to the solution of the transport equation.
37 //
38 // Description of ex9:
39 // This example code solves the time-dependent advection equation
40 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
41 // u0(x)=u(0,x) is a given initial condition.
42 //
43 // The example demonstrates the use of Discontinuous Galerkin (DG)
44 // bilinear forms in MFEM (face integrators), the use of explicit
45 // ODE time integrators, the definition of periodic boundary
46 // conditions through periodic meshes, as well as the use of GLVis
47 // for persistent visualization of a time-evolving solution. The
48 // saving of time-dependent data files for external visualization
49 // with VisIt (visit.llnl.gov) is also illustrated.
50 
51 #include "mfem.hpp"
52 #include <fstream>
53 #include <iostream>
54 
55 #ifndef MFEM_USE_HIOP
56 #error This example requires that MFEM is built with MFEM_USE_HIOP=YES
57 #endif
58 
59 using namespace std;
60 using namespace mfem;
61 
62 // Choice for the problem setup. The fluid velocity, initial condition and
63 // inflow boundary condition are chosen based on this parameter.
64 int problem;
65 
66 // Nonlinear optimizer.
68 
69 // Velocity coefficient
70 bool invert_velocity = false;
71 void velocity_function(const Vector &x, Vector &v);
72 
73 // Initial condition
74 double u0_function(const Vector &x);
75 
76 // Inflow boundary condition
77 double inflow_function(const Vector &x);
78 
79 // Mesh bounding box
81 
82 /// Computes C(x) = sum w_i x_i, where w is a given Vector.
83 class LinearScaleOperator : public Operator
84 {
85 private:
87  // Local weights.
88  const Vector &w;
89  // Gradient for the tdofs.
90  mutable DenseMatrix grad;
91 
92 public:
93  LinearScaleOperator(ParFiniteElementSpace &space, const Vector &weight)
94  : Operator(1, space.TrueVSize()),
95  pfes(space), w(weight), grad(1, width)
96  {
97  Vector w_glob(width);
98  pfes.Dof_TrueDof_Matrix()->MultTranspose(w, w_glob);
99  for (int i = 0; i < width; i++) { grad(0, i) = w_glob(i); }
100  }
101 
102  virtual void Mult(const Vector &x, Vector &y) const
103  {
104  Vector x_loc(w.Size());
105  pfes.GetProlongationMatrix()->Mult(x, x_loc);
106  const double loc_res = w * x_loc;
107  MPI_Allreduce(&loc_res, &y(0), 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
108  }
109 
110  virtual Operator &GetGradient(const Vector &x) const
111  {
112  return grad;
113  }
114 };
115 
116 /// Nonlinear monotone bounded operator to test nonlinear inequality constraints
117 /// Computes D(x) = tanh(sum(x_i)).
118 class TanhSumOperator : public Operator
119 {
120 private:
121  // Gradient for the tdofs.
122  mutable DenseMatrix grad;
123 
124 public:
125  TanhSumOperator(ParFiniteElementSpace &space)
126  : Operator(1, space.TrueVSize()), grad(1, width) { }
127 
128  virtual void Mult(const Vector &x, Vector &y) const
129  {
130  double sum_loc = x.Sum();
131  MPI_Allreduce(&sum_loc, &y(0), 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
132  y(0) = std::tanh(y(0));
133  }
134 
135  virtual Operator &GetGradient(const Vector &x) const
136  {
137  double sum_loc = x.Sum();
138  double dtanh;
139  MPI_Allreduce(&sum_loc, &dtanh, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
140  dtanh = 1.0 - pow(std::tanh(dtanh), 2);
141 
142  for (int i = 0; i < width; i++) { grad(0, i) = dtanh; }
143  return grad;
144  }
145 };
146 
147 
148 /** Monotone and conservative a posteriori correction for transport solutions:
149  * Find x that minimizes 0.5 || x - x_HO ||^2, subject to
150  * sum w_i x_i = mass,
151  * tanh(sum(x_i_min)) <= tanh(sum(x_i)) <= tanh(sum(x_i_max)),
152  * x_i_min <= x_i <= x_i_max,
153  */
154 class OptimizedTransportProblem : public OptimizationProblem
155 {
156 private:
157  const Vector &x_HO;
158  Vector massvec, d_lo, d_hi;
159  const LinearScaleOperator LSoper;
160  const TanhSumOperator TSoper;
161 
162 public:
163  OptimizedTransportProblem(ParFiniteElementSpace &space,
164  const Vector &xho, const Vector &w, double mass,
165  const Vector &xmin, const Vector &xmax)
166  : OptimizationProblem(xho.Size(), NULL, NULL),
167  x_HO(xho), massvec(1), d_lo(1), d_hi(1),
168  LSoper(space, w), TSoper(space)
169  {
170  C = &LSoper;
171  massvec(0) = mass;
172  SetEqualityConstraint(massvec);
173 
174  D = &TSoper;
175  double lsums[2], gsums[2];
176  lsums[0] = xmin.Sum();
177  lsums[1] = xmax.Sum();
178  MPI_Allreduce(lsums, gsums, 2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
179  d_lo(0) = std::tanh(gsums[0]);
180  d_hi(0) = std::tanh(gsums[1]);
181  MFEM_ASSERT(d_lo(0) < d_hi(0),
182  "The bounds produce an infeasible optimization problem");
183  SetInequalityConstraint(d_lo, d_hi);
184 
185  SetSolutionBounds(xmin, xmax);
186  }
187 
188  virtual double CalcObjective(const Vector &x) const
189  {
190  double loc_res = 0.0;
191  for (int i = 0; i < input_size; i++)
192  {
193  const double d = x(i) - x_HO(i);
194  loc_res += d * d;
195  }
196  loc_res *= 0.5;
197  double res;
198  MPI_Allreduce(&loc_res, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
199  return res;
200  }
201 
202  virtual void CalcObjectiveGrad(const Vector &x, Vector &grad) const
203  {
204  for (int i = 0; i < input_size; i++) { grad(i) = x(i) - x_HO(i); }
205  }
206 };
207 
208 
209 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
210  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
211  and advection matrices, and b describes the flow on the boundary. This can
212  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
213  used to evaluate the right-hand side. */
214 class FE_Evolution : public TimeDependentOperator
215 {
216 private:
217  HypreParMatrix &M, &K;
218  const Vector &b;
219  HypreSmoother M_prec;
220  CGSolver M_solver;
221 
222  mutable Vector z;
223 
224  double dt;
225  ParBilinearForm &pbf;
226  Vector &M_rowsums;
227 
228 public:
230  const Vector &b_, ParBilinearForm &pbf_, Vector &M_rs);
231 
232  void SetTimeStep(double dt_) { dt = dt_; }
233  void SetK(HypreParMatrix &K_) { K = K_; }
234  virtual void Mult(const Vector &x, Vector &y) const;
235 
236  virtual ~FE_Evolution() { }
237 };
238 
239 
240 int main(int argc, char *argv[])
241 {
242  // 1. Initialize MPI and HYPRE.
243  Mpi::Init(argc, argv);
244  int num_procs = Mpi::WorldSize();
245  int myid = Mpi::WorldRank();
246  Hypre::Init();
247 
248  // 2. Parse command-line options.
249  problem = 0;
250  optimizer_type = 2;
251  const char *mesh_file = "../../data/periodic-hexagon.mesh";
252  int ser_ref_levels = 2;
253  int par_ref_levels = 0;
254  int order = 3;
255  int ode_solver_type = 3;
256  double t_final = 1.0;
257  double dt = 0.01;
258  bool visualization = true;
259  bool visit = false;
260  bool binary = false;
261  int vis_steps = 5;
262 
263  int precision = 8;
264  cout.precision(precision);
265 
266  OptionsParser args(argc, argv);
267  args.AddOption(&mesh_file, "-m", "--mesh",
268  "Mesh file to use.");
269  args.AddOption(&problem, "-p", "--problem",
270  "Problem setup to use. See options in velocity_function().");
271  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
272  "Number of times to refine the mesh uniformly in serial.");
273  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
274  "Number of times to refine the mesh uniformly in parallel.");
275  args.AddOption(&order, "-o", "--order",
276  "Order (degree) of the finite elements.");
277  args.AddOption(&optimizer_type, "-opt", "--optimizer",
278  "Nonlinear optimizer: 1 - SLBQP,\n\t"
279  " 2 - HIOP.");
280  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
281  "ODE solver: 1 - Forward Euler,\n\t"
282  " 2 - RK2 SSP, 3 - RK3 SSP.");
283  args.AddOption(&t_final, "-tf", "--t-final",
284  "Final time; start time is 0.");
285  args.AddOption(&dt, "-dt", "--time-step",
286  "Time step.");
287  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
288  "--no-visualization",
289  "Enable or disable GLVis visualization.");
290  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
291  "--no-visit-datafiles",
292  "Save data files for VisIt (visit.llnl.gov) visualization.");
293  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
294  "--ascii-datafiles",
295  "Use binary (Sidre) or ascii format for VisIt data files.");
296  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
297  "Visualize every n-th timestep.");
298  args.Parse();
299  if (!args.Good())
300  {
301  if (myid == 0) { args.PrintUsage(cout); }
302  return 1;
303  }
304  if (myid == 0) { args.PrintOptions(cout); }
305 
306  // 3. Read the serial mesh from the given mesh file on all processors. We can
307  // handle geometrically periodic meshes in this code.
308  Mesh *mesh = new Mesh(mesh_file, 1, 1);
309  int dim = mesh->Dimension();
310 
311  // 4. Define the ODE solver used for time integration. Several explicit
312  // Runge-Kutta methods are available.
313  ODESolver *ode_solver = NULL;
314  switch (ode_solver_type)
315  {
316  case 1: ode_solver = new ForwardEulerSolver; break;
317  case 2: ode_solver = new RK2Solver(1.0); break;
318  case 3: ode_solver = new RK3SSPSolver; break;
319  case 4: ode_solver = new RK4Solver; break;
320  case 6: ode_solver = new RK6Solver; break;
321  default:
322  if (myid == 0)
323  {
324  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
325  }
326  delete mesh;
327  return 3;
328  }
329 
330  // 5. Refine the mesh in serial to increase the resolution. In this example
331  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
332  // a command-line parameter. If the mesh is of NURBS type, we convert it
333  // to a (piecewise-polynomial) high-order mesh.
334  for (int lev = 0; lev < ser_ref_levels; lev++)
335  {
336  mesh->UniformRefinement();
337  }
338  if (mesh->NURBSext)
339  {
340  mesh->SetCurvature(max(order, 1));
341  }
342  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
343 
344  // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
345  // this mesh further in parallel to increase the resolution. Once the
346  // parallel mesh is defined, the serial mesh can be deleted.
347  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
348  delete mesh;
349  for (int lev = 0; lev < par_ref_levels; lev++)
350  {
351  pmesh->UniformRefinement();
352  }
353 
354  // 7. Define the parallel discontinuous DG finite element space on the
355  // parallel refined mesh of the given polynomial order.
356  DG_FECollection fec(order, dim, BasisType::Positive);
357  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
358 
359  HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
360  if (myid == 0)
361  {
362  cout << "Number of unknowns: " << global_vSize << endl;
363  }
364 
365  // 8. Set up and assemble the parallel bilinear and linear forms (and the
366  // parallel hypre matrices) corresponding to the DG discretization. The
367  // DGTraceIntegrator involves integrals over mesh interior faces.
371 
372  ParBilinearForm *m = new ParBilinearForm(fes);
374  ParBilinearForm *k = new ParBilinearForm(fes);
375  k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
377  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
379  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
380 
381  ParLinearForm *b = new ParLinearForm(fes);
383  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
384 
385  m->Assemble();
386  m->Finalize();
387  int skip_zeros = 0;
388  k->Assemble(skip_zeros);
389  k->Finalize(skip_zeros);
390  b->Assemble();
391 
395 
396  // 9. Define the initial conditions, save the corresponding grid function to
397  // a file and (optionally) save data in the VisIt format and initialize
398  // GLVis visualization.
399  ParGridFunction *u = new ParGridFunction(fes);
400  u->ProjectCoefficient(u0);
401  HypreParVector *U = u->GetTrueDofs();
402 
403  {
404  ostringstream mesh_name, sol_name;
405  mesh_name << "ex9-mesh." << setfill('0') << setw(6) << myid;
406  sol_name << "ex9-init." << setfill('0') << setw(6) << myid;
407  ofstream omesh(mesh_name.str().c_str());
408  omesh.precision(precision);
409  pmesh->Print(omesh);
410  ofstream osol(sol_name.str().c_str());
411  osol.precision(precision);
412  u->Save(osol);
413  }
414 
415  // Create data collection for solution output: either VisItDataCollection for
416  // ascii data files, or SidreDataCollection for binary data files.
417  DataCollection *dc = NULL;
418  if (visit)
419  {
420  if (binary)
421  {
422 #ifdef MFEM_USE_SIDRE
423  dc = new SidreDataCollection("Example9-Parallel", pmesh);
424 #else
425  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
426 #endif
427  }
428  else
429  {
430  dc = new VisItDataCollection("Example9-Parallel", pmesh);
431  dc->SetPrecision(precision);
432  // To save the mesh using MFEM's parallel mesh format:
433  // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
434  }
435  dc->RegisterField("solution", u);
436  dc->SetCycle(0);
437  dc->SetTime(0.0);
438  dc->Save();
439  }
440 
441  socketstream sout;
442  if (visualization)
443  {
444  char vishost[] = "localhost";
445  int visport = 19916;
446  sout.open(vishost, visport);
447  if (!sout)
448  {
449  if (myid == 0)
450  cout << "Unable to connect to GLVis server at "
451  << vishost << ':' << visport << endl;
452  visualization = false;
453  if (myid == 0)
454  {
455  cout << "GLVis visualization disabled.\n";
456  }
457  }
458  else
459  {
460  sout << "parallel " << num_procs << " " << myid << "\n";
461  sout.precision(precision);
462  sout << "solution\n" << *pmesh << *u;
463  sout << "pause\n";
464  sout << flush;
465  if (myid == 0)
466  cout << "GLVis visualization paused."
467  << " Press space (in the GLVis window) to resume it.\n";
468  }
469  }
470 
471  Vector M_rowsums(m->Size());
472  m->SpMat().GetRowSums(M_rowsums);
473 
474  // 10. Define the time-dependent evolution operator describing the ODE
475  // right-hand side, and perform time-integration (looping over the time
476  // iterations, ti, with a time-step dt).
477  FE_Evolution adv(*M, *K, *B, *k, M_rowsums);
478 
479  double t = 0.0;
480  adv.SetTime(t);
481  ode_solver->Init(adv);
482 
483  *u = *U;
484 
485  // Compute initial volume.
486  const double vol0_loc = M_rowsums * (*u);
487  double vol0;
488  MPI_Allreduce(&vol0_loc, &vol0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
489 
490  bool done = false;
491  for (int ti = 0; !done; )
492  {
493  double dt_real = min(dt, t_final - t);
494  adv.SetTimeStep(dt_real);
495  ode_solver->Step(*U, t, dt_real);
496  ti++;
497 
498  done = (t >= t_final - 1e-8*dt);
499 
500  if (done || ti % vis_steps == 0)
501  {
502  if (myid == 0)
503  {
504  cout << "time step: " << ti << ", time: " << t << endl;
505  }
506 
507  // 11. Extract the parallel grid function corresponding to the finite
508  // element approximation U (the local solution on each processor).
509  *u = *U;
510 
511  if (visualization)
512  {
513  sout << "parallel " << num_procs << " " << myid << "\n";
514  sout << "solution\n" << *pmesh << *u << flush;
515  }
516 
517  if (visit)
518  {
519  dc->SetCycle(ti);
520  dc->SetTime(t);
521  dc->Save();
522  }
523  }
524  }
525 
526  // Print the error vs exact solution.
527  const double max_error = u->ComputeMaxError(u0),
528  l1_error = u->ComputeL1Error(u0),
529  l2_error = u->ComputeL2Error(u0);
530  if (myid == 0)
531  {
532  std::cout << "Linf error = " << max_error << endl
533  << "L1 error = " << l1_error << endl
534  << "L2 error = " << l2_error << endl;
535  }
536 
537  // Print error in volume.
538  const double vol_loc = M_rowsums * (*u);
539  double vol;
540  MPI_Allreduce(&vol_loc, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
541  if (myid == 0)
542  {
543  std::cout << "Vol error = " << vol - vol0 << endl;
544  }
545 
546  // 12. Save the final solution in parallel. This output can be viewed later
547  // using GLVis: "glvis -np <np> -m ex9-mesh -g ex9-final".
548  {
549  *u = *U;
550  ostringstream sol_name;
551  sol_name << "ex9-final." << setfill('0') << setw(6) << myid;
552  ofstream osol(sol_name.str().c_str());
553  osol.precision(precision);
554  u->Save(osol);
555  }
556 
557  // 13. Free the used memory.
558  delete U;
559  delete u;
560  delete B;
561  delete b;
562  delete K;
563  delete k;
564  delete M;
565  delete m;
566  delete fes;
567  delete pmesh;
568  delete ode_solver;
569  delete dc;
570 
571  return 0;
572 }
573 
574 
575 // Implementation of class FE_Evolution
577  const Vector &b_, ParBilinearForm &pbf_,
578  Vector &M_rs)
579  : TimeDependentOperator(M_.Height()),
580  M(M_), K(K_), b(b_), M_solver(M.GetComm()), z(M_.Height()),
581  pbf(pbf_), M_rowsums(M_rs)
582 {
583  M_prec.SetType(HypreSmoother::Jacobi);
584  M_solver.SetPreconditioner(M_prec);
585  M_solver.SetOperator(M);
586 
587  M_solver.iterative_mode = false;
588  M_solver.SetRelTol(1e-9);
589  M_solver.SetAbsTol(0.0);
590  M_solver.SetMaxIter(100);
591  M_solver.SetPrintLevel(0);
592 }
593 
594 void FE_Evolution::Mult(const Vector &x, Vector &y) const
595 {
596  // Get values on the ldofs.
597  ParFiniteElementSpace *pfes = pbf.ParFESpace();
598  ParGridFunction x_gf(pfes);
599  pfes->GetProlongationMatrix()->Mult(x, x_gf);
600 
601  // Compute bounds y_min, y_max for y from from x on the ldofs.
602  const int ldofs = x_gf.Size();
603  Vector y_min(ldofs), y_max(ldofs);
604  x_gf.ExchangeFaceNbrData();
605  Vector &x_nd = x_gf.FaceNbrData();
606  const int *In = pbf.SpMat().GetI(), *Jn = pbf.SpMat().GetJ();
607  for (int i = 0, k = 0; i < ldofs; i++)
608  {
609  double x_i_min = +std::numeric_limits<double>::infinity();
610  double x_i_max = -std::numeric_limits<double>::infinity();
611  for (int end = In[i+1]; k < end; k++)
612  {
613  const int j = Jn[k];
614  const double x_j = (j < ldofs) ? x(j): x_nd(j-ldofs);
615 
616  if (x_j > x_i_max) { x_i_max = x_j; }
617  if (x_j < x_i_min) { x_i_min = x_j; }
618  }
619  y_min(i) = x_i_min;
620  y_max(i) = x_i_max;
621  }
622  for (int i = 0; i < ldofs; i++)
623  {
624  y_min(i) = (y_min(i) - x_gf(i) ) / dt;
625  y_max(i) = (y_max(i) - x_gf(i) ) / dt;
626  }
627  Vector y_min_tdofs(y.Size()), y_max_tdofs(y.Size());
628  // Move the bounds to the tdofs.
629  pfes->GetRestrictionMatrix()->Mult(y_min, y_min_tdofs);
630  pfes->GetRestrictionMatrix()->Mult(y_max, y_max_tdofs);
631 
632  // Compute the high-order solution y = M^{-1} (K x + b) on the tdofs.
633  K.Mult(x, z);
634  z += b;
635  M_solver.Mult(z, y);
636 
637  // The solution y is an increment; it should not introduce new mass.
638  const double mass_y = 0.0;
639 
640  // Perform optimization on the tdofs.
641  Vector y_out(y.Size());
642  const int max_iter = 500;
643  const double rtol = 1.e-7;
644  double atol = 1.e-7;
645 
646  OptimizationSolver* optsolver = NULL;
647  if (optimizer_type == 2)
648  {
649 #ifdef MFEM_USE_HIOP
650  HiopNlpOptimizer *tmp_opt_ptr = new HiopNlpOptimizer(MPI_COMM_WORLD);
651  optsolver = tmp_opt_ptr;
652 #else
653  MFEM_ABORT("MFEM is not built with HiOp support!");
654 #endif
655  }
656  else
657  {
658  SLBQPOptimizer *slbqp = new SLBQPOptimizer(MPI_COMM_WORLD);
659  slbqp->SetBounds(y_min_tdofs, y_max_tdofs);
660  slbqp->SetLinearConstraint(M_rowsums, mass_y);
661  atol = 1.e-15;
662  optsolver = slbqp;
663  }
664 
665  OptimizedTransportProblem ot_prob(*pfes, y, M_rowsums, mass_y,
666  y_min_tdofs, y_max_tdofs);
667  optsolver->SetOptimizationProblem(ot_prob);
668 
669  optsolver->SetMaxIter(max_iter);
670  optsolver->SetAbsTol(atol);
671  optsolver->SetRelTol(rtol);
672  optsolver->SetPrintLevel(0);
673  optsolver->Mult(y, y_out);
674 
675  y = y_out;
676 
677  delete optsolver;
678 }
679 
680 
681 // Velocity coefficient
682 void velocity_function(const Vector &x, Vector &v)
683 {
684  int dim = x.Size();
685 
686  // map to the reference [-1,1] domain
687  Vector X(dim);
688  for (int i = 0; i < dim; i++)
689  {
690  double center = (bb_min[i] + bb_max[i]) * 0.5;
691  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
692  }
693 
694  switch (problem)
695  {
696  case 0:
697  {
698  // Translations in 1D, 2D, and 3D
699  switch (dim)
700  {
701  case 1: v(0) = (invert_velocity) ? -1.0 : 1.0; break;
702  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
703  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
704  break;
705  }
706  break;
707  }
708  case 1:
709  case 2:
710  {
711  // Clockwise rotation in 2D around the origin
712  const double w = M_PI/2;
713  switch (dim)
714  {
715  case 1: v(0) = 1.0; break;
716  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
717  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
718  }
719  break;
720  }
721  case 3:
722  {
723  // Clockwise twisting rotation in 2D around the origin
724  const double w = M_PI/2;
725  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
726  d = d*d;
727  switch (dim)
728  {
729  case 1: v(0) = 1.0; break;
730  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
731  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
732  }
733  break;
734  }
735  }
736 }
737 
738 // Initial condition
739 double u0_function(const Vector &x)
740 {
741  int dim = x.Size();
742 
743  // map to the reference [-1,1] domain
744  Vector X(dim);
745  for (int i = 0; i < dim; i++)
746  {
747  double center = (bb_min[i] + bb_max[i]) * 0.5;
748  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
749  }
750 
751  switch (problem)
752  {
753  case 0:
754  case 1:
755  {
756  switch (dim)
757  {
758  case 1:
759  return (X(0) > -0.15 && X(0) < 0.15) ? 1.0 : 0.0;
760  //return exp(-40.*pow(X(0)-0.0,2));
761  case 2:
762  case 3:
763  {
764  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
765  if (dim == 3)
766  {
767  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
768  rx *= s;
769  ry *= s;
770  }
771  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
772  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
773  }
774  }
775  }
776  case 2:
777  {
778  double x_ = X(0), y_ = X(1), rho, phi;
779  rho = hypot(x_, y_);
780  phi = atan2(y_, x_);
781  return pow(sin(M_PI*rho),2)*sin(3*phi);
782  }
783  case 3:
784  {
785  const double f = M_PI;
786  return sin(f*X(0))*sin(f*X(1));
787  }
788  }
789  return 0.0;
790 }
791 
792 // Inflow boundary condition (zero for the problems considered in this example)
793 double inflow_function(const Vector &x)
794 {
795  switch (problem)
796  {
797  case 0:
798  case 1:
799  case 2:
800  case 3: return 0.0;
801  }
802  return 0.0;
803 }
Vector bb_max
Definition: ex9.cpp:67
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
Conjugate gradient method.
Definition: solvers.hpp:465
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1153
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:434
Base abstract class for first order time dependent operators.
Definition: operator.hpp:289
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void SetTime(const double t_)
Set the current time.
Definition: operator.hpp:333
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:659
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
Data collection with Sidre routines following the Conduit mesh blueprint specification.
int optimizer_type
Definition: ex9.cpp:67
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
virtual void Save()
Save the collection to disk.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
constexpr int visport
int Size() const
Get the size of the BilinearForm as a square matrix.
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5343
void SetTimeStep(double dt_)
Definition: ex9p.cpp:232
void SetLinearConstraint(const Vector &w_, double a_)
Definition: solvers.cpp:2362
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetK(HypreParMatrix &K_)
Definition: ex9p.cpp:233
Parallel smoothers in hypre.
Definition: hypre.hpp:962
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:85
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:67
A general vector function coefficient.
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
void SetAbsTol(double atol)
Definition: solvers.hpp:199
void SetRelTol(double rtol)
Definition: solvers.hpp:198
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:501
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:729
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
int problem
Definition: ex15.cpp:62
string space
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:149
HYPRE_Int HYPRE_BigInt
virtual ~FE_Evolution()
Definition: ex9p.cpp:236
Adapts the HiOp functionality to the MFEM OptimizationSolver interface.
Definition: hiop.hpp:180
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: ex18.hpp:107
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
Definition: ex18.hpp:81
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
Abstract solver for OptimizationProblems.
Definition: solvers.hpp:820
bool invert_velocity
Definition: ex9.cpp:70
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:389
virtual void Mult(const Vector &xt, Vector &x) const =0
Operator application: y=A(x).
int dim
Definition: ex24.cpp:53
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
Class for parallel bilinear form.
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
double u0_function(const Vector &x)
Definition: ex9.cpp:558
RefCoord t[3]
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
RefCoord s[3]
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4770
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The classical forward Euler method.
Definition: ode.hpp:116
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition: solvers.cpp:2356
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.hpp:834
Class for parallel meshes.
Definition: pmesh.hpp:32
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:920
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:95
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
int main()
double inflow_function(const Vector &x)
Definition: ex9.cpp:611
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
void GetRowSums(Vector &x) const
For all i compute .
Definition: sparsemat.cpp:1210
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
alpha (q . grad u, v)