78 virtual void ImplicitSolve(
const double dt,
const Vector &u,
Vector &k);
81 void SetParameters(
const Vector &u);
83 virtual ~ConductionOperator();
88 int main(
int argc,
char *argv[])
91 const char *mesh_file =
"../data/star.mesh";
94 int ode_solver_type = 3;
97 double alpha = 1.0e-2;
99 bool visualization =
true;
104 cout.precision(precision);
107 args.
AddOption(&mesh_file,
"-m",
"--mesh",
108 "Mesh file to use.");
109 args.
AddOption(&ref_levels,
"-r",
"--refine",
110 "Number of times to refine the mesh uniformly.");
112 "Order (degree) of the finite elements.");
113 args.
AddOption(&ode_solver_type,
"-s",
"--ode-solver",
114 "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
115 "\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
116 args.
AddOption(&t_final,
"-tf",
"--t-final",
117 "Final time; start time is 0.");
118 args.
AddOption(&dt,
"-dt",
"--time-step",
121 "Alpha coefficient.");
123 "Kappa coefficient offset.");
124 args.
AddOption(&visualization,
"-vis",
"--visualization",
"-no-vis",
125 "--no-visualization",
126 "Enable or disable GLVis visualization.");
127 args.
AddOption(&visit,
"-visit",
"--visit-datafiles",
"-no-visit",
128 "--no-visit-datafiles",
129 "Save data files for VisIt (visit.llnl.gov) visualization.");
130 args.
AddOption(&vis_steps,
"-vs",
"--visualization-steps",
131 "Visualize every n-th timestep.");
142 Mesh *mesh =
new Mesh(mesh_file, 1, 1);
149 switch (ode_solver_type)
157 case 12: ode_solver =
new RK2Solver(0.5);
break;
159 case 14: ode_solver =
new RK4Solver;
break;
166 cout <<
"Unknown ODE solver type: " << ode_solver_type <<
'\n';
174 for (
int lev = 0; lev < ref_levels; lev++)
185 cout <<
"Number of temperature unknowns: " << fe_size << endl;
197 ConductionOperator oper(fespace, alpha, kappa, u);
201 ofstream omesh(
"ex16.mesh");
202 omesh.precision(precision);
204 ofstream osol(
"ex16-init.gf");
205 osol.precision(precision);
221 char vishost[] =
"localhost";
223 sout.
open(vishost, visport);
226 cout <<
"Unable to connect to GLVis server at "
227 << vishost <<
':' << visport << endl;
228 visualization =
false;
229 cout <<
"GLVis visualization disabled.\n";
233 sout.precision(precision);
234 sout <<
"solution\n" << *mesh << u_gf;
237 cout <<
"GLVis visualization paused."
238 <<
" Press space (in the GLVis window) to resume it.\n";
244 ode_solver->
Init(oper);
247 bool last_step =
false;
248 for (
int ti = 1; !last_step; ti++)
250 if (t + dt >= t_final - dt/2)
255 ode_solver->
Step(u, t, dt);
257 if (last_step || (ti % vis_steps) == 0)
259 cout <<
"step " << ti <<
", t = " << t << endl;
264 sout <<
"solution\n" << *mesh << u_gf << flush;
274 oper.SetParameters(u);
280 ofstream osol(
"ex16-final.gf");
281 osol.precision(precision);
293 double kap,
const Vector &u)
295 T(NULL), current_dt(0.0), z(height)
297 const double rel_tol = 1e-8;
302 M->FormSystemMatrix(ess_tdof_list, Mmat);
304 M_solver.iterative_mode =
false;
305 M_solver.SetRelTol(rel_tol);
306 M_solver.SetAbsTol(0.0);
307 M_solver.SetMaxIter(30);
308 M_solver.SetPrintLevel(0);
309 M_solver.SetPreconditioner(M_prec);
310 M_solver.SetOperator(Mmat);
315 T_solver.iterative_mode =
false;
316 T_solver.SetRelTol(rel_tol);
317 T_solver.SetAbsTol(0.0);
318 T_solver.SetMaxIter(100);
319 T_solver.SetPrintLevel(0);
320 T_solver.SetPreconditioner(T_prec);
332 M_solver.Mult(z, du_dt);
335 void ConductionOperator::ImplicitSolve(
const double dt,
343 T =
Add(1.0, Mmat, dt, Kmat);
345 T_solver.SetOperator(*T);
347 MFEM_VERIFY(dt == current_dt,
"");
350 T_solver.Mult(z, du_dt);
353 void ConductionOperator::SetParameters(
const Vector &u)
356 u_alpha_gf.SetFromTrueDofs(u);
357 for (
int i = 0; i < u_alpha_gf.Size(); i++)
369 K->FormSystemMatrix(ess_tdof_list, Kmat);
374 ConductionOperator::~ConductionOperator()
virtual void Print(std::ostream &out=mfem::out) const
Conjugate gradient method.
Class for grid function - Vector with associated FE space.
double InitialTemperature(const Vector &x)
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Base abstract class for time dependent operators.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
double Norml2() const
Returns the l2 norm of the vector.
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]...
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
virtual void Save()
Save the collection and a VisIt root file.
int main(int argc, char *argv[])
Backward Euler ODE solver. L-stable.
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 UniformRefinement(int i, const DSTable &, int *, int *, int *)
Data collection with VisIt I/O routines.
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction. If all dofs are true, then tv will be set to point to th...
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
void PrintUsage(std::ostream &out) const
void SetTime(double t)
Set physical time (for time-dependent simulations)
The classical explicit forth-order Runge-Kutta method, RK4.
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)
Third-order, strong stability preserving (SSP) Runge-Kutta method.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
Implicit midpoint method. A-stable, not L-stable.
void PrintOptions(std::ostream &out) const
virtual void ProjectCoefficient(Coefficient &coeff)
int open(const char hostname[], int port)
class for C-function coefficient
Arbitrary order H1-conforming (continuous) finite elements.
The classical forward Euler method.
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.