56 #error This example requires that MFEM is built with MFEM_USE_HIOP=YES
83 class LinearScaleOperator :
public Operator
90 LinearScaleOperator(
const Vector &weight)
93 for (
int i = 0; i < width; i++) { grad(0, i) = w(i); }
109 class TanhSumOperator :
public Operator
115 TanhSumOperator(
int size) :
Operator(1, size), grad(1, width) { }
125 const double dtanh = 1.0 - ts * ts;
126 for (
int i = 0; i < width; i++) { grad(0, i) = dtanh; }
141 Vector massvec, d_lo, d_hi;
142 const LinearScaleOperator LSoper;
143 const TanhSumOperator TSoper;
146 OptimizedTransportProblem(
const Vector &xho,
const Vector &w,
double mass,
149 x_HO(xho), massvec(1), d_lo(1), d_hi(1),
150 LSoper(w), TSoper(w.Size())
154 SetEqualityConstraint(massvec);
159 MFEM_ASSERT(d_lo(0) < d_hi(0),
160 "The bounds produce an infeasible optimization problem");
161 SetInequalityConstraint(d_lo, d_hi);
163 SetSolutionBounds(xmin, xmax);
166 virtual double CalcObjective(
const Vector &x)
const
169 for (
int i = 0; i < input_size; i++)
171 const double d = x(i) - x_HO(i);
177 virtual void CalcObjectiveGrad(
const Vector &x,
Vector &grad)
const
179 for (
int i = 0; i < input_size; i++) { grad(i) = x(i) - x_HO(i); }
215 int main(
int argc,
char *argv[])
220 const char *mesh_file =
"../../data/periodic-hexagon.mesh";
223 int ode_solver_type = 3;
224 double t_final = 1.0;
226 bool visualization =
true;
232 cout.precision(precision);
235 args.
AddOption(&mesh_file,
"-m",
"--mesh",
236 "Mesh file to use.");
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.");
242 "Order (degree) of the finite elements.");
244 "Nonlinear optimizer: 1 - SLBQP,\n\t"
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",
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",
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.");
274 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
280 switch (ode_solver_type)
283 case 2: ode_solver =
new RK2Solver(1.0);
break;
285 case 4: ode_solver =
new RK4Solver;
break;
286 case 6: ode_solver =
new RK6Solver;
break;
288 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
297 for (
int lev = 0; lev < ref_levels; lev++)
312 cout <<
"Number of unknowns: " << fes.
GetVSize() << endl;
348 ofstream omesh(
"ex9.mesh");
349 omesh.precision(precision);
351 ofstream osol(
"ex9-init.gf");
352 osol.precision(precision);
363 #ifdef MFEM_USE_SIDRE
366 MFEM_ABORT(
"Must build with MFEM_USE_SIDRE=YES for binary output.");
385 sout.
open(vishost, visport);
388 cout <<
"Unable to connect to GLVis server at "
389 << vishost <<
':' << visport << endl;
390 visualization =
false;
391 cout <<
"GLVis visualization disabled.\n";
395 sout.precision(precision);
396 sout <<
"solution\n" << *mesh <<
u;
399 cout <<
"GLVis visualization paused."
400 <<
" Press space (in the GLVis window) to resume it.\n";
414 ode_solver->
Init(adv);
417 const double vol0 = M_rowsums *
u;
420 for (
int ti = 0; !done; )
422 double dt_real = min(dt, t_final - t);
423 adv.SetTimeStep(dt_real);
424 ode_solver->
Step(u, t, dt_real);
427 done = (t >= t_final - 1e-8*dt);
429 if (done || ti % vis_steps == 0)
431 cout <<
"time step: " << ti <<
", time: " << t << endl;
435 sout <<
"solution\n" << *mesh << u << flush;
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;
456 const double vol = M_rowsums *
u;
457 std::cout <<
"Vol error = " << vol - vol0 << endl;
462 ofstream osol(
"ex9-final.gf");
463 osol.precision(precision);
480 M(M_), K(K_),
b(b_), M_prec(), M_solver(), z(M_.Size()),
481 bf(bf_), M_rowsums(M_rs)
496 const int dofs = x.
Size();
497 Vector y_min(dofs), y_max(dofs);
499 for (
int i = 0, k = 0; i < dofs; i++)
503 for (
int end = In[i+1]; k < end; 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); }
512 for (
int i = 0; i < dofs; i++)
514 y_min(i) = (y_min(i) - x(i) ) / dt;
515 y_max(i) = (y_max(i) - x(i) ) / dt;
524 const double mass_y = 0.0;
528 const int max_iter = 500;
529 const double rtol = 1.e-7;
537 optsolver = tmp_opt_ptr;
539 MFEM_ABORT(
"MFEM is not built with HiOp support!");
551 OptimizedTransportProblem ot_prob(y, M_rowsums, mass_y, y_min, y_max);
558 optsolver->
Mult(y, y_out);
573 for (
int i = 0; i <
dim; i++)
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.);
597 const double w = M_PI/2;
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;
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.);
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;
630 for (
int i = 0; i <
dim; i++)
644 return (X(0) > -0.15 && X(0) < 0.15) ? 1.0 : 0.0;
649 double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
652 const double s = (1. + 0.25*
cos(2*M_PI*X(2)));
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;
663 double x_ = X(0), y_ = X(1), rho, phi;
670 const double f = M_PI;
671 return sin(f*X(0))*
sin(f*X(1));
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().
Conjugate gradient method.
Class for grid function - Vector with associated FE space.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
FDualNumber< tbase > tanh(const FDualNumber< tbase > &f)
tanh([dual number])
Data type for scaled Jacobi-type smoother of sparse matrix.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Base abstract class for first order time dependent operators.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
int * GetJ()
Return the array J.
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]...
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
int * GetI()
Return the array I.
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
void SetK(SparseMatrix &K_)
Data collection with Sidre routines following the Conduit mesh blueprint specification.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
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.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
virtual void Save()
Save the collection to disk.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
void SetTimeStep(double dt_)
void SetLinearConstraint(const Vector &w_, double a_)
void PrintUsage(std::ostream &out) const
Print the usage message.
void SetTime(double t)
Set physical time (for time-dependent simulations)
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
A general vector function coefficient.
The classical explicit forth-order Runge-Kutta method, RK4.
void SetAbsTol(double atol)
void SetRelTol(double rtol)
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
void velocity_function(const Vector &x, Vector &v)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Adapts the HiOp functionality to the MFEM OptimizationSolver interface.
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...
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Abstract solver for OptimizationProblems.
virtual void Print(std::ostream &os=mfem::out) const
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
virtual void Mult(const Vector &xt, Vector &x) const =0
Operator application: y=A(x).
void PrintOptions(std::ostream &out) const
Print the options.
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
double u0_function(const Vector &x)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
A general function coefficient.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
double u(const Vector &xvec)
The classical forward Euler method.
void SetBounds(const Vector &lo_, const Vector &hi_)
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
double Sum() const
Return the sum of the vector entries.
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
double inflow_function(const Vector &x)
void GetRowSums(Vector &x) const
For all i compute .
Arbitrary order "L2-conforming" discontinuous finite elements.
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.