MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex9.cpp
Go to the documentation of this file.
1 // MFEM Example 9
2 // Nonlinear Constrained Optimization Modification
3 //
4 // Compile with: make ex9
5 //
6 // Sample runs:
7 // ex9 -m ../../data/periodic-segment.mesh -r 3 -p 0 -o 2 -dt 0.002 -opt 1
8 // ex9 -m ../../data/periodic-segment.mesh -r 3 -p 0 -o 2 -dt 0.002 -opt 2
9 //
10 // ex9 -m ../../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10 -opt 1
11 // ex9 -m ../../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10 -opt 2
12 //
13 // ex9 -m ../../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 1
14 // ex9 -m ../../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 2
15 //
16 // ex9 -m ../../data/amr-quad.mesh -p 1 -r 1 -dt 0.002 -tf 9 -opt 1
17 // ex9 -m ../../data/amr-quad.mesh -p 1 -r 1 -dt 0.002 -tf 9 -opt 2
18 //
19 // ex9 -m ../../data/disc-nurbs.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 1
20 // ex9 -m ../../data/disc-nurbs.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 2
21 //
22 // ex9 -m ../../data/disc-nurbs.mesh -p 2 -r 2 -dt 0.01 -tf 9 -opt 1
23 // ex9 -m ../../data/disc-nurbs.mesh -p 2 -r 2 -dt 0.01 -tf 9 -opt 2
24 //
25 // ex9 -m ../../data/periodic-square.mesh -p 3 -r 3 -dt 0.0025 -tf 9 -opt 1
26 // ex9 -m ../../data/periodic-square.mesh -p 3 -r 3 -dt 0.0025 -tf 9 -opt 2
27 //
28 // ex9 -m ../../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8 -opt 1
29 // ex9 -m ../../data/periodic-cube.mesh -p 0 -r 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:
86  const Vector &w;
87  mutable DenseMatrix grad;
88 
89 public:
90  LinearScaleOperator(const Vector &weight)
91  : Operator(1, weight.Size()), w(weight), grad(1, width)
92  {
93  for (int i = 0; i < width; i++) { grad(0, i) = w(i); }
94  }
95 
96  virtual void Mult(const Vector &x, Vector &y) const
97  {
98  y(0) = w * x;
99  }
100 
101  virtual Operator &GetGradient(const Vector &x) const
102  {
103  return grad;
104  }
105 };
106 
107 /// Nonlinear monotone bounded operator to test nonlinear inequality constraints
108 /// Computes D(x) = tanh(sum(x_i)).
109 class TanhSumOperator : public Operator
110 {
111 private:
112  mutable DenseMatrix grad;
113 
114 public:
115  TanhSumOperator(int size) : Operator(1, size), grad(1, width) { }
116 
117  virtual void Mult(const Vector &x, Vector &y) const
118  {
119  y(0) = std::tanh(x.Sum());
120  }
121 
122  virtual Operator &GetGradient(const Vector &x) const
123  {
124  const double ts = std::tanh(x.Sum());
125  const double dtanh = 1.0 - ts * ts;
126  for (int i = 0; i < width; i++) { grad(0, i) = dtanh; }
127  return grad;
128  }
129 };
130 
131 /** Monotone and conservative a posteriori correction for transport solutions:
132  * Find x that minimizes 0.5 || x - x_HO ||^2, subject to
133  * sum w_i x_i = mass,
134  * tanh(sum(x_i_min)) <= tanh(sum(x_i)) <= tanh(sum(x_i_max)),
135  * x_i_min <= x_i <= x_i_max,
136  */
137 class OptimizedTransportProblem : public OptimizationProblem
138 {
139 private:
140  const Vector &x_HO;
141  Vector massvec, d_lo, d_hi;
142  const LinearScaleOperator LSoper;
143  const TanhSumOperator TSoper;
144 
145 public:
146  OptimizedTransportProblem(const Vector &xho, const Vector &w, double mass,
147  const Vector &xmin, const Vector &xmax)
148  : OptimizationProblem(xho.Size(), NULL, NULL),
149  x_HO(xho), massvec(1), d_lo(1), d_hi(1),
150  LSoper(w), TSoper(w.Size())
151  {
152  C = &LSoper;
153  massvec(0) = mass;
154  SetEqualityConstraint(massvec);
155 
156  D = &TSoper;
157  d_lo(0) = std::tanh(xmin.Sum());
158  d_hi(0) = std::tanh(xmax.Sum());
159  MFEM_ASSERT(d_lo(0) < d_hi(0),
160  "The bounds produce an infeasible optimization problem");
161  SetInequalityConstraint(d_lo, d_hi);
162 
163  SetSolutionBounds(xmin, xmax);
164  }
165 
166  virtual double CalcObjective(const Vector &x) const
167  {
168  double res = 0.0;
169  for (int i = 0; i < input_size; i++)
170  {
171  const double d = x(i) - x_HO(i);
172  res += d * d;
173  }
174  return 0.5 * res;
175  }
176 
177  virtual void CalcObjectiveGrad(const Vector &x, Vector &grad) const
178  {
179  for (int i = 0; i < input_size; i++) { grad(i) = x(i) - x_HO(i); }
180  }
181 };
182 
183 
184 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
185  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
186  and advection matrices, and b describes the flow on the boundary. This can
187  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
188  used to evaluate the right-hand side. */
189 class FE_Evolution : public TimeDependentOperator
190 {
191 private:
192  SparseMatrix &M, &K;
193  const Vector &b;
194  DSmoother M_prec;
195  CGSolver M_solver;
196 
197  mutable Vector z;
198 
199  double dt;
200  BilinearForm &bf;
201  Vector &M_rowsums;
202 
203 public:
204  FE_Evolution(SparseMatrix &_M, SparseMatrix &_K, const Vector &_b,
205  BilinearForm &_bf, Vector &M_rs);
206 
207  void SetTimeStep(double _dt) { dt = _dt; }
208  void SetK(SparseMatrix &_K) { K = _K; }
209  virtual void Mult(const Vector &x, Vector &y) const;
210 
211  virtual ~FE_Evolution() { }
212 };
213 
214 
215 int main(int argc, char *argv[])
216 {
217  // 1. Parse command-line options.
218  problem = 0;
219  optimizer_type = 1;
220  const char *mesh_file = "../../data/periodic-hexagon.mesh";
221  int ref_levels = 2;
222  int order = 3;
223  int ode_solver_type = 3;
224  double t_final = 1.0;
225  double dt = 0.01;
226  bool visualization = true;
227  bool visit = false;
228  bool binary = false;
229  int vis_steps = 5;
230 
231  int precision = 8;
232  cout.precision(precision);
233 
234  OptionsParser args(argc, argv);
235  args.AddOption(&mesh_file, "-m", "--mesh",
236  "Mesh file to use.");
237  args.AddOption(&problem, "-p", "--problem",
238  "Problem setup to use. See options in velocity_function().");
239  args.AddOption(&ref_levels, "-r", "--refine",
240  "Number of times to refine the mesh uniformly.");
241  args.AddOption(&order, "-o", "--order",
242  "Order (degree) of the finite elements.");
243  args.AddOption(&optimizer_type, "-opt", "--optimizer",
244  "Nonlinear optimizer: 1 - SLBQP,\n\t"
245  " 2 - HIOP.");
246  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
247  "ODE solver: 1 - Forward Euler,\n\t"
248  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
249  args.AddOption(&t_final, "-tf", "--t-final",
250  "Final time; start time is 0.");
251  args.AddOption(&dt, "-dt", "--time-step",
252  "Time step.");
253  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
254  "--no-visualization",
255  "Enable or disable GLVis visualization.");
256  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
257  "--no-visit-datafiles",
258  "Save data files for VisIt (visit.llnl.gov) visualization.");
259  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
260  "--ascii-datafiles",
261  "Use binary (Sidre) or ascii format for VisIt data files.");
262  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
263  "Visualize every n-th timestep.");
264  args.Parse();
265  if (!args.Good())
266  {
267  args.PrintUsage(cout);
268  return 1;
269  }
270  args.PrintOptions(cout);
271 
272  // 2. Read the mesh from the given mesh file. We can handle geometrically
273  // periodic meshes in this code.
274  Mesh *mesh = new Mesh(mesh_file, 1, 1);
275  int dim = mesh->Dimension();
276 
277  // 3. Define the ODE solver used for time integration. Several explicit
278  // Runge-Kutta methods are available.
279  ODESolver *ode_solver = NULL;
280  switch (ode_solver_type)
281  {
282  case 1: ode_solver = new ForwardEulerSolver; break;
283  case 2: ode_solver = new RK2Solver(1.0); break;
284  case 3: ode_solver = new RK3SSPSolver; break;
285  case 4: ode_solver = new RK4Solver; break;
286  case 6: ode_solver = new RK6Solver; break;
287  default:
288  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
289  delete mesh;
290  return 3;
291  }
292 
293  // 4. Refine the mesh to increase the resolution. In this example we do
294  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
295  // command-line parameter. If the mesh is of NURBS type, we convert it to
296  // a (piecewise-polynomial) high-order mesh.
297  for (int lev = 0; lev < ref_levels; lev++)
298  {
299  mesh->UniformRefinement();
300  }
301  if (mesh->NURBSext)
302  {
303  mesh->SetCurvature(max(order, 1));
304  }
305  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
306 
307  // 5. Define the discontinuous DG finite element space of the given
308  // polynomial order on the refined mesh.
309  DG_FECollection fec(order, dim, BasisType::Positive);
310  FiniteElementSpace fes(mesh, &fec);
311 
312  cout << "Number of unknowns: " << fes.GetVSize() << endl;
313 
314  // 6. Set up and assemble the bilinear and linear forms corresponding to the
315  // DG discretization. The DGTraceIntegrator involves integrals over mesh
316  // interior faces.
320 
321  BilinearForm m(&fes);
323  BilinearForm k(&fes);
324  k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
326  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
328  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
329 
330  LinearForm b(&fes);
332  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
333 
334  m.Assemble();
335  m.Finalize();
336  int skip_zeros = 0;
337  k.Assemble(skip_zeros);
338  k.Finalize(skip_zeros);
339  b.Assemble();
340 
341  // 7. Define the initial conditions, save the corresponding grid function to
342  // a file and (optionally) save data in the VisIt format and initialize
343  // GLVis visualization.
344  GridFunction u(&fes);
345  u.ProjectCoefficient(u0);
346 
347  {
348  ofstream omesh("ex9.mesh");
349  omesh.precision(precision);
350  mesh->Print(omesh);
351  ofstream osol("ex9-init.gf");
352  osol.precision(precision);
353  u.Save(osol);
354  }
355 
356  // Create data collection for solution output: either VisItDataCollection for
357  // ascii data files, or SidreDataCollection for binary data files.
358  DataCollection *dc = NULL;
359  if (visit)
360  {
361  if (binary)
362  {
363 #ifdef MFEM_USE_SIDRE
364  dc = new SidreDataCollection("Example9", mesh);
365 #else
366  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
367 #endif
368  }
369  else
370  {
371  dc = new VisItDataCollection("Example9", mesh);
372  dc->SetPrecision(precision);
373  }
374  dc->RegisterField("solution", &u);
375  dc->SetCycle(0);
376  dc->SetTime(0.0);
377  dc->Save();
378  }
379 
380  socketstream sout;
381  if (visualization)
382  {
383  char vishost[] = "localhost";
384  int visport = 19916;
385  sout.open(vishost, visport);
386  if (!sout)
387  {
388  cout << "Unable to connect to GLVis server at "
389  << vishost << ':' << visport << endl;
390  visualization = false;
391  cout << "GLVis visualization disabled.\n";
392  }
393  else
394  {
395  sout.precision(precision);
396  sout << "solution\n" << *mesh << u;
397  sout << "pause\n";
398  sout << flush;
399  cout << "GLVis visualization paused."
400  << " Press space (in the GLVis window) to resume it.\n";
401  }
402  }
403 
404  Vector M_rowsums(m.Size());
405  m.SpMat().GetRowSums(M_rowsums);
406 
407  // 8. Define the time-dependent evolution operator describing the ODE
408  // right-hand side, and perform time-integration (looping over the time
409  // iterations, ti, with a time-step dt).
410  FE_Evolution adv(m.SpMat(), k.SpMat(), b, k, M_rowsums);
411 
412  double t = 0.0;
413  adv.SetTime(t);
414  ode_solver->Init(adv);
415 
416  // Compute initial volume.
417  const double vol0 = M_rowsums * u;
418 
419  bool done = false;
420  for (int ti = 0; !done; )
421  {
422  double dt_real = min(dt, t_final - t);
423  adv.SetTimeStep(dt_real);
424  ode_solver->Step(u, t, dt_real);
425  ti++;
426 
427  done = (t >= t_final - 1e-8*dt);
428 
429  if (done || ti % vis_steps == 0)
430  {
431  cout << "time step: " << ti << ", time: " << t << endl;
432 
433  if (visualization)
434  {
435  sout << "solution\n" << *mesh << u << flush;
436  }
437 
438  if (visit)
439  {
440  dc->SetCycle(ti);
441  dc->SetTime(t);
442  dc->Save();
443  }
444  }
445  }
446 
447  // Print the error vs exact solution.
448  const double max_error = u.ComputeMaxError(u0),
449  l1_error = u.ComputeL1Error(u0),
450  l2_error = u.ComputeL2Error(u0);
451  std::cout << "Linf error = " << max_error << endl
452  << "L1 error = " << l1_error << endl
453  << "L2 error = " << l2_error << endl;
454 
455  // Print error in volume.
456  const double vol = M_rowsums * u;
457  std::cout << "Vol error = " << vol - vol0 << endl;
458 
459  // 9. Save the final solution. This output can be viewed later using GLVis:
460  // "glvis -m ex9.mesh -g ex9-final.gf".
461  {
462  ofstream osol("ex9-final.gf");
463  osol.precision(precision);
464  u.Save(osol);
465  }
466 
467  // 10. Free the used memory.
468  delete ode_solver;
469  delete dc;
470  delete mesh;
471 
472  return 0;
473 }
474 
475 
476 // Implementation of class FE_Evolution
478  const Vector &_b, BilinearForm &_bf, Vector &M_rs)
479  : TimeDependentOperator(_M.Size()),
480  M(_M), K(_K), b(_b), M_prec(), M_solver(), z(_M.Size()),
481  bf(_bf), M_rowsums(M_rs)
482 {
483  M_solver.SetPreconditioner(M_prec);
484  M_solver.SetOperator(M);
485 
486  M_solver.iterative_mode = false;
487  M_solver.SetRelTol(1e-9);
488  M_solver.SetAbsTol(0.0);
489  M_solver.SetMaxIter(100);
490  M_solver.SetPrintLevel(0);
491 }
492 
493 void FE_Evolution::Mult(const Vector &x, Vector &y) const
494 {
495  // Compute bounds y_min, y_max for y from x on the ldofs.
496  const int dofs = x.Size();
497  Vector y_min(dofs), y_max(dofs);
498  const int *In = bf.SpMat().GetI(), *Jn = bf.SpMat().GetJ();
499  for (int i = 0, k = 0; i < dofs; i++)
500  {
501  double x_i_min = +std::numeric_limits<double>::infinity();
502  double x_i_max = -std::numeric_limits<double>::infinity();
503  for (int end = In[i+1]; k < end; k++)
504  {
505  const int j = Jn[k];
506  if (x(j) > x_i_max) { x_i_max = x(j); }
507  if (x(j) < x_i_min) { x_i_min = x(j); }
508  }
509  y_min(i) = x_i_min;
510  y_max(i) = x_i_max;
511  }
512  for (int i = 0; i < dofs; i++)
513  {
514  y_min(i) = (y_min(i) - x(i) ) / dt;
515  y_max(i) = (y_max(i) - x(i) ) / dt;
516  }
517 
518  // Compute the high-order solution y = M^{-1} (K x + b).
519  K.Mult(x, z);
520  z += b;
521  M_solver.Mult(z, y);
522 
523  // The solution y is an increment; it should not introduce new mass.
524  const double mass_y = 0.0;
525 
526  // Perform optimization.
527  Vector y_out(dofs);
528  const int max_iter = 500;
529  const double rtol = 1.e-7;
530  double atol = 1.e-7;
531 
532  OptimizationSolver *optsolver = NULL;
533  if (optimizer_type == 2)
534  {
535 #ifdef MFEM_USE_HIOP
536  HiopNlpOptimizer *tmp_opt_ptr = new HiopNlpOptimizer();
537  optsolver = tmp_opt_ptr;
538 #else
539  MFEM_ABORT("MFEM is not built with HiOp support!");
540 #endif
541  }
542  else
543  {
544  SLBQPOptimizer *slbqp = new SLBQPOptimizer();
545  slbqp->SetBounds(y_min, y_max);
546  slbqp->SetLinearConstraint(M_rowsums, mass_y);
547  atol = 1.e-15;
548  optsolver = slbqp;
549  }
550 
551  OptimizedTransportProblem ot_prob(y, M_rowsums, mass_y, y_min, y_max);
552  optsolver->SetOptimizationProblem(ot_prob);
553 
554  optsolver->SetMaxIter(max_iter);
555  optsolver->SetAbsTol(atol);
556  optsolver->SetRelTol(rtol);
557  optsolver->SetPrintLevel(0);
558  optsolver->Mult(y, y_out);
559 
560  y = y_out;
561 
562  delete optsolver;
563 }
564 
565 
566 // Velocity coefficient
567 void velocity_function(const Vector &x, Vector &v)
568 {
569  int dim = x.Size();
570 
571  // map to the reference [-1,1] domain
572  Vector X(dim);
573  for (int i = 0; i < dim; i++)
574  {
575  double center = (bb_min[i] + bb_max[i]) * 0.5;
576  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
577  }
578 
579  switch (problem)
580  {
581  case 0:
582  {
583  // Translations in 1D, 2D, and 3D
584  switch (dim)
585  {
586  case 1: v(0) = (invert_velocity) ? -1.0 : 1.0; break;
587  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
588  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
589  break;
590  }
591  break;
592  }
593  case 1:
594  case 2:
595  {
596  // Clockwise rotation in 2D around the origin
597  const double w = M_PI/2;
598  switch (dim)
599  {
600  case 1: v(0) = 1.0; break;
601  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
602  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
603  }
604  break;
605  }
606  case 3:
607  {
608  // Clockwise twisting rotation in 2D around the origin
609  const double w = M_PI/2;
610  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
611  d = d*d;
612  switch (dim)
613  {
614  case 1: v(0) = 1.0; break;
615  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
616  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
617  }
618  break;
619  }
620  }
621 }
622 
623 // Initial condition
624 double u0_function(const Vector &x)
625 {
626  int dim = x.Size();
627 
628  // map to the reference [-1,1] domain
629  Vector X(dim);
630  for (int i = 0; i < dim; i++)
631  {
632  double center = (bb_min[i] + bb_max[i]) * 0.5;
633  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
634  }
635 
636  switch (problem)
637  {
638  case 0:
639  case 1:
640  {
641  switch (dim)
642  {
643  case 1:
644  return (X(0) > -0.15 && X(0) < 0.15) ? 1.0 : 0.0;
645  //return exp(-40.*pow(X(0)-0.0,2));
646  case 2:
647  case 3:
648  {
649  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
650  if (dim == 3)
651  {
652  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
653  rx *= s;
654  ry *= s;
655  }
656  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
657  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
658  }
659  }
660  }
661  case 2:
662  {
663  double x_ = X(0), y_ = X(1), rho, phi;
664  rho = hypot(x_, y_);
665  phi = atan2(y_, x_);
666  return pow(sin(M_PI*rho),2)*sin(3*phi);
667  }
668  case 3:
669  {
670  const double f = M_PI;
671  return sin(f*X(0))*sin(f*X(1));
672  }
673  }
674  return 0.0;
675 }
676 
677 // Inflow boundary condition (zero for the problems considered in this example)
678 double inflow_function(const Vector &x)
679 {
680  switch (problem)
681  {
682  case 0:
683  case 1:
684  case 2:
685  case 3: return 0.0;
686  }
687  return 0.0;
688 }
Vector bb_max
Definition: ex9.cpp:64
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:400
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1234
Conjugate gradient method.
Definition: solvers.hpp:258
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:535
Base abstract class for first order time dependent operators.
Definition: operator.hpp:267
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:118
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:79
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:178
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]...
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:311
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
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:173
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:638
int main(int argc, char *argv[])
Definition: ex1.cpp:66
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:102
Data collection with Sidre routines following the Conduit mesh blueprint specification.
int optimizer_type
Definition: ex9.cpp:67
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3417
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetLinearConstraint(const Vector &_w, double _a)
Definition: solvers.cpp:1991
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
Data type sparse matrix.
Definition: sparsemat.hpp:46
constexpr char vishost[]
virtual void Save()
Save the collection to disk.
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
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:98
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4165
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:66
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
void SetTimeStep(double _dt)
Definition: ex9.cpp:207
Vector bb_min
Definition: ex9.cpp:64
A general vector function coefficient.
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
void SetAbsTol(double atol)
Definition: solvers.hpp:97
void SetRelTol(double rtol)
Definition: solvers.hpp:96
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:497
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
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:54
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:149
virtual ~FE_Evolution()
Definition: ex9.cpp:211
Adapts the HiOp functionality to the MFEM OptimizationSolver interface.
Definition: hiop.hpp:172
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:128
void SetBounds(const Vector &_lo, const Vector &_hi)
Definition: solvers.cpp:1985
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:203
Abstract solver for OptimizationProblems.
Definition: solvers.hpp:526
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
bool invert_velocity
Definition: ex9.cpp:70
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
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:304
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:45
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:2252
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:554
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:272
A general function coefficient.
void SetK(SparseMatrix &_K)
Definition: ex9.cpp:208
Vector data type.
Definition: vector.hpp:51
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix.
The classical forward Euler method.
Definition: ode.hpp:116
Abstract operator.
Definition: operator.hpp:24
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: solvers.hpp:540
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:849
double inflow_function(const Vector &x)
Definition: ex9.cpp:607
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:1020
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:221
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145
alpha (q . grad u, v)