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.