Dear Dealii group,
I am currently learning Dealii (and i really like this lib). I would
like to bench it in order to compare its performance with an in house
FEM code at my lab.
I want to solve the stokes equation with a static condensation of the
pressure over the velocity dof. I started from the step-22.cc example.
It's working fine. I also computed the "pressure", to within one
constant. I put the code in attachment.
However, i don't know how to compute the residual ! It should be easy,
but i still haven't succeeded. If I extract the different submatrices of
the system matrix at the element level, i don't have the right type of
{U} to compute {P}:
{P} = -[Kpp]^{-1}*[Kpu]*{U}
The matrices are in FullMatrix format, and {u} is of type
std::vector<Tensor<1, dim>> if I'm using "get_function_values". Thus i
can't compute [Kpu]*{U}. I may have the nodal velocity values, but they
have to be ordered as [Kup], [Kpu], ...
Can you give me some advice on the correct way to solve this problem?
Thanks,
Yann
PS :
1) {U} : column vector U
2) I'm using a penalisation method for the pressure, by defining div u =
\lambda p, \lambda the penalisation parameter. Thus the weak form looks
like (v shape function of u, and phi shape function of p),
2*a(grad v, grad u)-a(div v,p)-a(phi,div u)+a(phi,\lambda p) = l(v*f)
then
{P}=-[Kpp]^{-1}[Kpu]*{U}
and finally
([Kuu]-[Kup][Kpp]^-1[Kpu]){U}={f}
--
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 on the web visit
https://groups.google.com/d/msgid/dealii/ab97be61-3dff-4781-a17b-4ceba41a87eb%40gmail.com.
/* ---------------------------------------------------------------------
*
* Test de la solution exacte pour Stokes.
* Attention : raffine sur l'erreur calcul� sur une energie.
*
* ---------------------------------------------------------------------
*/
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_ilu.h>
// Finally, include the header file that declares the Timer class that we will
// use to determine how much time each of the operations of our program takes:
#include <deal.II/base/timer.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <memory>
namespace Step22 {
using namespace dealii;
template<int dim>
class StokesProblem {
public:
StokesProblem(const unsigned int degree);
void run();
private:
void make_grid();
void setup_dofs();
void assemble_system();
void computeP();
void determine_component_extractors();
void solve_sc(bool computePress);
void checkSolution();
void output_results(const unsigned int refinement_cycle) const;
void refine_mesh(bool refineOnError);
const unsigned int degree;
const double lambda = 1e-8;
enum {
u_dof = 0,
p_dof = 1
};
std::vector<types::global_dof_index> element_indices_u;
std::vector<types::global_dof_index> element_indices_p;
unsigned int nbDofPress;
Triangulation<dim> triangulation;
FESystem<dim> velocity_fe;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
AffineConstraints<double> zero_constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockVector<double> solution;
BlockVector<double> system_rhs;
TimerOutput computing_timer;
};
template<int dim>
StokesProblem<dim>::StokesProblem(const unsigned int degree)
: degree(degree),
triangulation(Triangulation<dim>::maximum_smoothing),
velocity_fe(FE_Q<dim>(degree), dim),
fe(velocity_fe, 1, FE_Q<dim>(degree - 1), 1),
dof_handler(triangulation),
computing_timer(std::cout, TimerOutput::summary,
TimerOutput::wall_times) {}
template<int dim>
class RightHandSide : public Function<dim> {
public:
RightHandSide()
: Function<dim>(dim + 1) {}
virtual void vector_value(const Point<dim> &p,
Vector<double> &value) const override;
};
template<int dim>
void RightHandSide<dim>::vector_value(const Point<dim> &p,
Vector<double> &values) const {
const double Px = p[0];
const double Py = p[1];
const double mu = 1.;
constexpr double pi = numbers::PI;
constexpr double pi2 = numbers::PI * numbers::PI;
constexpr double pi3 = numbers::PI * numbers::PI * numbers::PI;
values[0] = 2 * pi * std::cos(2 * pi * Px) + mu * pi2 * (std::sin(pi * Px)
+ std::sin(pi * Py));
values[1] = 2 * pi * std::cos(2 * pi * Py) - mu * pi3 * std::cos(pi * Px) *
Py;
values[2] = 0;
}
template<int dim>
class ExactSolution : public Function<dim> {
public:
ExactSolution()
: Function<dim>(dim + 1) {}
virtual void vector_value(const Point<dim> &p,
Vector<double> &values) const override;
};
template<int dim>
void ExactSolution<dim>::vector_value(const Point<dim> &p,
Vector<double> &values) const {
const double Px = p[0];
const double Py = p[1];
constexpr double pi = numbers::PI;
values[0] = std::sin(pi * Px) + std::sin(pi * Py);
values[1] = -pi * std::cos(pi * Px) * Py;
values[2] = std::sin(2 * pi * Px) + std::sin(2 * pi * Py);
}
template<int dim>
void StokesProblem<dim>::setup_dofs() {
system_matrix.clear();
dof_handler.distribute_dofs(fe);
DoFRenumbering::Cuthill_McKee(dof_handler);
std::vector<unsigned int> block_component(2);
block_component[0] = 0;
block_component[1] = 1;
DoFRenumbering::block_wise(dof_handler);
const FEValuesExtractors::Vector velocities(0);
{
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
ExactSolution<dim>(),
constraints,
fe.component_mask(velocities));
}
constraints.close();
{
zero_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<dim>(
dim + 1),
zero_constraints,
fe.component_mask(velocities));
}
zero_constraints.close();
const std::vector<types::global_dof_index> dofs_per_block =
DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
const types::global_dof_index n_u = dofs_per_block[0];
const types::global_dof_index n_p = dofs_per_block[1];
nbDofPress = n_p;
std::cout << " Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << n_u << '+' << n_p << ')' << std::endl;
{
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from(dsp);
}
system_matrix.reinit(sparsity_pattern);
solution.reinit(dofs_per_block);
system_rhs.reinit(dofs_per_block);
}
template<int dim>
void StokesProblem<dim>::determine_component_extractors() {
element_indices_u.clear();
element_indices_p.clear();
for (unsigned int k = 0; k < fe.n_dofs_per_cell(); ++k) {
const unsigned int k_group = fe.system_to_base_index(k).first.first;
if (k_group == u_dof)
element_indices_u.push_back(k);
else if (k_group == p_dof)
element_indices_p.push_back(k);
else {
Assert(k_group <= p_dof, ExcInternalError());
}
}
}
template<int dim>
void StokesProblem<dim>::assemble_system() {
int debug = 0;
system_matrix = 0;
system_rhs = 0;
QGauss<dim> quadrature_formula(degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_quadrature_points |
update_JxW_values | update_gradients);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
std::vector<SymmetricTensor<2, dim>> symgrad_phi_u(dofs_per_cell);
std::vector<double> div_phi_u(dofs_per_cell);
std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
std::vector<double> phi_p(dofs_per_cell);
// Pour la condensation statique
const unsigned int n_u = element_indices_u.size();
const unsigned int n_p = element_indices_p.size();
FullMatrix<double> k_orig(dofs_per_cell, dofs_per_cell);
FullMatrix<double> k_uu(n_u, n_u);
FullMatrix<double> k_up(n_u, n_p);
FullMatrix<double> k_pu(n_p, n_u);
FullMatrix<double> k_pp(n_p, n_p);
FullMatrix<double> k_pp_inv(n_p, n_p);
FullMatrix<double> A(n_p, n_u);
FullMatrix<double> B(n_u, n_u);
FullMatrix<double> C(n_p, n_p);
FullMatrix<double> local_matrix_sc(dofs_per_cell, dofs_per_cell);
for (const auto &cell: dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
local_matrix = 0;
local_rhs = 0;
right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
rhs_values);
for (unsigned int q = 0; q < n_q_points; ++q) {
for (unsigned int k = 0; k < dofs_per_cell; ++k) {
symgrad_phi_u[k] =
fe_values[velocities].symmetric_gradient(k, q);
div_phi_u[k] = fe_values[velocities].divergence(k, q);
phi_u[k] = fe_values[velocities].value(k, q);
phi_p[k] = fe_values[pressure].value(k, q);
}
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
for (unsigned int j = 0; j <= i; ++j) {
local_matrix(i, j) +=
(2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) // (1)
- div_phi_u[i] * phi_p[j] // (2)
- phi_p[i] * div_phi_u[j] // (3)
+ lambda * phi_p[i] * phi_p[j]
) * fe_values.JxW(q); // * dx
}
const unsigned int component_i =
fe.system_to_component_index(i).first;
local_rhs(i) += fe_values.shape_value(i, q) // phi_u_i(x_q)
* rhs_values[q](component_i) // * f(x_q)
* fe_values.JxW(q); // * dx
}
}
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = i + 1; j < dofs_per_cell; ++j) {
local_matrix(i, j) = local_matrix(j, i);
}
//Static condentation of the P into U, block(0,0)
{
k_uu = 0;
k_up = 0;
k_pu = 0;
k_pp = 0;
k_pp_inv = 0;
A = 0;
B = 0;
local_matrix_sc = 0;
k_uu.extract_submatrix_from(local_matrix,
element_indices_u,
element_indices_u);
k_pu.extract_submatrix_from(local_matrix,
element_indices_p,
element_indices_u);
k_up.extract_submatrix_from(local_matrix,
element_indices_u,
element_indices_p);
k_pp.extract_submatrix_from(local_matrix,
element_indices_p,
element_indices_p);
if (debug) {
std::ofstream out0("local_matrix.txt");
local_matrix.print_formatted(out0, 3, true, 0, "0", 1., 1e-12);
std::ofstream out1("k_uu.txt");
k_uu.print_formatted(out1, 3, true, 0, "0", 1., 1e-12);
std::ofstream out2("k_pu.txt");
k_pu.print_formatted(out2, 3, true, 0, "0", 1., 1e-12);
std::ofstream out3("k_up.txt");
k_up.print_formatted(out3, 3, true, 0, "0", 1., 1e-12);
std::ofstream out4("k_pp.txt");
k_pp.print_formatted(out4, 3, true, 0, "0", 1., 1e-12);
}
/*
/ \
| Kuu Kup | [Kpu]U+[Kpp]P=0 => P=-[Kpp]^-1[Kpu]U
| Kpu Kpp | [Kuu]U+[Kup]P=f => ([Kuu]-[Kup][Kpp]^-1[Kpu])U=f
\ /
*/
k_pp_inv.invert(k_pp);
k_pp_inv.mmult(A, k_pu); // A = [K_pp]^-1*[K_pu]
k_up.mmult(B, A); // B = [K_up][A]
B.scatter_matrix_to(element_indices_u,
element_indices_u,
local_matrix_sc);
local_matrix.add(-1.0, local_matrix_sc);
if (debug) {
std::ofstream out5("B.txt");
B.print_formatted(out5, 3, true, 0, "0", 1., 1e-12);
std::ofstream out6("local_matrix_sc.txt");
local_matrix.print_formatted(out6, 3, true, 0, "0", 1., 1e-12);
std::ofstream out7("A.txt");
A.print_formatted(out7, 3, true, 0, "0", 1., 1e-12);
std::ofstream out8("k_pp_inv.txt");
k_pp_inv.print_formatted(out8, 3, true, 0, "0", 1., 1e-12);
C = 0;
k_pp_inv.mmult(C, k_pp); // A = [K_pp]^-1*[K_pu]
std::ofstream out9("Identite.txt");
C.print_formatted(out9, 3, true, 0, "0", 1., 1e-12);
}
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(local_matrix,
local_rhs,
local_dof_indices,
system_matrix,
system_rhs);
}
}
template<int dim>
void StokesProblem<dim>::computeP() {
const QGauss<dim> quadrature_formula(degree + 1);
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_quadrature_points |
update_JxW_values | update_gradients);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
const double unSurLambda = 1. / lambda;
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
Vector<double> local_pression(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
std::vector<Tensor<2, dim>> present_velocity_gradients(n_q_points);
double phi_p = 0.;
double present_velocity_divergence = 0.;
for (const auto &cell: dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
local_pression = 0;
fe_values[velocities].get_function_gradients(
solution, present_velocity_gradients);
for (unsigned int q = 0; q < n_q_points; ++q) {
present_velocity_divergence = trace(present_velocity_gradients[q]);
for (unsigned int i = 0; i < dofs_per_cell; ++i) {
phi_p = fe_values[pressure].value(i, q);
local_pression(i) += unSurLambda * phi_p *
present_velocity_divergence *
fe_values.JxW(q);
}
}
cell->get_dof_indices(local_dof_indices);
zero_constraints.distribute_local_to_global(local_pression,
local_dof_indices, solution);
}
}
template<int dim>
void StokesProblem<dim>::solve_sc(bool computePress) {
std::cout << "Solving linear system... ";
Timer timer;
SparseDirectUMFPACK A_direct;
solution.block(0) = system_rhs.block(0);
//solution.block(1) = 0.;
A_direct.solve(system_matrix.block(0, 0), solution.block(0));
constraints.distribute(solution);
// On recup le champs de pression
if (computePress) {
//solution.block(1) = 0;
constraints.distribute(solution);
computeP();
}
constraints.distribute(solution);
timer.stop();
std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
}
template<int dim>
void
StokesProblem<dim>::checkSolution() {
{
const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
dim + 1);
Vector<double> cellwise_errors(triangulation.n_active_cells());
QGauss<dim> quadrature(degree + 2);
VectorTools::integrate_difference(dof_handler,
solution,
ExactSolution<dim>(),
cellwise_errors,
quadrature,
VectorTools::L2_norm,
&velocity_mask);
const double error_u_l2 =
VectorTools::compute_global_error(triangulation,
cellwise_errors,
VectorTools::L2_norm);
VectorTools::integrate_difference(dof_handler,
solution,
ExactSolution<dim>(),
cellwise_errors,
quadrature,
VectorTools::L2_norm,
&pressure_mask);
const double error_p_l2 =
VectorTools::compute_global_error(triangulation,
cellwise_errors,
VectorTools::L2_norm);
std::cout << "error: u_0: " << error_u_l2 << " p_0: " << error_p_l2
<< std::endl;
}
}
template<int dim>
void
StokesProblem<dim>::output_results(const unsigned int refinement_cycle) const
{
std::vector<std::string> solution_names(dim, "velocity");
solution_names.emplace_back("pressure");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation(
dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation.push_back(
DataComponentInterpretation::component_is_scalar);
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution,
solution_names,
DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.build_patches();
std::ofstream output(
"solution-" + Utilities::int_to_string(refinement_cycle, 2) +
".vtk");
data_out.write_vtk(output);
}
template<int dim>
void StokesProblem<dim>::refine_mesh(bool refineOnError) {
if (refineOnError) {
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
const FEValuesExtractors::Vector velocity(0);
KellyErrorEstimator<dim>::estimate(
dof_handler,
QGauss<dim - 1>(degree + 1),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell,
fe.component_mask(velocity));
GridRefinement::refine_and_coarsen_fixed_number(triangulation,
estimated_error_per_cell,
0.3,
0.0);
triangulation.execute_coarsening_and_refinement();
} else {
triangulation.refine_global(1);
}
}
template<int dim>
void StokesProblem<dim>::make_grid() {
GridGenerator::hyper_cube(triangulation, 0, 1);
triangulation.refine_global(3);
}
template<int dim>
void StokesProblem<dim>::run() {
double current_res = 1.0;
make_grid();
for (unsigned int refinement_cycle = 0; refinement_cycle < 4;
++refinement_cycle) {
std::cout << "Refinement cycle " << refinement_cycle << std::endl;
if (refinement_cycle > 0)
refine_mesh(false);
setup_dofs();
determine_component_extractors();
std::cout << " Assembling..." << std::endl << std::flush;
assemble_system();
std::cout << " Solving..." << std::flush;
solve_sc(true);
checkSolution();
compute_residual();
system_rhs.block(1)=0;
current_res = system_rhs.l2_norm();
std::cout << "The residual is " << current_res
<< std::endl;
output_results(refinement_cycle);
std::cout << std::endl;
}
}
} // namespace Step22
int main()
{
try
{
using namespace Step22;
StokesProblem<2> flow_problem(2);
flow_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;
}
