Hi Timo and Wolfgang
I wrapped up my A preconditioner with an inner CG solver (I am using IMEX
scheme which moves the convection term to the RHS so A is symmetric) but
the SparseILU still seems not to work. I attached my code here, if you
could take a look, it is very easy to reproduce the issue:
1. compile and run it, you will see everything works as expected
2. comment out line 221 and uncomment line 222, the program just hangs.
Wolfgang does have a point: if I were able to use a direct solver on A then
I might as well apply it on the whole system. I think I should get
SparseILU working at least.
Thank you again for your helpful comments.
Jie
--
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.
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/timer.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_ilu.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/grid/grid_tools.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_accessor.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/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
// To transfer solutions between meshes, this file is included:
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/physics/elasticity/standard_tensors.h>
#include <fstream>
#include <iostream>
#include <sstream>
namespace fluid
{
using namespace dealii;
void create_triangulation_2d(Triangulation<2> &tria)
{
HyperBallBoundary<2> boundary(Point<2>(0.5, 0.2), 0.05);
Triangulation<2> left, middle, right, tmp, tmp2;
GridGenerator::subdivided_hyper_rectangle(
left,
std::vector<unsigned int>({3U, 4U}),
Point<2>(),
Point<2>(0.3, 0.41),
false);
GridGenerator::subdivided_hyper_rectangle(
right,
std::vector<unsigned int>({18U, 4U}),
Point<2>(0.7, 0),
Point<2>(2.5, 0.41),
false);
// create middle part first as a hyper shell
GridGenerator::hyper_shell(middle, Point<2>(0.5, 0.2), 0.05, 0.2, 4, true);
middle.set_manifold(0, boundary);
middle.refine_global(1);
// then move the vertices to the points where we want them to be to create a
// slightly asymmetric cube with a hole
for (Triangulation<2>::cell_iterator cell = middle.begin();
cell != middle.end();
++cell)
for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
{
Point<2> &vertex = cell->vertex(v);
if (std::abs(vertex[0] - 0.7) < 1e-10 &&
std::abs(vertex[1] - 0.2) < 1e-10)
vertex = Point<2>(0.7, 0.205);
else if (std::abs(vertex[0] - 0.6) < 1e-10 &&
std::abs(vertex[1] - 0.3) < 1e-10)
vertex = Point<2>(0.7, 0.41);
else if (std::abs(vertex[0] - 0.6) < 1e-10 &&
std::abs(vertex[1] - 0.1) < 1e-10)
vertex = Point<2>(0.7, 0);
else if (std::abs(vertex[0] - 0.5) < 1e-10 &&
std::abs(vertex[1] - 0.4) < 1e-10)
vertex = Point<2>(0.5, 0.41);
else if (std::abs(vertex[0] - 0.5) < 1e-10 &&
std::abs(vertex[1] - 0.0) < 1e-10)
vertex = Point<2>(0.5, 0.0);
else if (std::abs(vertex[0] - 0.4) < 1e-10 &&
std::abs(vertex[1] - 0.3) < 1e-10)
vertex = Point<2>(0.3, 0.41);
else if (std::abs(vertex[0] - 0.4) < 1e-10 &&
std::abs(vertex[1] - 0.1) < 1e-10)
vertex = Point<2>(0.3, 0);
else if (std::abs(vertex[0] - 0.3) < 1e-10 &&
std::abs(vertex[1] - 0.2) < 1e-10)
vertex = Point<2>(0.3, 0.205);
else if (std::abs(vertex[0] - 0.56379) < 1e-4 &&
std::abs(vertex[1] - 0.13621) < 1e-4)
vertex = Point<2>(0.59, 0.11);
else if (std::abs(vertex[0] - 0.56379) < 1e-4 &&
std::abs(vertex[1] - 0.26379) < 1e-4)
vertex = Point<2>(0.59, 0.29);
else if (std::abs(vertex[0] - 0.43621) < 1e-4 &&
std::abs(vertex[1] - 0.13621) < 1e-4)
vertex = Point<2>(0.41, 0.11);
else if (std::abs(vertex[0] - 0.43621) < 1e-4 &&
std::abs(vertex[1] - 0.26379) < 1e-4)
vertex = Point<2>(0.41, 0.29);
}
// refine once to create the same level of refinement as in the
// neighboring domains
middle.refine_global(1);
// must copy the triangulation because we cannot merge triangulations with
// refinement...
GridGenerator::flatten_triangulation(middle, tmp2);
GridGenerator::merge_triangulations(tmp2, right, tria);
// Set the left boundary (inflow) to 1, the right to 2, the rest to 0.
for (Triangulation<2>::active_cell_iterator cell = tria.begin();
cell != tria.end();
++cell)
for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary())
{
if (std::abs(cell->face(f)->center()[0] - 0.3) < 1e-12)
cell->face(f)->set_all_boundary_ids(0);
else if (std::abs(cell->face(f)->center()[0] - 2.5) < 1e-12)
cell->face(f)->set_all_boundary_ids(2); // outlet is 2
else if (Point<2>(0.5, 0.2).distance(cell->face(f)->center()) <=
0.05)
{
cell->face(f)->set_all_manifold_ids(10);
cell->face(f)->set_all_boundary_ids(1);
}
else
cell->face(f)->set_all_boundary_ids(0);
}
}
class Time
{
public:
Time(const double time_end, const double delta_t)
: timestep(0), time_current(0.0), time_end(time_end), delta_t(delta_t)
{
}
virtual ~Time() {}
double current() const { return time_current; }
double end() const { return time_end; }
double get_delta_t() const { return delta_t; }
unsigned int get_timestep() const { return timestep; }
void increment()
{
time_current += delta_t;
++timestep;
}
private:
unsigned int timestep;
double time_current;
const double time_end;
const double delta_t;
};
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues() : Function<dim>(dim + 1) {}
virtual double value(const Point<dim> &p,
const unsigned int component) const;
virtual void vector_value(const Point<dim> &p,
Vector<double> &values) const;
};
template <int dim>
double BoundaryValues<dim>::value(const Point<dim> &p,
const unsigned int component) const
{
Assert(component < this->n_components,
ExcIndexRange(component, 0, this->n_components));
if (component == 0 && std::abs(p[0] - 0.3) < 1e-10)
{
double U = 1.5;
double y = p[1];
return 4 * U * y * (0.41 - y) / (0.41 * 0.41);
}
return 0;
}
template <int dim>
void BoundaryValues<dim>::vector_value(const Point<dim> &p,
Vector<double> &values) const
{
for (unsigned int c = 0; c < this->n_components; ++c)
values(c) = BoundaryValues<dim>::value(p, c);
}
template <int dim>
struct InnerPreconditioner;
template <>
struct InnerPreconditioner<2>
{
typedef SparseDirectUMFPACK type;
//typedef SparseILU<double> type;
};
template <>
struct InnerPreconditioner<3>
{
typedef SparseILU<double> type;
};
// This is used for S inverse and A inverse, which are symmetric so we use CG
// inside
template <class MatrixType, class PreconditionerType>
class InverseMatrix : public Subscriptor
{
public:
InverseMatrix(const MatrixType &m,
const PreconditionerType &preconditioner);
void vmult(Vector<double> &dst, const Vector<double> &src) const;
private:
const SmartPointer<const MatrixType> matrix;
const SmartPointer<const PreconditionerType> preconditioner;
};
template <class MatrixType, class PreconditionerType>
InverseMatrix<MatrixType, PreconditionerType>::InverseMatrix(
const MatrixType &m, const PreconditionerType &preconditioner)
: matrix(&m), preconditioner(&preconditioner)
{
}
template <class MatrixType, class PreconditionerType>
void InverseMatrix<MatrixType, PreconditionerType>::vmult(
Vector<double> &dst, const Vector<double> &src) const
{
SolverControl solver_control(src.size(), 1e-6 * src.l2_norm());
SolverCG<> cg(solver_control);
dst = 0;
cg.solve(*matrix, dst, src, *preconditioner);
}
// Schur complement for the velocity mass matrix,
// Sm = B Mu^{-1} B^T, we approximate Sm with B (diag Mu)^{-1} B^T
class ApproximateMassSchur : public Subscriptor
{
public:
ApproximateMassSchur(const BlockSparseMatrix<double> &M);
void vmult(Vector<double> &dst, const Vector<double> &src) const;
private:
const SmartPointer<const BlockSparseMatrix<double>> mass_matrix;
mutable Vector<double> tmp1, tmp2;
};
ApproximateMassSchur::ApproximateMassSchur(
const BlockSparseMatrix<double> &M)
: mass_matrix(&M), tmp1(M.block(0, 0).m()), tmp2(M.block(0, 0).m())
{
}
void ApproximateMassSchur::vmult(Vector<double> &dst,
const Vector<double> &src) const
{
mass_matrix->block(0, 1).vmult(tmp1, src);
mass_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
mass_matrix->block(1, 0).vmult(dst, tmp2);
}
// Inverse of the total Schur complement which is the sum of the inverse of
// diffusion, Grad-Div term, velocity mass Schur complements
// S^{-1} = dt*(nu + gamma)Mp^{-1} + {B diag{Mu}^{-1} B^T}^{-1}
template <class PreconditionerSm, class PreconditionerMp>
class SchurComplementInverse : public Subscriptor
{
public:
SchurComplementInverse(
double gamma, double viscosity, double dt,
const InverseMatrix<ApproximateMassSchur, PreconditionerSm> &Sm_inv,
const InverseMatrix<SparseMatrix<double>, PreconditionerMp> &Mp_inv);
void vmult(Vector<double> &dst, const Vector<double> &src) const;
private:
const double gamma;
const double viscosity;
const double dt;
const SmartPointer<const InverseMatrix<ApproximateMassSchur,
PreconditionerSm>> Sm_inverse;
const SmartPointer<const InverseMatrix<SparseMatrix<double>,
PreconditionerMp>> Mp_inverse;
};
template <class PreconditionerSm, class PreconditionerMp>
SchurComplementInverse<PreconditionerSm, PreconditionerMp>::SchurComplementInverse(
double gamma, double viscosity, double dt,
const InverseMatrix<ApproximateMassSchur, PreconditionerSm> &Sm_inv,
const InverseMatrix<SparseMatrix<double>, PreconditionerMp> &Mp_inv) :
gamma(gamma), viscosity(viscosity), dt(dt), Sm_inverse(&Sm_inv), Mp_inverse(&Mp_inv)
{
}
template <class PreconditionerSm, class PreconditionerMp>
void SchurComplementInverse<PreconditionerSm, PreconditionerMp>::vmult(
Vector<double> &dst, const Vector<double> &src) const
{
Vector<double> tmp(src.size());
Sm_inverse->vmult(dst, src);
Mp_inverse->vmult(tmp, src);
tmp *= (viscosity + gamma) * dt;
dst += tmp;
}
// PreconditionerA is needed for {\tilde{A}}^{-1}, PreconditionerSm and PreconditionerMp are for {\tilde{S}}^{-1}
template <class PreconditionerA, class PreconditionerSm, class PreconditionerMp>
class BlockSchurPreconditioner : public Subscriptor
{
public:
BlockSchurPreconditioner(
const BlockSparseMatrix<double> &system_m,
const InverseMatrix<SparseMatrix<double>, PreconditionerA> &A_inv,
//const PreconditionerA &A_inv,
const SchurComplementInverse<PreconditionerSm, PreconditionerMp> &S_inv);
void vmult(BlockVector<double> &dst, const BlockVector<double> &src) const;
private:
const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
const SmartPointer<
const InverseMatrix<SparseMatrix<double>, PreconditionerA>> A_inverse;
//const SmartPointer<const PreconditionerA> A_inverse;
const SmartPointer<
const SchurComplementInverse<PreconditionerSm, PreconditionerMp>> S_inverse;
mutable Vector<double> tmp;
};
template <class PreconditionerA, class PreconditionerSm, class PreconditionerMp>
BlockSchurPreconditioner<PreconditionerA, PreconditionerSm, PreconditionerMp>::
BlockSchurPreconditioner(
const BlockSparseMatrix<double> &system_m,
const InverseMatrix<SparseMatrix<double>, PreconditionerA> &A_inv,
//const PreconditionerA &A_inv,
const SchurComplementInverse<PreconditionerSm, PreconditionerMp> &S_inv)
: system_matrix(&system_m), A_inverse(&A_inv), S_inverse(&S_inv),
tmp(system_matrix->block(1, 1).m())
{
}
template <class PreconditionerA, class PreconditionerSm, class PreconditionerMp>
void BlockSchurPreconditioner<PreconditionerA, PreconditionerSm, PreconditionerMp>::vmult(
BlockVector<double> &dst, const BlockVector<double> &src) const
{
A_inverse->vmult(dst.block(0), src.block(0));
system_matrix->block(1, 0).residual(tmp, dst.block(0), src.block(1));
tmp *= -1;
S_inverse->vmult(dst.block(1), tmp);
}
template <int dim>
class NavierStokes
{
public:
NavierStokes(const unsigned int degree);
void run();
private:
void setup();
void assemble(bool assemble_lhs);
std::pair<unsigned int, double> solve_linear_system(bool update_preconditioner);
void output_results(const unsigned int index) const;
void process_solution(std::ofstream& out) const;
const ConstraintMatrix &get_constraints() const;
double viscosity;
double gamma;
const unsigned int degree;
std::vector<types::global_dof_index> dofs_per_block;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
QGauss<dim> quadrature_formula;
QGauss<dim-1> face_quadrature_formula;
ConstraintMatrix zero_constraints;
ConstraintMatrix nonzero_constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockSparseMatrix<double> mass_matrix; // both velocity mass and pressure mass
BlockVector<double> solution;
BlockVector<double> solution_increment;
BlockVector<double> system_rhs;
Time time;
mutable TimerOutput timer;
std::shared_ptr<ApproximateMassSchur> approximate_Sm;
std::shared_ptr<PreconditionIdentity> preconditioner_Sm; // try Jacobi later
std::shared_ptr<InverseMatrix<ApproximateMassSchur, PreconditionIdentity>> Sm_inverse;
std::shared_ptr<SparseILU<double>> preconditioner_Mp; // try ILU later
std::shared_ptr<InverseMatrix<SparseMatrix<double>, SparseILU<double>>> Mp_inverse;
std::shared_ptr<SchurComplementInverse<PreconditionIdentity,
SparseILU<double>>> S_inverse;
std::shared_ptr<typename InnerPreconditioner<dim>::type> preconditioner_A;
std::shared_ptr<InverseMatrix<SparseMatrix<double>,
typename InnerPreconditioner<dim>::type>> A_inverse;
std::shared_ptr<BlockSchurPreconditioner
<typename InnerPreconditioner<dim>::type, PreconditionIdentity, SparseILU<double>>> preconditioner;
};
template <int dim>
NavierStokes<dim>::NavierStokes(const unsigned int degree)
: viscosity(0.001),
gamma(1),
degree(degree),
triangulation(Triangulation<dim>::maximum_smoothing),
fe(FE_Q<dim>(degree + 1), dim, FE_Q<dim>(degree), 1),
dof_handler(triangulation),
quadrature_formula(degree+2),
face_quadrature_formula(degree+2),
time(3e-3, 1e-3),
timer(std::cout, TimerOutput::summary, TimerOutput::wall_times)
{
}
template <int dim>
void NavierStokes<dim>::setup()
{
timer.enter_subsection("Setup system");
// The first step is to associate DoFs with a given mesh.
dof_handler.distribute_dofs(fe);
DoFRenumbering::Cuthill_McKee(dof_handler);
// We renumber the components to have all velocity DoFs come before
// the pressure DoFs to be able to split the solution vector in two blocks
// which are separately accessed
std::vector<unsigned int> block_component(dim + 1, 0);
block_component[dim] = 1;
DoFRenumbering::component_wise(dof_handler, block_component);
dofs_per_block.resize(2);
DoFTools::count_dofs_per_block(
dof_handler, dofs_per_block, block_component);
unsigned int dof_u = dofs_per_block[0];
unsigned int dof_p = dofs_per_block[1];
FEValuesExtractors::Vector velocities(0);
{
nonzero_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, nonzero_constraints);
VectorTools::interpolate_boundary_values(dof_handler,
0,
BoundaryValues<dim>(),
nonzero_constraints,
fe.component_mask(velocities));
VectorTools::interpolate_boundary_values(dof_handler,
1,
BoundaryValues<dim>(),
nonzero_constraints,
fe.component_mask(velocities));
}
nonzero_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));
VectorTools::interpolate_boundary_values(
dof_handler,
1,
Functions::ZeroFunction<dim>(dim + 1),
zero_constraints,
fe.component_mask(velocities));
}
zero_constraints.close();
std::cout << " Number of active cells: " << triangulation.n_active_cells()
<< std::endl
<< " Number of vertices: " << triangulation.n_vertices()
<< std::endl
<< " Number of degrees of freedom: " << dof_handler.n_dofs()
<< " (" << dof_u << '+' << dof_p << ')' << std::endl;
BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
DoFTools::make_sparsity_pattern(dof_handler, dsp, nonzero_constraints);
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
mass_matrix.reinit(sparsity_pattern);
solution.reinit(dofs_per_block);
solution_increment.reinit(dofs_per_block);
system_rhs.reinit(dofs_per_block);
timer.leave_subsection();
}
template <int dim>
const ConstraintMatrix &NavierStokes<dim>::get_constraints() const
{
return time.get_timestep() == 0 ? nonzero_constraints : zero_constraints;
}
template <int dim>
void NavierStokes<dim>::assemble(bool assemble_lhs)
{
timer.enter_subsection("Assemble system");
if (assemble_lhs)
{
system_matrix = 0;
mass_matrix = 0;
}
system_rhs = 0;
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_quadrature_points |
update_JxW_values | update_gradients);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
FullMatrix<double> local_mass_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);
std::vector<Tensor<1, dim>> current_velocity_values(n_q_points);
std::vector<Tensor<2, dim>> current_velocity_gradients(n_q_points);
std::vector<double> current_velocity_divergences(n_q_points);
std::vector<double> current_pressure_values(n_q_points);
std::vector<double> div_phi_u(dofs_per_cell);
std::vector<Tensor<1, dim>> phi_u(dofs_per_cell);
std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
std::vector<double> phi_p(dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator cell =
dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell != endc; ++cell)
{
fe_values.reinit(cell);
local_matrix = 0;
local_rhs = 0;
local_mass_matrix = 0;
fe_values[velocities].get_function_values(solution,
current_velocity_values);
fe_values[velocities].get_function_gradients(
solution, current_velocity_gradients);
fe_values[velocities].get_function_divergences(
solution, current_velocity_divergences);
fe_values[pressure].get_function_values(solution,
current_pressure_values);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
div_phi_u[k] = fe_values[velocities].divergence(k, q);
grad_phi_u[k] = fe_values[velocities].gradient(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)
{
if (assemble_lhs)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
// LHS = a((u, p), (v, q))*dt + m(u, v)
// = (grad_v, nu*grad_u) - (div_v, p) - (q, div_u) +
// m(u, v)
local_matrix(i, j) +=
((viscosity *
scalar_product(grad_phi_u[j], grad_phi_u[i]) -
div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j] +
gamma*div_phi_u[j]*div_phi_u[i]) *
time.get_delta_t() +
phi_u[i] * phi_u[j]) *
fe_values.JxW(q);
// M = [Mu, B^T; B, Mp]
local_mass_matrix(i, j) +=
(phi_u[i] * phi_u[j] + phi_p[i] * phi_p[j] -
div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
fe_values.JxW(q);
}
}
// RHS = - dt*[ a((u_prev, p_prev), (v, q)) + c(u_prev; u_prev,
// v)]
local_rhs(i) -=
(viscosity * scalar_product(current_velocity_gradients[q],
grad_phi_u[i]) -
current_velocity_divergences[q] * phi_p[i] -
current_pressure_values[q] * div_phi_u[i] +
current_velocity_gradients[q] * current_velocity_values[q] *
phi_u[i] +
gamma * current_velocity_divergences[q] * div_phi_u[i]) *
fe_values.JxW(q) * time.get_delta_t();
}
}
cell->get_dof_indices(local_dof_indices);
const ConstraintMatrix &constraints_used = get_constraints();
if (assemble_lhs)
{
constraints_used.distribute_local_to_global(local_matrix,
local_rhs,
local_dof_indices,
system_matrix,
system_rhs);
constraints_used.distribute_local_to_global(local_mass_matrix,
local_dof_indices,
mass_matrix);
}
else
{
constraints_used.distribute_local_to_global(
local_rhs, local_dof_indices, system_rhs);
}
}
timer.leave_subsection();
}
template <int dim>
std::pair<unsigned int, double> NavierStokes<dim>::solve_linear_system(bool update_precondition)
{
const ConstraintMatrix &constraints_used = get_constraints();
if (update_precondition)
{
timer.enter_subsection("Precondition linear system");
preconditioner.reset();
A_inverse.reset();
preconditioner_A.reset();
S_inverse.reset();
Mp_inverse.reset();
preconditioner_Mp.reset();
Sm_inverse.reset();
preconditioner_Sm.reset();
approximate_Sm.reset();
approximate_Sm.reset(new ApproximateMassSchur(mass_matrix));
preconditioner_Sm.reset(new PreconditionIdentity());
Sm_inverse.reset(new InverseMatrix<ApproximateMassSchur, PreconditionIdentity>
(*approximate_Sm, *preconditioner_Sm));
preconditioner_Mp.reset(new SparseILU<double>());
preconditioner_Mp->initialize(mass_matrix.block(1,1));
Mp_inverse.reset(new InverseMatrix<SparseMatrix<double>, SparseILU<double>>
(mass_matrix.block(1,1), *preconditioner_Mp));
S_inverse.reset(new SchurComplementInverse<PreconditionIdentity,
SparseILU<double>>(gamma, viscosity, time.get_delta_t(), *Sm_inverse, *Mp_inverse));
preconditioner_A.reset(new typename InnerPreconditioner<dim>::type());
preconditioner_A->initialize(system_matrix.block(0,0),
typename InnerPreconditioner<dim>::type::AdditionalData());
A_inverse.reset(new InverseMatrix<SparseMatrix<double>,
typename InnerPreconditioner<dim>::type>(system_matrix.block(0,0), *preconditioner_A));
preconditioner.reset(new BlockSchurPreconditioner<
typename InnerPreconditioner<dim>::type, PreconditionIdentity,
SparseILU<double>>(system_matrix, *A_inverse, *S_inverse));
timer.leave_subsection();
}
// Solve
timer.enter_subsection("Solve linear system");
SolverControl solver_control(system_matrix.m(),
1e-8 * system_rhs.l2_norm());
GrowingVectorMemory<BlockVector<double>> vector_memory;
SolverGMRES<BlockVector<double>>::AdditionalData gmres_data;
gmres_data.max_n_tmp_vectors = 100;
SolverGMRES<BlockVector<double>> gmres(
solver_control, vector_memory, gmres_data);
gmres.solve(system_matrix, solution_increment, system_rhs, *preconditioner);
constraints_used.distribute(solution_increment);
timer.leave_subsection();
return {solver_control.last_step(), solver_control.last_value()};
}
template <int dim>
void NavierStokes<dim>::run()
{
create_triangulation_2d(triangulation);
triangulation.refine_global(2);
setup();
std::ofstream out("grid.eps");
GridOut grid_out;
grid_out.write_eps(triangulation, out);
std::ofstream out2("force.txt");
out2 << std::setw(13) << std::left << "Time/s"
<< std::setw(13) << std::left << " Drag" << std::setw(13)
<< std::left << " Lift" << std::endl;
// Without newton's method, this function is simple:
// 1. Solve for the increment;
// 2. Update the solution.
output_results(time.get_timestep());
while (time.current() <= time.end())
{
std::cout << "*****************************************" << std::endl;
std::cout << "Time = " << time.current() << std::endl;
assemble(time.get_timestep() < 2);
auto state = solve_linear_system(time.get_timestep() < 2);
solution.add(1.0, solution_increment);
// solution needs to be distributed with nonzero_constraints all the
// time
nonzero_constraints.distribute(solution);
solution_increment = 0;
std::cout << " FGMRES steps = " << state.first
<< " residual = " << std::setw(6) << state.second << std::endl;
time.increment();
if (time.get_timestep() % 1 == 0)
{
output_results(time.get_timestep());
process_solution(out2);
}
}
out2.close();
}
template <int dim>
void NavierStokes<dim>::output_results(const unsigned int output_index) const
{
timer.enter_subsection("Output");
std::cout << " Writing results..." << std::endl;
std::vector<std::string> solution_names(dim, "velocity");
solution_names.push_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::ostringstream filename;
filename << "Re100-"
<< Utilities::int_to_string(output_index, 6) << ".vtu";
std::ofstream output(filename.str().c_str());
data_out.write_vtu(output);
timer.leave_subsection();
}
template <int dim>
void NavierStokes<dim>::process_solution(std::ofstream& out) const
{
// Compute the drag and lift coefficients of the cylinder
timer.enter_subsection("Process solution");
Tensor<1, dim> force;
FEFaceValues<dim> fe_face_values(fe,
face_quadrature_formula,
update_values | update_quadrature_points |
update_JxW_values | update_normal_vectors |
update_gradients);
const unsigned int n_q_points = face_quadrature_formula.size();
const FEValuesExtractors::Vector velocities(0);
const FEValuesExtractors::Scalar pressure(dim);
std::vector<double> p(n_q_points);
std::vector<SymmetricTensor<2, dim>> grad_sym_v(n_q_points);
for (auto cell = dof_handler.begin_active(); cell != dof_handler.end(); ++cell)
{
for (unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; ++f)
{
if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1)
{
fe_face_values.reinit(cell, f);
fe_face_values[pressure].get_function_values(solution, p);
fe_face_values[velocities].get_function_symmetric_gradients(solution, grad_sym_v);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const Tensor<1, dim> &N = fe_face_values.normal_vector(q);
SymmetricTensor<2, dim> stress = -p[q]*Physics::Elasticity::StandardTensors<dim>::I
+ viscosity*grad_sym_v[q];
force -= stress*N*fe_face_values.JxW(q);
}
}
}
}
double drag_coef = 2*force[0]/(0.1);
double lift_coef = 2*force[dim-1]/(0.1);
out.precision(6);
out.width(12);
out << std::scientific << std::left <<
time.current() << " " << drag_coef << " " << lift_coef << std::endl;
timer.leave_subsection();
}
}
int main()
{
try
{
using namespace dealii;
using namespace fluid;
NavierStokes<2> flow(/* degree = */ 1);
flow.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;
}
