NOTE - THIS EXAMPLE HAS NOT BEEN UPGRADED TO NEW WEAK FORMS YET
Git reference: Example 05-trilinos-coupled.
The purpose of this example is to show how to use Trilinos for nonlinear time-dependent coupled PDE systems. Solved by NOX solver via Newton or JFNK, with or without preconditioning.
We solve the simplified flame propagation problem from P03-timedep/flame.
The code is the same as in example 19 until the definition of the weak formulation, where we use diagonal blocks of the Jacobian for preconditioning:
// Initialize weak formulation.
WeakForm wf(2, JFNK ? true : false);
if (!JFNK || (JFNK && PRECOND == 1))
{
wf.add_matrix_form(0, 0, callback(newton_bilinear_form_0_0), HERMES_NONSYM, HERMES_ANY, &omega_dt);
wf.add_matrix_form_surf(0, 0, callback(newton_bilinear_form_0_0_surf), 3);
wf.add_matrix_form(1, 1, callback(newton_bilinear_form_1_1), HERMES_NONSYM, HERMES_ANY, &omega_dc);
wf.add_matrix_form(0, 1, callback(newton_bilinear_form_0_1), HERMES_NONSYM, HERMES_ANY, &omega_dc);
wf.add_matrix_form(1, 0, callback(newton_bilinear_form_1_0), HERMES_NONSYM, HERMES_ANY, &omega_dt);
}
else if (PRECOND == 2)
{
wf.add_matrix_form(0, 0, callback(precond_0_0));
wf.add_matrix_form(1, 1, callback(precond_1_1));
}
wf.add_vector_form(0, callback(newton_linear_form_0), HERMES_ANY,
Hermes::Tuple<MeshFunction*>(&t_prev_time_1, &t_prev_time_2, &omega));
wf.add_vector_form_surf(0, callback(newton_linear_form_0_surf), 3);
wf.add_vector_form(1, callback(newton_linear_form_1), HERMES_ANY,
Hermes::Tuple<MeshFunction*>(&c_prev_time_1, &c_prev_time_2, &omega));
Next we project the initial conditions to obtain an initial coefficient vector for NOX:
// Project the functions "t_prev_time_1" and "c_prev_time_1" on the FE space
// in order to obtain initial vector for NOX.
info("Projecting initial solutions on the FE meshes.");
scalar* coeff_vec = new scalar[ndof];
OGProjection::project_global(Hermes::Tuple<Space *>(t_space, c_space),
Hermes::Tuple<MeshFunction*>(&t_prev_time_1, &c_prev_time_1),
coeff_vec);
Then we initialize the DiscreteProblem class, NOX solver, and preconditioner:
// Initialize finite element problem.
DiscreteProblem dp(&wf, Hermes::Tuple<Space*>(t_space, c_space));
// Initialize NOX solver and preconditioner.
NoxSolver solver(&dp);
RCP<Precond> pc = rcp(new MlPrecond("sa"));
if (PRECOND)
{
if (JFNK) solver.set_precond(pc);
else solver.set_precond("Ifpack");
}
if (TRILINOS_OUTPUT)
solver.set_output_flags(NOX::Utils::Error | NOX::Utils::OuterIteration |
NOX::Utils::OuterIterationStatusTest |
NOX::Utils::LinearSolverDetails);
Output flags are set as follows:
if (TRILINOS_OUTPUT)
solver.set_output_flags(NOX::Utils::Error | NOX::Utils::OuterIteration |
NOX::Utils::OuterIterationStatusTest |
NOX::Utils::LinearSolverDetails);
The time stepping loop is as usual:
// Time stepping loop:
double total_time = 0.0;
cpu_time.tick_reset();
for (int ts = 1; total_time <= T_FINAL; ts++)
{
info("---- Time step %d, t = %g s", ts, total_time + TAU);
cpu_time.tick(HERMES_SKIP);
solver.set_init_sln(coeff_vec);
if (solver.solve())
{
Solution::vector_to_solutions(solver.get_solution(), Hermes::Tuple<Space *>(t_space, c_space),
Hermes::Tuple<Solution *>(&t_prev_newton, &c_prev_newton));
cpu_time.tick();
info("Number of nonlin iterations: %d (norm of residual: %g)",
solver.get_num_iters(), solver.get_residual());
info("Total number of iterations in linsolver: %d (achieved tolerance in the last step: %g)",
solver.get_num_lin_iters(), solver.get_achieved_tol());
// Time measurement.
cpu_time.tick(HERMES_SKIP);
// Visualization.
DXDYFilter omega_view(omega_fn, Hermes::Tuple<MeshFunction*>(&t_prev_newton, &c_prev_newton));
rview.set_min_max_range(0.0,2.0);
rview.show(&omega_view);
cpu_time.tick(HERMES_SKIP);
// Skip visualization time.
cpu_time.tick(HERMES_SKIP);
// Update global time.
total_time += TAU;
// Saving solutions for the next time step.
t_prev_time_2.copy(&t_prev_time_1);
c_prev_time_2.copy(&c_prev_time_1);
t_prev_time_1 = t_prev_newton;
c_prev_time_1 = c_prev_newton;
}
else
error("NOX failed.");