I would like to solve the simple heat equation u_t = \nabla^2 u (and to
compare, in addition the equation -\nabla^2 u = c). For that I am using the
following assemble_matrix-function:
template <int dim> void MinimalSurfaceProblem<dim>::assemble_system ()
{
const QGauss<dim> quadrature_formula(fe.degree+1);
const FEValuesExtractors::Vector surface_N(0);
const FEValuesExtractors::Vector surface_TE(dim);
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();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
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)
{
if(cell->is_locally_owned())
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);
cell->get_dof_indices(local_dof_indices);
Table<3, Sacado::Fad::DFad<double> > grad_N (n_q_points, dim, dim),
grad_N_old(n_q_points, dim, dim);
Table<3, Sacado::Fad::DFad<double> > grad_TE (n_q_points, dim, dim),
grad_TE_old(n_q_points, dim, dim);
Table<2, Sacado::Fad::DFad<double> > val_N (n_q_points, dim), val_N_old(
n_q_points, dim);
Table<2, Sacado::Fad::DFad<double> > val_TE (n_q_points, dim), val_TE_old(
n_q_points, dim);
std::vector<Sacado::Fad::DFad<double> > ind_local_dof_values(dofs_per_cell
);
std::vector<double> local_dof_values(cell->get_fe().dofs_per_cell);
std::vector<Tensor<1, dim, Sacado::Fad::DFad<double>>> grad_N_AD(n_q_points
);
cell->get_dof_values(present_solution, local_dof_values.begin(),
local_dof_values.end());
std::vector<Sacado::Fad::DFad<double>> local_dof_values_AD(cell->get_fe().
dofs_per_cell);
for (unsigned int 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);
local_dof_values_AD[i] = present_solution(local_dof_indices[i]);
local_dof_values_AD[i].diff(i, dofs_per_cell);
}
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int d=0; d<dim; ++d)
{
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
grad_N[q][d][local_dim] = 0;
grad_N_old[q][d][local_dim] = 0;
grad_TE[q][d][local_dim] = 0;
grad_TE_old[q][d][local_dim] = 0;
}
}
}
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int d=0; d<dim; ++d)
{
val_N[q][d] = 0;
val_N_old[q][d] = 0;
val_TE[q][d] = 0;
val_TE_old[q][d] = 0;
}
}
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int d = 0; d < dim; d++)
{
val_N[q][d] += ind_local_dof_values[i] * fe_values[surface_N].value(i, q)[d
];
val_N_old[q][d] += old_solution(local_dof_indices[i]) * fe_values[surface_N
].value(i, q)[d];
val_TE[q][d] += ind_local_dof_values[i] * fe_values[surface_TE].value(i, q
)[d];
val_TE_old[q][d] += old_solution(local_dof_indices[i]) * fe_values[
surface_TE].value(i, q)[d];
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
grad_N[q][d][local_dim] += ind_local_dof_values[i] * fe_values[surface_N].
gradient(i, q)[d][local_dim];
grad_N_old[q][d][local_dim] += old_solution(local_dof_indices[i]) *
fe_values[surface_N].gradient(i, q)[d][local_dim];
grad_TE[q][d][local_dim] += ind_local_dof_values[i] * fe_values[surface_TE
].gradient(i, q)[d][local_dim];
grad_TE_old[q][d][local_dim] += old_solution(local_dof_indices[i]) *
fe_values[surface_TE].gradient(i, q)[d][local_dim];
}
}
}
}
//fe_values[surface_N].get_function_gradients_from_local_dof_values(local_dof_values_AD,
grad_N_AD);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
Sacado::Fad::DFad<double> R_i = 0;
for(unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int d = 0; d < dim; d++)
{
R_i += (((val_TE[q][d] - val_TE_old[q][d]) / time_step) * fe_values[
surface_TE].value(i, q)[d]/* - 5 * fe_values[surface_TE].value(i, q)[d]*/) *
fe_values.JxW(q);
for(size_t local_dim = 0; local_dim < dim; ++local_dim)
{
R_i += ((theta * grad_TE[q][d][local_dim] + (1 - theta) * grad_TE_old[q][d
][local_dim]) * fe_values[surface_TE].gradient(i, q)[d][local_dim])
* fe_values.JxW(q);
R_i += -(grad_N[q][d][local_dim] * fe_values[surface_N].gradient(i, q)[d][
local_dim] + 5 * fe_values[surface_N].value(i, q)[d]) * fe_values.JxW(q);
}
}
}
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
residual_derivatives[j] = R_i.fastAccessDx(j);
}
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix(i, j) += residual_derivatives[j];
}
cell_rhs(i) -= R_i.val();
}
//cell->get_dof_indices (local_dof_indices);
newton_constraints.distribute_local_to_global(cell_matrix, cell_rhs,
local_dof_indices, system_matrix, system_rhs);
}
}
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
}
Now, when setting theta = 1, i.e. without using the old gradient values, I
receive the correct solution for the first step:
<https://lh3.googleusercontent.com/-du9C_uqA52Y/WecHoXLuPHI/AAAAAAAACUs/X-JG35LSYlUOCbRxGeRfbmRR1ZWvFPdcgCLcBGAs/s1600/visit0000.png>
When setting theta = 0.5, i.e. including the old gradient values, the
result is obviously incorrect:
<https://lh3.googleusercontent.com/-WuRgvrrFJ2U/WecHsqMT4aI/AAAAAAAACUw/g4SmcVztZyYR2EUSx2GItYTwDArA9YZ7wCLcBGAs/s1600/visit0001.png>
Apparently my old values are correct, but my old gradients are not. Thus I
am wondering if i getting the gradient values from the old solution in an
incorrect way, or is it correct?
--
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 [email protected].
For more options, visit https://groups.google.com/d/optout.
