48 #ifndef MFEM_USE_SUNDIALS
49 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
55 class ReducedSystemOperator;
82 ReducedSystemOperator *reduced_oper;
101 enum NonlinearSolverType
108 double visc,
double mu,
double K,
109 NonlinearSolverType nls_type);
116 virtual void ImplicitSolve(
const double dt,
const Vector &x,
Vector &k);
142 virtual int SUNImplicitSetup(
const Vector &y,
const Vector &fy,
143 int jok,
int *jcur,
double gamma);
147 virtual int SUNImplicitSolve(
const Vector &
b,
Vector &x,
double tol);
149 double ElasticEnergy(
const Vector &x)
const;
150 double KineticEnergy(
const Vector &v)
const;
153 virtual ~HyperelasticOperator();
160 class ReducedSystemOperator :
public Operator
174 void SetParameters(
double dt_,
const Vector *v_,
const Vector *x_);
182 virtual ~ReducedSystemOperator();
188 class ElasticEnergyCoefficient :
public Coefficient
197 : model(m), x(x_) { }
199 virtual ~ElasticEnergyCoefficient() { }
208 bool init_vis =
false);
211 int main(
int argc,
char *argv[])
214 const char *mesh_file =
"../../data/beam-quad.mesh";
217 int ode_solver_type = 3;
218 double t_final = 300.0;
223 bool visualization =
true;
224 const char *nls =
"newton";
228 const double reltol = 1e-1, abstol = 1e-1;
232 const double cvode_eps_lin = 1e-4;
234 const double arkode_eps_nonlin = 1e-6;
237 args.
AddOption(&mesh_file,
"-m",
"--mesh",
238 "Mesh file to use.");
239 args.
AddOption(&ref_levels,
"-r",
"--refine",
240 "Number of times to refine the mesh uniformly.");
242 "Order (degree) of the finite elements.");
243 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
245 "1 - Backward Euler,\n\t"
246 "2 - SDIRK2, L-stable\n\t"
247 "3 - SDIRK3, L-stable\n\t"
248 "4 - Implicit Midpoint,\n\t"
249 "5 - SDIRK2, A-stable,\n\t"
250 "6 - SDIRK3, A-stable,\n\t"
251 "7 - Forward Euler,\n\t"
255 "11 - CVODE implicit BDF, approximate Jacobian,\n\t"
256 "12 - CVODE implicit BDF, specified Jacobian,\n\t"
257 "13 - CVODE implicit ADAMS, approximate Jacobian,\n\t"
258 "14 - CVODE implicit ADAMS, specified Jacobian,\n\t"
259 "15 - ARKODE implicit, approximate Jacobian,\n\t"
260 "16 - ARKODE implicit, specified Jacobian,\n\t"
261 "17 - ARKODE explicit, 4th order.");
262 args.
AddOption(&nls,
"-nls",
"--nonlinear-solver",
263 "Nonlinear systems solver: "
264 "\"newton\" (plain Newton) or \"kinsol\" (KINSOL).");
265 args.
AddOption(&t_final,
"-tf",
"--t-final",
266 "Final time; start time is 0.");
267 args.
AddOption(&dt,
"-dt",
"--time-step",
269 args.
AddOption(&visc,
"-v",
"--viscosity",
270 "Viscosity coefficient.");
271 args.
AddOption(&mu,
"-mu",
"--shear-modulus",
272 "Shear modulus in the Neo-Hookean hyperelastic model.");
273 args.
AddOption(&K,
"-K",
"--bulk-modulus",
274 "Bulk modulus in the Neo-Hookean hyperelastic model.");
275 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
276 "--no-visualization",
277 "Enable or disable GLVis visualization.");
278 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
279 "Visualize every n-th timestep.");
289 if (ode_solver_type < 1 || ode_solver_type > 17)
291 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
297 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
301 map<string,HyperelasticOperator::NonlinearSolverType> nls_map;
302 nls_map[
"newton"] = HyperelasticOperator::NEWTON;
303 nls_map[
"kinsol"] = HyperelasticOperator::KINSOL;
304 if (nls_map.find(nls) == nls_map.end())
306 cout <<
"Unknown type of nonlinear solver: " << nls << endl;
313 for (
int lev = 0; lev < ref_levels; lev++)
328 cout <<
"Number of velocity/deformation unknowns: " << fe_size << endl;
331 fe_offset[1] = fe_size;
332 fe_offset[2] = 2*fe_size;
361 HyperelasticOperator oper(fespace, ess_bdr, visc, mu, K, nls_map[nls]);
366 char vishost[] =
"localhost";
368 vis_v.
open(vishost, visport);
371 visualize(vis_v, mesh, &x, &v,
"Velocity",
true);
372 vis_w.
open(vishost, visport);
375 oper.GetElasticEnergyDensity(x, w);
377 visualize(vis_w, mesh, &x, &w,
"Elastic energy density",
true);
383 cout <<
"initial elastic energy (EE) = " << ee0 << endl;
384 cout <<
"initial kinetic energy (KE) = " << ke0 << endl;
385 cout <<
"initial total energy (TE) = " << (ee0 + ke0) << endl;
396 switch (ode_solver_type)
408 case 8: ode_solver =
new RK2Solver(0.5);
break;
410 case 10: ode_solver =
new RK4Solver;
break;
417 CVodeSetEpsLin(cvode->
GetMem(), cvode_eps_lin);
419 if (ode_solver_type == 11)
423 ode_solver = cvode;
break;
430 CVodeSetEpsLin(cvode->
GetMem(), cvode_eps_lin);
432 if (ode_solver_type == 13)
436 ode_solver = cvode;
break;
443 ARKStepSetNonlinConvCoef(arkode->
GetMem(), arkode_eps_nonlin);
445 if (ode_solver_type == 15)
449 ode_solver = arkode;
break;
456 ode_solver = arkode;
break;
460 if (ode_solver_type < 11) { ode_solver->
Init(oper); }
464 bool last_step =
false;
465 for (
int ti = 1; !last_step; ti++)
467 double dt_real = min(dt, t_final - t);
469 ode_solver->
Step(vx, t, dt_real);
471 last_step = (t >= t_final - 1e-8*dt);
473 if (last_step || (ti % vis_steps) == 0)
478 cout <<
"step " << ti <<
", t = " << t <<
", EE = " << ee <<
", KE = "
479 << ke <<
", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
482 else if (arkode) { arkode->
PrintInfo(); }
490 oper.GetElasticEnergyDensity(x, w);
503 ofstream mesh_ofs(
"deformed.mesh");
504 mesh_ofs.precision(8);
505 mesh->
Print(mesh_ofs);
507 ofstream velo_ofs(
"velocity.sol");
508 velo_ofs.precision(8);
510 ofstream ee_ofs(
"elastic_energy.sol");
512 oper.GetElasticEnergyDensity(x, w);
525 GridFunction *field,
const char *field_name,
bool init_vis)
537 out <<
"solution\n" << *mesh << *field;
543 out <<
"window_size 800 800\n";
544 out <<
"window_title '" << field_name <<
"'\n";
551 out <<
"autoscale value\n";
558 ReducedSystemOperator::ReducedSystemOperator(
560 :
Operator(M_->Height()), M(M_), S(S_), H(H_), Jacobian(NULL),
561 dt(0.0), v(NULL), x(NULL), w(height), z(height)
564 void ReducedSystemOperator::SetParameters(
double dt_,
const Vector *v_,
567 dt = dt_; v = v_; x = x_;
580 Operator &ReducedSystemOperator::GetGradient(
const Vector &k)
const
583 Jacobian =
Add(1.0, M->SpMat(), dt, S->SpMat());
587 Jacobian->
Add(dt*dt, *grad_H);
591 ReducedSystemOperator::~ReducedSystemOperator()
600 NonlinearSolverType nls_type)
602 M(&fespace), S(&fespace), H(&fespace),
603 viscosity(visc), z(height/2),
604 grad_H(NULL), Jacobian(NULL)
606 const double rel_tol = 1e-8;
607 const int skip_zero_entries = 0;
609 const double ref_density = 1.0;
612 M.Assemble(skip_zero_entries);
614 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
616 M.FormSystemMatrix(ess_tdof_list, tmp);
618 M_solver.iterative_mode =
false;
619 M_solver.SetRelTol(rel_tol);
620 M_solver.SetAbsTol(0.0);
621 M_solver.SetMaxIter(30);
622 M_solver.SetPrintLevel(0);
623 M_solver.SetPreconditioner(M_prec);
624 M_solver.SetOperator(M.SpMat());
628 H.SetEssentialTrueDofs(ess_tdof_list);
632 S.Assemble(skip_zero_entries);
633 S.FormSystemMatrix(ess_tdof_list, tmp);
635 reduced_oper =
new ReducedSystemOperator(&M, &S, &H);
637 #ifndef MFEM_USE_SUITESPARSE
651 if (nls_type == KINSOL)
654 newton_solver = kinsolver;
656 newton_solver->SetMaxIter(200);
657 newton_solver->SetRelTol(rel_tol);
658 newton_solver->SetPrintLevel(0);
664 newton_solver->SetOperator(*reduced_oper);
665 newton_solver->SetMaxIter(10);
666 newton_solver->SetRelTol(rel_tol);
667 newton_solver->SetPrintLevel(-1);
669 newton_solver->SetSolver(*J_solver);
670 newton_solver->iterative_mode =
false;
683 if (viscosity != 0.0)
688 M_solver.Mult(z, dv_dt);
693 void HyperelasticOperator::ImplicitSolve(
const double dt,
708 reduced_oper->SetParameters(dt, &v, &x);
710 newton_solver->Mult(zero, dv_dt);
711 MFEM_VERIFY(newton_solver->GetConverged(),
712 "Nonlinear solver did not converge.");
714 cout <<
" num nonlin sol iters = " << newton_solver->GetNumIterations()
715 <<
", final norm = " << newton_solver->GetFinalNorm() <<
'\n';
717 add(v, dt, dv_dt, dx_dt);
720 int HyperelasticOperator::SUNImplicitSetup(
const Vector &y,
721 const Vector &fy,
int jok,
int *jcur,
724 int sc = y.
Size() / 2;
728 if (Jacobian) {
delete Jacobian; }
729 Jacobian =
Add(1.0, M.SpMat(), gamma, S.SpMat());
730 grad_H =
dynamic_cast<SparseMatrix *
>(&H.GetGradient(x));
731 Jacobian->
Add(gamma * gamma, *grad_H);
734 J_solver->SetOperator(*Jacobian);
746 int HyperelasticOperator::SUNImplicitSolve(
const Vector &
b,
Vector &x,
749 int sc = b.
Size() / 2;
757 grad_H->Mult(b_x, rhs);
761 J_solver->iterative_mode =
false;
762 J_solver->Mult(rhs, x_v);
764 add(b_x, saved_gamma, x_v, x_x);
769 double HyperelasticOperator::ElasticEnergy(
const Vector &x)
const
771 return H.GetEnergy(x);
774 double HyperelasticOperator::KineticEnergy(
const Vector &v)
const
776 return 0.5*M.InnerProduct(v, v);
779 void HyperelasticOperator::GetElasticEnergyDensity(
782 ElasticEnergyCoefficient w_coeff(*model, x);
786 HyperelasticOperator::~HyperelasticOperator()
789 delete newton_solver;
800 model.SetTransformation(T);
803 return model.EvalW(J)/J.Det();
817 const double s = 0.1/64.;
820 v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
void Init(TimeDependentOperator &f_)
Initialize CVODE: calls CVodeCreate() to create the CVODE memory and set some defaults.
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
double Eval(ElementTransformation &T, const IntegrationPoint &ip, double t)
Evaluate the coefficient in the element described by T at the point ip at time t. ...
void Add(const int i, const int j, const double a)
void InitialDeformation(const Vector &x, Vector &y)
virtual void Print(std::ostream &out=mfem::out) const
Conjugate gradient method.
Class for grid function - Vector with associated FE space.
Data type for scaled Jacobi-type smoother of sparse matrix.
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
A class to handle Vectors in a block fashion.
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Subclass constant coefficient.
Base abstract class for first order time dependent operators.
void SwapNodes(GridFunction *&nodes, int &own_nodes_)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
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.
void * GetMem() const
Access the SUNDIALS memory structure.
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)
void SetMaxStep(double dt_max)
Set the maximum time step.
int main(int argc, char *argv[])
Backward Euler ODE solver. L-stable.
double * GetData() const
Return a pointer to the beginning of the Vector data.
void add(const Vector &v1, const Vector &v2, Vector &v)
void SetMaxSetupCalls(int max_calls)
Set maximum number of nonlinear iterations without a Jacobian update.
void InitialVelocity(const Vector &x, Vector &v)
Interface to ARKode's ARKStep module – additive Runge-Kutta methods.
Direct sparse solver using UMFPACK.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
void SetPrintLevel(int print_lvl)
Mesh * GetMesh() const
Returns the mesh.
Interface to the CVODE library – linear multi-step methods.
Interface to the KINSOL library – nonlinear solver methods.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
void SetMaxIter(int max_it)
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
void SetMaxStep(double dt_max)
Set the maximum time step.
Newton's method for solving F(x)=b for a given operator F.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
void PrintUsage(std::ostream &out) const
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
int SpaceDimension() const
The classical explicit forth-order Runge-Kutta method, RK4.
void SetAbsTol(double atol)
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
void SetRelTol(double rtol)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
virtual void SetOperator(const Operator &op)
Set the nonlinear Operator of the system and initialize KINSOL.
Base class Coefficient that may optionally depend on time.
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)
Third-order, strong stability preserving (SSP) Runge-Kutta method.
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Implicit midpoint method. A-stable, not L-stable.
Class for integration point with weight.
void PrintInfo() const
Print various ARKStep statistics.
void PrintOptions(std::ostream &out) const
virtual void ProjectCoefficient(Coefficient &coeff)
Abstract class for hyperelastic models.
int open(const char hostname[], int port)
void GetNodes(Vector &node_coord) const
Arbitrary order H1-conforming (continuous) finite elements.
void PrintInfo() const
Print various CVODE statistics.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
The classical forward Euler method.
void Init(TimeDependentOperator &f_)
Initialize ARKode: calls ARKStepCreate() to create the ARKStep memory and set some defaults...
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
void UseSundialsLinearSolver()
Attach a SUNDIALS GMRES linear solver to ARKode.
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Arbitrary order "L2-conforming" discontinuous finite elements.
void Neg()
(*this) = -(*this)