Based on the answers on my last question (Using sacado in example 15) I
thougth that I understood enough in order to extend it to the next step,
using Sacado for solving the time-dependent heat-equation. Thus, the
function F() is
F(u_{n+1}, u_n) = u_{n+1}-u_n - dt*theta*nabla^2 u_{n+1} -
dt*(1-theta)*nabla^2u_n
and when integrating with a test function, I get
(F, z)=(u_{n+1}, z)-(u_n, z) + dt*theta*(nabla u_{n+1}, nabla z) +
dt*(1-theta)*(nabla u_n, nabla z)
When putting that in code I obtain
template <int dim>
void HeatEquation<dim>::assemble_system(const double &local_time)
{
const QGauss<dim> quadrature_formula(fe.degree+1);
system_matrix = 0;
system_rhs = 0;
FEValues<dim> fe_values(fe, quadrature_formula, update_values |
update_q_points | 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();
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell
);
std::vector<double> residual_derivatives(dofs_per_cell);
for(auto cell = dof_handler.begin_active(); cell != dof_handler.end
(); ++cell)
{
fe_values.reinit(cell);
cell->get_dof_indices(local_dof_indices);
Table<2, Sacado::Fad::DFad<double>> value_u(n_q_points, dim),
value_u_old(n_q_points, dim);
Table<2, Sacado::Fad::DFad<double>> gradient_u(n_q_points, dim),
gradient_u_old(n_q_points, dim);
std::vector<Sacado::Fad::DFad<double>> ind_local_dof_values(
dofs_per_cell);
for(size_t i = 0; i < dofs_per_cell; ++i)
{
ind_local_dof_values[i] = present_solution(local_dof_indices
[i]);
ind_local_dof_values[i].diff(i, dofs_per_cell);
}
for(size_t i = 0; i < n_q_points; ++i)
for(size_t j = 0; j < dim; ++j)
{
value_u[i][j] = 0;
gradient_u[i][j] = 0;
value_u_old[i][j] = 0;
gradient_u_old[i][j] = 0;
}
for(size_t q = 0; q < n_q_points; ++q)
{
for(size_t i = 0; i < dofs_per_cell; ++i)
for(size_t d = 0; d < dim; ++d)
{
value_u[q][d] += ind_local_dof_values[i] * fe_values
.shape_value(i, q);
value_u_old[q][d] += old_solution(local_dof_indices[
i]) * fe_values.shape_value(i, q);
gradient_u[q][d] += ind_local_dof_values[i] *
fe_values.shape_grad(i, q)[d] * time_step * theta;
gradient_u_old[q][d] += old_solution(
local_dof_indices[i]) * fe_values.shape_grad(i, q)[d] * time_step * (1 -
theta);
}
}
for(size_t i = 0; i < dofs_per_cell; ++i)
{
Sacado::Fad::DFad<double> R_i = 0;
for(size_t q = 0; q < n_q_points; ++q)
for(size_t d = 0; d < dim; ++d)
R_i += (value_u[q][d] * fe_values.shape_value(i, q)
- value_u_old[q][d] * fe_values.shape_value(i, q) + gradient_u[q][d] *
fe_values.shape_grad(i, q)[d] + gradient_u_old[q][d] * fe_values.shape_grad(
i, q)[d])*fe_values.JxW(q);
for(size_t j = 0; j < dofs_per_cell; ++j)
residual_derivatives[j] = R_i.fastAccessDx(j);
system_matrix.add(local_dof_indices[i], local_dof_indices,
residual_derivatives);
system_rhs(local_dof_indices[i]) -= R_i.val();
}
}
constraint_matrix.condense(system_matrix);
constraint_matrix.condense(system_rhs);
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler, 0,
ZeroFunction<dim>(), boundary_values);
MatrixTools::apply_boundary_values(boundary_values, system_matrix,
newton_update, system_rhs);
}
and for the residual I have
R = u_{n+1}-u_n-theta*dt*\nabla^2u_{n+1}-(1-theta)*dt*nabla^2u_n
as in code
template <int dim>
double HeatEquation<dim>::compute_residual(const double alpha) const
{
Vector<double> residual(dof_handler.n_dofs());
Vector<double> evaluation_point(dof_handler.n_dofs());
evaluation_point = present_solution;
evaluation_point.add(alpha, newton_update);
const QGauss<dim> quadrature_formula(fe.degree+1);
FEValues<dim> fe_values(fe, quadrature_formula, 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();
Vector<double> cell_rhs(dofs_per_cell);
std::vector<Tensor<1, dim> > gradients(n_q_points);
std::vector<double> values(n_q_points);
std::vector<Tensor<1, dim> > gradients_old(n_q_points);
std::vector<double> values_old(n_q_points);
std::vector<types::global_dof_index> local_dof_indices (
dofs_per_cell);
for(auto cell = dof_handler.begin_active(); cell != dof_handler.end
(); ++cell)
{
cell_rhs = 0;
fe_values.reinit(cell);
fe_values.get_function_gradients(evaluation_point, gradients);
fe_values.get_function_values(evaluation_point, values);
fe_values.get_function_gradients(old_solution, gradients_old);
fe_values.get_function_values(old_solution, values_old);
for(size_t q_point = 0; q_point < n_q_points; ++q_point)
{
for(size_t i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) -= (values[q_point] - values_old[q_point] -
fe_values.shape_grad(i, q_point) * time_step * theta * gradients[q_point] -
fe_values.shape_grad(i, q_point) * time_step * (1 - theta) * gradients_old[
q_point]) * fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
for(size_t i = 0; i < dofs_per_cell; ++i)
residual(local_dof_indices[i]) += cell_rhs(i);
}
constraint_matrix.condense(residual);
std::vector<bool> boundary_dofs(dof_handler.n_dofs());
DoFTools::extract_boundary_dofs(dof_handler, ComponentMask(),
boundary_dofs);
for(size_t i = 0; i < dof_handler.n_dofs(); ++i)
if(boundary_dofs[i])
residual(i) = 0;
return residual.l2_norm();
}
Nevertheless, my solver fails in the second iteration with -nan as
residual. What would be the best approach to debug now? Is my code wrong,
or already the equation?
Thanks!
--
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+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
