Thanks for giving feedback as to what the source of the issue was. I'm glad
that you resolved your problem!
On Tuesday, August 29, 2017 at 10:09:57 AM UTC+2, Maxi Miller wrote:
>
> Found my problem, it was the remeshing. Whenever I was remeshing, I
> basically killed the old solution, which then resulted in nan-values for
> it. That in turn generated -nan-values for the system matrix.
>
> Am Montag, 28. August 2017 09:52:45 UTC+2 schrieb Uwe Köcher:
>>
>> looks like your system matrix might be wrong. I think this line
>> system_matrix.add(local_dof_indices[i], local_dof_indices,
>> residual_derivatives);
>> has a missing [i] or [j] or similar on the second local_dof_indices.
>>
>> Getting nan is typical if you have a wrong matrix or right hand side
>> (e.g. forgetting boundary conditions).
>>
>> On Friday, 25 August 2017 17:01:02 UTC+2, Maxi Miller wrote:
>>>
>>> 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.
