Git reference: Tutorial example 02-system.
We consider a simplified version of the Fitzhugh-Nagumo equation. This equation is a prominent example of activator-inhibitor systems in two-component reaction-diffusion equations, It describes a prototype of an excitable system (e.g., a neuron) and its stationary form is
-d^2_u \Delta u - f(u) + \sigma v - g_1 = 0,\\ -d^2_v \Delta v - u + v - g_2 = 0.
Here the unknowns u, v are the voltage and v-gate, respectively. The nonlinear function
f(u) = \lambda u - u^3 - \kappa
describes how an action potential travels through a nerve. Obviously this system is nonlinear.
In order to make it simpler for this tutorial, we replace the function f(u) with just u:
f(u) = u.
Our computational domain is the square (-1,1)^2 and we consider zero Dirichlet conditions for both u and v. In order to enable fair convergence comparisons, we will use the following functions as the exact solution:
u(x,y) = \cos\left(\frac{\pi}{2}x\right) \cos\left(\frac{\pi}{2}y\right),\\ v(x,y) = \hat u(x) \hat u(y)
where
\hat u(x) = 1 - \frac{e^{kx} + e^{-kx}}{e^k + e^{-k}}
is the exact solution of the one-dimensional singularly perturbed problem
-u'' + k^2 u - k^2 = 0
in (-1,1), equipped with zero Dirichlet boundary conditions.
The following two figures show the solutions u and v. Notice their large qualitative differences: While u is smooth in the entire domain, v has a thin boundary layer along the boundary:
The source functions g_1 and g_2 are obtained by inserting u and v into the PDE system. These functions are not extremely pretty, but they are not too bad either:
class CustomExactFunction1
{
public:
CustomExactFunction1()
{
};
double val(double x)
{
return cos(M_PI*x/2);
}
double dx(double x)
{
return -sin(M_PI*x/2)*(M_PI/2.);
}
double ddxx(double x)
{
return -cos(M_PI*x/2)*(M_PI/2.)*(M_PI/2.);
}
};
class CustomExactFunction2
{
public:
CustomExactFunction2(double K) : K(K)
{
};
double val(double x)
{
return 1. - (exp(K*x) + exp(-K*x))/(exp(K) + exp(-K));
}
double dx(double x)
{
return -K*(exp(K*x) - exp(-K*x))/(exp(K) + exp(-K));
}
double ddxx(double x)
{
return -K*K*(exp(K*x) + exp(-K*x))/(exp(K) + exp(-K));
}
double K;
};
The weak forms can be found in the file definitions.cpp and they are registered as follows:
Weak formulation comprises default and custom forms:
class CustomWeakForm : public WeakForm
{
public:
CustomWeakForm(CustomRightHandSide1* g1, CustomRightHandSide2* g2) : WeakForm(2)
{
// Jacobian.
add_matrix_form(new WeakFormsH1::DefaultJacobianDiffusion(0, 0, HERMES_ANY, new HermesFunction(D_u * D_u)));
add_matrix_form(new WeakFormsH1::DefaultMatrixFormVol(0, 0, HERMES_ANY, new HermesFunction(-1.0)));
add_matrix_form(new WeakFormsH1::DefaultMatrixFormVol(0, 1, HERMES_ANY, new HermesFunction(g1->sigma), HERMES_NONSYM));
add_matrix_form(new WeakFormsH1::DefaultMatrixFormVol(1, 0, HERMES_ANY, new HermesFunction(-1.0), HERMES_NONSYM));
add_matrix_form(new WeakFormsH1::DefaultJacobianDiffusion(1, 1, HERMES_ANY, new HermesFunction(D_v * D_v)));
add_matrix_form(new WeakFormsH1::DefaultMatrixFormVol(1, 1, HERMES_ANY, new HermesFunction(1.0)));
// Residual.
add_vector_form(new CustomResidual1(D_u, g1->sigma, g1));
add_vector_form(new CustomResidual2(D_v, g2));
}
};
Beware that although each of the forms is actually symmetric, one cannot use the HERMES_SYM flag as in the elasticity equations, since it has a slightly different meaning (see example P01-linear/08-system).
The adaptivity workflow is standard, first we construct the reference spaces:
// Construct globally refined reference mesh and setup reference space.
Hermes::vector<Space *>* ref_spaces =
Space::construct_refined_spaces(Hermes::vector<Space *>(&u_space, &v_space));
Then we initialize matrix solver:
// Initialize matrix solver.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);
Solve the discrete problem using the Newton’s method:
// Perform Newton's iteration.
bool jacobian_changed = true;
bool verbose = true;
if (!hermes2d.solve_newton(coeff_vec, &dp, solver, matrix, rhs, jacobian_changed,
1e-8, 100, verbose)) error("Newton's iteration failed.");
Translate the coefficient vector into the two Solutions:
// Translate the resulting coefficient vector into Solutions.
Solution::vector_to_solutions(coeff_vec, *ref_spaces, Hermes::vector<Solution *>(&u_ref_sln, &v_ref_sln));
Project reference solutions to the coarse meshes:
// Project the fine mesh solution onto the coarse mesh.
info("Projecting reference solution on coarse mesh.");
OGProjection::project_global(Hermes::vector<Space *>(&u_space, &v_space),
Hermes::vector<Solution *>(&u_ref_sln, &v_ref_sln),
Hermes::vector<Solution *>(&u_sln, &v_sln), matrix_solver);
Calculate error estimates:
// Calculate error estimate for each solution component and the total error estimate.
Hermes::vector<double> err_est_rel;
double err_est_rel_total = adaptivity->calc_err_est(Hermes::vector<Solution *>(&u_sln, &v_sln),
Hermes::vector<Solution *>(&u_ref_sln, &v_ref_sln),
&err_est_rel) * 100;
Calculate exact errors (optional):
// Calculate exact error for each solution component and the total exact error.
Hermes::vector<double> err_exact_rel;
bool solutions_for_adapt = false;
double err_exact_rel_total = adaptivity->calc_err_exact(Hermes::vector<Solution *>(&u_sln, &v_sln),
Hermes::vector<Solution *>(&exact_u, &exact_v),
&err_exact_rel, solutions_for_adapt) * 100;
Adapt the coarse meshes:
// If err_est too large, adapt the mesh.
if (err_est_rel_total < ERR_STOP)
done = true;
else
{
info("Adapting coarse mesh.");
done = adaptivity->adapt(Hermes::vector<RefinementSelectors::Selector *>(&selector, &selector),
THRESHOLD, STRATEGY, MESH_REGULARITY);
}
if (Space::get_num_dofs(Hermes::vector<Space *>(&u_space, &v_space)) >= NDOF_STOP) done = true;
Clean up:
// Clean up.
delete solver;
delete matrix;
delete rhs;
delete adaptivity;
for(int i = 0; i < ref_spaces->size(); i++)
delete (*ref_spaces)[i]->get_mesh();
delete ref_spaces;
delete dp;
// Increase counter.
as++;
Now we can show some numerical results. First let us show the resulting meshes for u and v obtained using conventional (single-mesh) hp-FEM: 9,330 DOF (4665 for each solution component).
Next we show the resulting meshes for u and v obtained using the multimesh hp-FEM: 1,723 DOF (49 DOF for u and 1,673 for v).
Finally let us compare the DOF and CPU convergence graphs for both cases:
DOF convergence graphs:
CPU time convergence graphs: