Dear all,
I have a question regarding an *extension* of step.74.
Eventually, I would like to solve a non-linear Poisson equation using DG
elements. I would also like to use automatic differentiation to obtain the
Jacobian from the Residual as in step-72.
To start, I wanted to solve the linear problem as if it was non-linear
using a single Newton iteration.
For the cell_worker and boundary_worker, I was able to use calls like the
following
std::vector<ADNumberType> old_solution_values(fe_fv.n_quadrature_points);
fe_fv[u_fe].get_function_values_from_local_dof_values(dof_values_ad,
old_solution_values);
std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(fe_fv.
n_quadrature_points);
fe_fv[u_fe].get_function_gradients_from_local_dof_values(dof_values_ad,
old_solution_gradients);
to then assemble the residual and automatically compute the Jacobian. The
full code is appended.
For the face_worker, however, I think I would need functions like
std::vector<Tensor<1, dim, ADNumberType>> old_average_of_gradients(fe_fv.
n_quadrature_points);
fe_iv[u_fe].get_average_of_gradients_from_local_dof_values(dof_values_ad,
old_average_of_gradients);
Since I didn't find these functions, I think they might not exist yet.
Therefore my question.
- If these functions don't exist, what alternative approach do you see to
get the desired average of gradients of the old solution vector?
Thank you very much in advance for any help and best regards,
Nik
--
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].
To view this discussion visit
https://groups.google.com/d/msgid/dealii/3ca01bd8-4f57-4cdb-8570-16f09c7d7f0dn%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_lib.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_interface_values.h>
#include <deal.II/numerics/derivative_approximation.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/meshworker/copy_data.h>
#include <deal.II/meshworker/mesh_loop.h>
#include <deal.II/meshworker/scratch_data.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/differentiation/ad.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/fe/fe_q.h>
namespace Step74 {
using namespace dealii;
double get_penalty_factor(const unsigned int fe_degree, const double cell_extent_left,
const double cell_extent_right) {
const unsigned int degree = std::max(1U, fe_degree);
return degree * (degree + 1.) * 0.5 * (1. / cell_extent_left + 1. / cell_extent_right);
}
struct CopyDataFace {
FullMatrix<double> cell_matrix;
Vector<double> cell_rhs;
std::vector<types::global_dof_index> joint_dof_indices;
std::array<double, 2> values;
std::array<unsigned int, 2> cell_indices;
};
struct CopyData {
FullMatrix<double> cell_matrix;
Vector<double> cell_rhs;
std::vector<types::global_dof_index> local_dof_indices;
std::vector<CopyDataFace> face_data;
double value;
unsigned int cell_index;
template <class Iterator>
void reinit(const Iterator &cell, const unsigned int dofs_per_cell) {
cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
cell_rhs.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
}
};
template <int dim>
class SIPGLaplace {
public:
SIPGLaplace();
void run();
private:
void setup_system();
void assemble_system();
void solve();
void refine_grid();
void output_results(const unsigned int cycle) const;
void compute_errors();
void compute_error_estimate();
double compute_energy_norm_error();
Triangulation<dim> triangulation;
const unsigned int degree;
const QGauss<dim> quadrature;
const QGauss<dim - 1> face_quadrature;
const MappingQ1<dim> mapping;
using ScratchData = MeshWorker::ScratchData<dim>;
const FE_DGQ<dim> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution;
Vector<double> current_solution;
Vector<double> system_rhs;
ConvergenceTable convergence_table;
const double diffusion_coefficient = 1.;
};
template <int dim>
SIPGLaplace<dim>::SIPGLaplace()
: degree(1),
quadrature(degree + 1),
face_quadrature(degree + 1),
mapping(),
fe(degree),
dof_handler(triangulation) {}
template <int dim>
void SIPGLaplace<dim>::setup_system() {
dof_handler.distribute_dofs(fe);
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
solution.reinit(dof_handler.n_dofs());
current_solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void SIPGLaplace<dim>::assemble_system() {
using ADHelper =
Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,
double>;
using ADNumberType = typename ADHelper::ad_type;
const FEValuesExtractors::Scalar u_fe(0);
const auto cell_worker = [&u_fe, this](const auto &cell, auto &scratch_data, auto ©_data) {
const FEValues<dim> &fe_v = scratch_data.reinit(cell);
const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
copy_data.reinit(cell, dofs_per_cell);
FullMatrix<double> local_cell_matrix;
Vector<double> local_rhs;
local_cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
local_rhs.reinit(dofs_per_cell);
std::vector<types::global_dof_index> &local_dof_indices = copy_data.local_dof_indices;
cell->get_dof_indices(local_dof_indices);
const unsigned int n_independent_variables = local_dof_indices.size();
const unsigned int n_dependent_variables = dofs_per_cell;
ADHelper ad_helper(n_independent_variables, n_dependent_variables);
ad_helper.register_dof_values(current_solution, local_dof_indices);
const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(fe_v.n_quadrature_points);
fe_v[u_fe].get_function_gradients_from_local_dof_values(dof_values_ad,
old_solution_gradients);
std::vector<ADNumberType> residual_ad(n_dependent_variables, ADNumberType(0.0));
for (const unsigned int q : fe_v.quadrature_point_indices()) {
for (const unsigned int i : fe_v.dof_indices()) {
residual_ad[i] +=
diffusion_coefficient *
(fe_v.shape_grad(i, q) // \nabla \phi_i // * a_n
* old_solution_gradients[q]) // * \nabla u_n
* fe_v.JxW(q); // * dx
}
}
ad_helper.register_residual_vector(residual_ad);
ad_helper.compute_residual(local_rhs);
local_rhs *= -1.0;
ad_helper.compute_linearization(local_cell_matrix);
// const auto &q_points = scratch_data.get_quadrature_points();
// const unsigned int n_q_points = q_points.size();
// const std::vector<double> &JxW = scratch_data.get_JxW_values();
//
// for (unsigned int point = 0; point < n_q_points; ++point)
// for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i) {
// for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
// copy_data.cell_matrix(i, j) += diffusion_coefficient * // nu
// fe_v.shape_grad(i, point) * // grad v_h
// fe_v.shape_grad(j, point) * // grad u_h
// JxW[point]; // dx
//
// copy_data.cell_rhs(i) += 0.0;
// }
// transfer to copy data
for (int i = 0; i < dofs_per_cell; i++) {
for (int j = 0; j < dofs_per_cell; j++) {
copy_data.cell_matrix(i, j) += local_cell_matrix(i, j);
}
copy_data.cell_rhs(i) += local_rhs(i);
}
};
const auto boundary_worker = [&u_fe, this](const auto &cell, const unsigned int &face_no,
auto &scratch_data, auto ©_data) {
const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
FullMatrix<double> local_cell_matrix;
Vector<double> local_rhs;
local_cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
local_rhs.reinit(dofs_per_cell);
std::vector<types::global_dof_index> &local_dof_indices = copy_data.local_dof_indices;
cell->get_dof_indices(local_dof_indices);
const unsigned int n_independent_variables = local_dof_indices.size();
const unsigned int n_dependent_variables = dofs_per_cell;
ADHelper ad_helper(n_independent_variables, n_dependent_variables);
ad_helper.register_dof_values(current_solution, local_dof_indices);
const std::vector<ADNumberType> &dof_values_ad = ad_helper.get_sensitive_dof_values();
std::vector<ADNumberType> old_solution_values(fe_fv.n_quadrature_points);
fe_fv[u_fe].get_function_values_from_local_dof_values(dof_values_ad, old_solution_values);
std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(fe_fv.n_quadrature_points);
fe_fv[u_fe].get_function_gradients_from_local_dof_values(dof_values_ad,
old_solution_gradients);
std::vector<ADNumberType> residual_ad(n_dependent_variables, ADNumberType(0.0));
const auto &q_points = scratch_data.get_quadrature_points();
const unsigned int n_q_points = q_points.size();
// const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
const std::vector<double> &JxW = scratch_data.get_JxW_values();
const std::vector<Tensor<1, dim>> &normals = scratch_data.get_normal_vectors();
const double extent1 = cell->measure() / cell->face(face_no)->measure();
const double penalty = get_penalty_factor(degree, extent1, extent1);
// Dirichlet Condition
if (cell->face(face_no)->boundary_id() == 0) {
for (unsigned int point = 0; point < n_q_points; ++point) {
/* for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
copy_data.cell_matrix(i, j) +=
(-diffusion_coefficient * // - nu
fe_fv.shape_value(i, point) * // v_h
(fe_fv.shape_grad(j, point) * // (grad u_h .
normals[point]) // n)
- diffusion_coefficient * // - nu
(fe_fv.shape_grad(i, point) * // (grad v_h .
normals[point]) * // n)
fe_fv.shape_value(j, point) // u_h
+ diffusion_coefficient * penalty * // + nu sigma
fe_fv.shape_value(i, point) * // v_h
fe_fv.shape_value(j, point) // u_h
) *
JxW[point]; // dx
*/
for (unsigned int i = 0; i < dofs_per_cell; ++i)
residual_ad[i] += (-diffusion_coefficient * // - nu
fe_fv.shape_value(i, point) * // v_h
(old_solution_gradients[point] * // (grad u_h .
normals[point]) // n)
- diffusion_coefficient * // - nu
(fe_fv.shape_grad(i, point) * // (grad v_h .
normals[point]) * // n)
old_solution_values[point] // u_h
+ diffusion_coefficient * penalty * // + nu sigma
fe_fv.shape_value(i, point) * // v_h
old_solution_values[point] // u_h
- diffusion_coefficient * // - nu
(fe_fv.shape_grad(i, point) * // (grad v_h .
normals[point]) * // n)
(cell->face(face_no)->boundary_id()) // g
+ diffusion_coefficient * penalty * // + nu sigma
fe_fv.shape_value(i, point) *
(cell->face(face_no)->boundary_id()) // v_h g
) *
JxW[point]; // dx
}
}
// Neumann Condition
if (cell->face(face_no)->boundary_id() == 1) {
for (unsigned int point = 0; point < n_q_points; ++point) {
// The Matrix remains uncahnged
// Modify the RHS
for (unsigned int i = 0; i < dofs_per_cell; ++i)
residual_ad[i] += (-diffusion_coefficient * // - nu
fe_fv.shape_value(i, point) * // v_h .
(3) // g
) *
JxW[point]; // dx
}
}
ad_helper.register_residual_vector(residual_ad);
ad_helper.compute_residual(local_rhs);
local_rhs *= -1.0;
ad_helper.compute_linearization(local_cell_matrix);
// cell_rhs.print(std::cout);
// copy_data.cell_rhs.print(std::cout);
for (int i = 0; i < dofs_per_cell; i++) {
for (int j = 0; j < dofs_per_cell; j++) {
copy_data.cell_matrix(i, j) += local_cell_matrix(i, j);
}
copy_data.cell_rhs(i) += local_rhs(i);
}
};
const auto face_worker = [&](const auto &cell, const unsigned int &f, const unsigned int &sf,
const auto &ncell, const unsigned int &nf, const unsigned int &nsf,
auto &scratch_data, auto ©_data) {
const FEInterfaceValues<dim> &fe_iv = scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
copy_data.face_data.emplace_back();
CopyDataFace ©_data_face = copy_data.face_data.back();
const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs();
copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
copy_data_face.cell_rhs.reinit(n_dofs_face, n_dofs_face);
FullMatrix<double> local_cell_matrix;
Vector<double> local_rhs;
local_cell_matrix.reinit(n_dofs_face, n_dofs_face);
local_rhs.reinit(n_dofs_face);
const unsigned int n_independent_variables = copy_data_face.joint_dof_indices.size();
const unsigned int n_dependent_variables = copy_data_face.joint_dof_indices.size();
ADHelper ad_helper(n_independent_variables, n_dependent_variables);
const std::vector<double> &JxW = fe_iv.get_JxW_values();
const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
const double extent1 = cell->measure() / cell->face(f)->measure();
const double extent2 = ncell->measure() / ncell->face(nf)->measure();
const double penalty = get_penalty_factor(degree, extent1, extent2);
for (const unsigned int point : fe_iv.quadrature_point_indices()) {
for (const unsigned int i : fe_iv.dof_indices())
for (const unsigned int j : fe_iv.dof_indices())
copy_data_face.cell_matrix(i, j) +=
(-diffusion_coefficient * // - nu
fe_iv.jump_in_shape_values(i, point) * // [v_h]
(fe_iv.average_of_shape_gradients(j,
point) * // ({grad u_h} .
normals[point]) // n)
- diffusion_coefficient * // - nu
(fe_iv.average_of_shape_gradients(i, point) * // (grad v_h .
normals[point]) * // n)
fe_iv.jump_in_shape_values(j, point) // [u_h]
+ diffusion_coefficient * penalty * // + nu sigma
fe_iv.jump_in_shape_values(i, point) * // [v_h]
fe_iv.jump_in_shape_values(j, point) // [u_h]
) *
JxW[point]; // dx
}
};
AffineConstraints<double> constraints;
constraints.close();
const auto copier = [&](const auto &c) {
constraints.distribute_local_to_global(c.cell_matrix, c.cell_rhs, c.local_dof_indices,
system_matrix, system_rhs);
for (auto &cdf : c.face_data) {
constraints.distribute_local_to_global(cdf.cell_matrix, cdf.joint_dof_indices,
system_matrix);
}
};
const UpdateFlags cell_flags =
update_values | update_gradients | update_quadrature_points | update_JxW_values;
const UpdateFlags face_flags = update_values | update_gradients | update_quadrature_points |
update_normal_vectors | update_JxW_values;
ScratchData scratch_data(mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
CopyData copy_data;
MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier,
scratch_data, copy_data,
MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces |
MeshWorker::assemble_own_interior_faces_once,
boundary_worker, face_worker);
}
template <int dim>
void SIPGLaplace<dim>::solve() {
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix);
A_direct.vmult(solution, system_rhs);
}
template <int dim>
void SIPGLaplace<dim>::output_results(const unsigned int cycle) const {
const std::string filename = "sol_Q" + Utilities::int_to_string(degree, 1) + "-" +
Utilities::int_to_string(cycle, 2) + ".vtu";
std::ofstream output(filename);
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
data_out.build_patches(mapping);
data_out.write_vtu(output);
}
template <int dim>
void SIPGLaplace<dim>::run() {
GridGenerator::hyper_cube(triangulation, 0, 1, true);
triangulation.refine_global(4);
{
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell != endc; ++cell) {
if ((cell->center()[1]) > 0.9) {
if ((cell->center()[0] > 0.9) || (cell->center()[0] < 0.1)) cell->set_refine_flag();
}
}
}
{
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
endc = dof_handler.end();
triangulation.execute_coarsening_and_refinement();
for (; cell != endc; ++cell) {
if ((cell->center()[1]) > 0.9) {
if ((cell->center()[0] > 0.9) || (cell->center()[0] < 0.1)) cell->set_refine_flag();
}
}
triangulation.execute_coarsening_and_refinement();
}
std::cout << " Number of active cells : " << triangulation.n_active_cells() << std::endl;
setup_system();
std::cout << " Number of degrees of freedom : " << dof_handler.n_dofs() << std::endl;
assemble_system();
solve();
output_results(1);
std::cout << std::endl;
}
} // namespace Step74
int main(int argc, char **argv) {
using namespace dealii;
try {
Step74::SIPGLaplace<2> problem;
problem.run();
}
catch (std::exception &exc) {
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
return 1;
} catch (...) {
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------" << std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------" << std::endl;
return 1;
}
return 0;
}
