Incorrect results in solving the heat equation by DG in time and there is no obvious reason for that?
Mohammad Sabawi
Dear Deal.II Developers and Users,
I am using Deal.II library for solving the linear heat equation by using cG(1)dG(1)(standard Galerkin finite element in space and discontinuous Galerkin method in time) :
u_t − D∆u=f,
u = 0,x∈∂Ω, ( 1)
u(0) = u0 ,
where ∆ is the laplacian operator, D is the diffusion coefficient and Ω ⊂ Rd, where d = 1,2,3, f is righthand side function and u0 is the initial condition function.
I have followed the same approach used in steps 20, 21, and 22 for solving the block systems arising from the discretisation of the mixed Laplace equation in solving my problem, since the solution function is discontinuous in time at each time node then on each time interval I have two unknowns vectors to be determined the first one represents the righthand side limit U+,n−1 and the other one represents the lefthand side limit U-,n.
I have 2 × 2 block matrix system A, 2 x 1 block righthand side vector F and 2 x 1 block solution vector S. I assembled the block system matrix as one local matrix and then added this local matrix to the global block system matrix but the problem is the solution is not convergent in the L2 and L infinity norms I have a solution but the results seem to be erroneous and I am confused where is the error. I think the source of the error maybe resulted from the setup_system() and assemble_system() because I modified step 20 for solving this problem and the difference is in step 20 we have two solution functions the pressure p which is scalar function and the velocity u which a vector valued function while in my example I have two scalar functions U+,n−1 = U1 and U-,n = U2. I used CG solver and Schur complement in solving the resulted block system. I am confused and I do not exactly what is the source of incorrectness of the results and I am waiting for your help.
Note: The comments inside the code will be written in purple colour, and the important parts with red and blue colours. I have included the constructor, setup_system (), assemble_system () and compute_errors() and I attached the block system which I want to be solved.
The code with comments inside the code:
The Constructor
template<int dim>
HeatEquation<dim>::HeatEquation ():
fe(FE_Q<dim>(1), 1,
FE_Q<dim>(1), 1),
dof_handler(triangulation),
fe_nodes(FE_Q<dim>(1)),
dof_handler_nodes(triangulation),
time_step(1./100)
{}
The Setup System Function
template<int dim>
void HeatEquation<dim>::setup_system()
{
GridGenerator::hyper_cube (triangulation, -1,1);
triangulation.refine_global (4);
dof_handler.distribute_dofs(fe);
dof_handler_nodes.distribute_dofs(fe_nodes);
DoFRenumbering::component_wise (dof_handler);
std::vector<types::global_dof_index> dofs_per_component (2);
DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
const unsigned int n_U1 = dofs_per_component[0],
n_U2 = dofs_per_component[1];
std::cout << "\nNumber of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< "Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< " (" << n_U1 << '+' << n_U2 << ")\n"
<< std::endl;
const unsigned int
n_couplings = dof_handler.max_couplings_between_dofs();
sparsity_pattern.reinit (2,2);
sparsity_pattern.block(0,0).reinit (n_U1, n_U1, n_couplings) ;
sparsity_pattern.block(1,0).reinit (n_U2, n_U1, n_couplings) ;
sparsity_pattern.block(0,1).reinit (n_U1, n_U2, n_couplings) ;
sparsity_pattern.block(1,1).reinit (n_U2, n_U2, n_couplings) ;
sparsity_pattern .collect_sizes ();
DoFTools::make_sparsity_pattern(dof_handler, sparsity_pattern);
sparsity_pattern.compress();
system_matrix.reinit (sparsity_pattern);
solution.reinit (2);
solution.block(0).reinit (n_U1);
solution.block(1).reinit (n_U2);
solution.collect_sizes ();
old_solution.reinit (2);
old_solution.block(0).reinit (n_U1);
old_solution.block(1).reinit (n_U2);
old_solution.collect_sizes ();
system_rhs.reinit (2);
system_rhs.block(0).reinit (n_U1);
system_rhs.block(1).reinit (n_U2);
system_rhs.collect_sizes ();
}
The Assemble System Function
template<int dim>
void HeatEquation<dim>::assemble_system()
{
I initialised here the block system matrix and the block righthand side.
system_matrix = 0;
system_rhs = 0;
QGauss<dim> quadrature_formula(fe.degree+1);
QGauss<1> quadrature_time(fe.degree+1);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const unsigned int n_time_q_points = quadrature_time.size();
FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
RightHandSide<dim> right_hand_side;
right_hand_side.set_time(time);
const BoundaryValues<dim> boundary_values_function;
std::vector<double> rhs_values (n_q_points);
std::vector<double> boundary_values (n_q_points);
std::vector<double> old_solution_values(n_q_points);
const FEValuesExtractors::Scalar U1 (0);
const FEValuesExtractors::Scalar U2 (1);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
I initialised here the local matrix and the local righthand side which I have declared above.
local_matrix = 0;
local_rhs = 0;
fe_values[U2].get_function_values (old_solution, old_solution_values);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const double phi_i_U1 = fe_values[U1].value (i, q_point);
const double phi_i_U2 = fe_values[U2].value (i, q_point);
const Tensor<1,dim> grad_phi_i_U1 = fe_values[U1].gradient (i, q_point);
const Tensor<1,dim> grad_phi_i_U2 = fe_values[U2].gradient (i, q_point);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const double phi_j_U1 = fe_values[U1].value (j, q_point);
const double phi_j_U2 = fe_values[U2].value (j, q_point);
const Tensor<1,dim> grad_phi_j_U1 = fe_values[U1].gradient (j, q_point);
const Tensor<1,dim> grad_phi_j_U2 = fe_values[U2].gradient (j, q_point);
I assembled here the local matrix.
local_matrix(i,j) += ((phi_i_U1 * phi_j_U1 +
time_step * grad_phi_i_U1 * grad_phi_j_U1) +
(0.5 * time_step * grad_phi_i_U2 * grad_phi_j_U1) +
(phi_i_U1 * phi_j_U2 +
0.5 * time_step * grad_phi_i_U1 * grad_phi_j_U2) +
(0.5 * phi_i_U2 * phi_j_U2 +
(1./3) * time_step * grad_phi_i_U2 * grad_phi_j_U2)) *
fe_values.JxW (q_point);
for (unsigned int time_q_point=0; time_q_point<n_time_q_points; ++time_q_point)
{
const double old_s = old_solution_values[q_point];
double quadrature_point_1d = quadrature_time.point(time_q_point)[0];
\\ p1 is the quadrature point in the time interval [t_(n-1), t_n].
double p1 = quadrature_point_1d * time_step + time;
double quadrature_weight_1d = quadrature_time.weight(time_q_point);
right_hand_side.set_time(p1);
In this part here I assembled the local righthand side vector which I have declared above.
local_rhs(i) += ( (phi_i_U1 * phi_i_U2 * old_s)
+ (phi_i_U1
* right_hand_side.value(fe_values.quadrature_point(q_point))
* fe_values.JxW (q_point) * time_step * quadrature_weight_1d))
+ ((1./time_step) * phi_i_U2 * (p1 - time)
* right_hand_side.value(fe_values.quadrature_point(q_point))
* fe_values.JxW (q_point) * time_step * quadrature_weight_1d);
}
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j
I added here the local matrix to the global block system matrix.
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
local_matrix(i,j));
for (unsigned int i=0; i<dofs_per_cell; ++i)
I added here the local righthand side vector to the global block righthand side vector.
system_rhs(local_dof_indices[i]) += local_rhs(i);
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
boundary_values_function,
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
}
The Compute Errors Function
template <int dim>
void HeatEquation<dim>::compute_errors() const
{
const ComponentSelectFunction<dim>
U1_mask (dim-1, dim);
const ComponentSelectFunction<dim>
U2_mask (dim-1, dim);
ExactSolution<dim> exact_solution;
Vector<double> cellwise_errors (triangulation.n_active_cells());
QTrapez<1> q_trapez;
QIterated<dim> quadrature (q_trapez, fe.degree+1);
VectorTools::integrate_difference (dof_handler, solution, exact_solution,
cellwise_errors, quadrature,
VectorTools::L2_norm,
&U1_mask);
const double U1_l2_error = cellwise_errors.l2_norm();
VectorTools::integrate_difference (dof_handler, solution, exact_solution,
cellwise_errors, quadrature,
VectorTools::L2_norm,
&U2_mask);
const double U2_l2_error = cellwise_errors.l2_norm();
std::cout << "Errors: ||e_U1||_L2 = " << U1_l2_error
<< ", ||e_U2||_L2 = " << U2_l2_error << "."
<< std::endl;
VectorTools::integrate_difference (dof_handler, solution, exact_solution,
cellwise_errors, quadrature,
VectorTools::H1_seminorm,
&U1_mask);
const double U1_H1_error = cellwise_errors.l2_norm();
VectorTools::integrate_difference (dof_handler, solution, exact_solution,
cellwise_errors, quadrature,
VectorTools::H1_seminorm,
&U2_mask);
const double U2_H1_error = cellwise_errors.l2_norm();
std::cout << "Errors: ||e_U1||_H1_seminorm = " << U1_l2_error
<< ", ||e_U2||_H1_seminorm = " << U2_l2_error << "."
<< std::endl;
VectorTools::integrate_difference (dof_handler, solution, exact_solution,
cellwise_errors, quadrature,
VectorTools::Linfty_norm,
&U1_mask);
const double U1_Linfty_error = cellwise_errors.linfty_norm();
VectorTools::integrate_difference (dof_handler, solution, exact_solution,
cellwise_errors, quadrature,
VectorTools::Linfty_norm,
&U2_mask);
const double U2_Linfty_error = cellwise_errors.linfty_norm();
std::cout << "Errors: ||e_U1||_L_infinity = " << U1_Linfty_error
<< ", ||e_U2||_L_infinity = " << U2_Linfty_error << "."
<< std::endl;
}
Wolfgang Bangerth
Mohammad,
> I have 2 × 2 block matrix system A, 2 x 1 block righthand side vector F and 2
> x 1 block solution vector S. I assembled the block system matrix as one local
> matrix and then added this local matrix to the global block system matrix but
> the problem is the solution is not convergent in the L2 and L infinity norms I
> have a solution but the results seem to be erroneous and I am confused where
> is the error.
first have to work on this some more yourself and see whether you can find out
where the problem is. For example, you don't say what you mean by erroneous:
(i) does the solution *look* correctly? (ii) if you simplify the problem by
not using a right hand side, does it work? (iii) if you simplify the problem
by setting the initial conditions to zero but use a nonzero right hand side,
does it work? (iv) if you run the program in 1d, does it work? Using a
sequence of questions such as this, you can often already narrow down where
the problem may lie.
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/
Guido Kanschat
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
--- You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
--
Interdisziplinäres Zentrum für Wissenschaftliches Rechnen
Universität Heidelberg
Im Neuenheimer Feld 368, 69120 Heidelberg
Mohammad Sabawi
Dear Bangerth and Guido,
...