i am using the implimenation by Jaekwang Kim to model polymer flows. The
major change that I want to do is to output convergence tables like in
Step-7 and Step-22, which means I have to have an ExactSolution.
The code is encountering an error because the dimension of the solutions do
not match at a some point when generating output:
*An error occurred in line <820> of file
</home/jurombo/Downloads/dealii-9.3.2/include/deal.II/numerics/vector_tools_integrate_difference.templates.h>
in function void
dealii::VectorTools::internal::do_integrate_difference(const
dealii::hp::MappingCollection<dim, spacedim>&, const
dealii::DoFHandler<dim, spacedim>&, const InVector&, const
dealii::Function<spacedim, typename InVector::value_type>&, OutVector&,
const dealii::hp::QCollection<dim>&, const dealii::VectorTools::NormType&,
const dealii::Function<spacedim>*, double) [with int dim = 2; int spacedim
= 2; InVector = dealii::BlockVector<double>; OutVector =
dealii::Vector<double>; typename InVector::value_type = double] The
violated condition was: exact_solution.n_components == n_components
Additional information: Dimension 3 not equal to 7. Stacktrace:
----------- #0 /usr/local/share/dealii/lib/libdeal_II.g.so.9.3.2: #1
/usr/local/share/dealii/lib/libdeal_II.g.so.9.3.2: void
dealii::VectorTools::integrate_difference<2, dealii::BlockVector<double>,
dealii::Vector<double>, 2>(dealii::Mapping<2, 2> const&,
dealii::DoFHandler<2, 2> const&, dealii::BlockVector<double> const&,
dealii::Function<2, dealii::BlockVector<double>::value_type> const&,
dealii::Vector<double>&, dealii::Quadrature<2> const&,
dealii::VectorTools::NormType const&, dealii::Function<2, double> const*,
double) #2 ./oldroyd_bn:
Viscoelastic::StokesProblem<2>::output_results(unsigned int) #3
./oldroyd_bn: Viscoelastic::StokesProblem<2>::run() #4 ./oldroyd_bn: main
--------------------------------------------------------*
I have checked the ExactSolution in the code to try and figure wher the
error is occurring without success.
Can anyone point me the direction to look.
--
The information in this message is confidential and legally
privileged.
It is intended solely for the addressee(s). Access to this message
by
anyone else is unauthorized. If received in error, please accept our
apologies
and notify the sender immediately. You must also delete the
original message
from your machine. If you are not the intended recipient,
any use, disclosure,
copying, distribution or action taken in reliance of
it, is prohibited and may
be unlawful. The information, attachments,
opinions or advice contained in this
email are not the views or opinions of
Harare Institute of Technology, its
subsidiaries or affiliates. Although
this email and any attachments are
believed to be free of any virus or
other defects which might affect any
computer or IT system into which they
are received, no responsibility is
accepted by Harare Institute of
Technology and/or its subsidiaries for any loss
or damage arising in any
way from the receipt or use thereof.
--
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/a7464d4b-87df-4e49-80cf-7a0900fe1a70n%40googlegroups.com.
// Developed by Jaekwang Kim
// Welcome to the world of viscoelastic fluids
// Coding from 18th May ~
// Viscoelastic Models
// Oldroyd B model
// Upper Convective Maxwel polymer(UCM) + a Newtonian solvent
// Mixed Finite Element (also called viscous formulation) will be considered
// (a) Starting from large viscous contribution from a Newtonian Solvent
// (b) Low Wissenberg Number
// (c) Stabilization Method
// (d) Be careful on the LBB condition; Continuous Discretization of T_tensor?
//what is the convention of grad_u_[0][1] ?
#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/convergence_table.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> // CG used for viscous flow
#include <deal.II/lac/solver_gmres.h> //Matrix solver for transport equation, unsymmetric matrix
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.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/fe/mapping_q.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/base/parameter_handler.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include "headers/post_processor.h"
#include "headers/precondition.h"
#include "headers/bcfunctions.h"
#include "headers/fluidclass.h"
#include "headers/flowclass.h"
namespace Viscoelastic
{
using namespace dealii;
template <int dim>
class StokesProblem
{
public:
StokesProblem (const unsigned int degree);
void run ();
private:
void read_mesh ();
void setup_dofs ();
void setup_bc_constraints ();
void assemble_stokes (); // Flow system
void solve_stokes (); // Solve with CG
void assemble_polymer_system (); // UCM advection
void solve_polymer_system (); // Solve with GMRES
void refine_mesh (const unsigned int refinement_cycle);
void output_results (const unsigned int refinement_cycle);
const unsigned int degree;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
MappingQ<dim> mapping;
AffineConstraints<double> constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockVector<double> solution;
// BlockVector<double> exact_solution;
BlockVector<double> interpolated_exact_solution;
BlockVector<double> system_rhs;
// Previous solution for Piccard iterations
BlockVector<double> previous_solution;
AdvectionBoundaryValues<dim> advection_boundary_values;
ConstantU<dim> Uplate; // BC Function
Oldroyd_FluidClass fluid;
double lam1 = fluid.lam1;
double lam2 = fluid.lam2;
double etap = fluid.etap;
double etas = fluid.etas;
Oldroyd_FlowClass flow;
double Ux = flow.Ux;
double h = flow.h;
double alpha = flow.alpha;
double Wi = flow.Wi;
std::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
ConvergenceTable convergence_table;
};
// @sect3{Boundary values and right hand side}
template <int dim>
class BoundaryValues : public Function<dim>
{
public:
BoundaryValues () : Function<dim>(dim+1) {}
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
const double U = 1.6;
// const double Ux = flow.Ux;
};
template <int dim>
void
BoundaryValues<dim>::vector_value (const Point<dim> & p,
Vector<double> &values) const
{
const double rad = 1.0 - p.square();
const double Ux = .0125;
values(0) = 0.0;
values(1) = 0.0;
values(2) = 0.0;
values(3) = Ux*(1 - rad);
values(4) = 0.0;
values(5) = 0.0;
values(6) = 0.0;
}
// the right hand side
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;
};
template <int dim>
void
RightHandSide<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
values(0) = -8.0*p[1] - 2.0*p[0];
values(1) = 8.0*p[0] - 2.0*p[1];
values(2) = 0.0;
values(3) = 0.0;
values(4) = 0.0;
values(5) = 0.0;
values(6) = 0.0;
}
// ExactSolution
template <int dim>
class ExactSolution : public Function<dim>
{
public:
ExactSolution () : Function<dim>(dim+1) {}
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
Tensor<1,dim> gradient (const Point<dim> &p,
const unsigned int component) const;
const double R = 1.0;
const double alpha = alpha = etas/(etap+etas);
const double Wi = lam1*(Ux/h);
//const double U = 1.6;
//const double Ux = .125;
};
template <int dim>
void
ExactSolution<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
Assert (values.size() == dim+1,
ExcDimensionMismatch (values.size(), dim+1));
const double rad = 1.0 - p.square();
values(0) = -p[1]*rad;
values(1) = p[0]*rad;
values(2) = rad;
values(3) = -2*p[1]*(1 - alpha)*U/R*R;
values(4) = 8*Wi*(1 - alpha)*p.square()*Ux*Ux/R*R;
values(5) = 8*Wi*(1 - alpha)*p.square()*Ux*Ux/R*R;
values(6) = 0.0;
}
template <int dim>
Tensor<1,dim>
ExactSolution<dim>::gradient (const Point<dim> &p,
const unsigned int component) const
{
Assert (component < this->n_components,
ExcIndexRange (component, 0, this->n_components));
double x = p[0];
double y = p[1];
Tensor<1,dim> gradient;
switch (component)
{
//velocity
case 0:
gradient[0]=2*x*y;
gradient[1]=-1+x*x+3*y*y;
break;
case 1:
gradient[0]=1-3*x*x-y*y;
gradient[1]=-2*x*y;
break;
//pressure
case 2:
gradient[0]=-2*x;
gradient[1]=-2*y;
break;
// stress
case 3:
gradient[0]=-2*(1-alpha)*U/R*R;
gradient[1]=16*Wi*(1-alpha)*p[1]*Ux*Ux/R*R;
gradient[2]=16*Wi*(1-alpha)*p[1]*Ux*Ux/R*R;
gradient[3]=0.0;
break;
}
return gradient;
}
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int degree)
:
degree (degree),
triangulation (Triangulation<dim>::allow_anisotropic_smoothing),
fe (FE_Q<dim>(degree+1), dim, // u
FE_Q<dim>(degree) , 1 , // p
FE_Q<dim>(degree+1) , dim*dim ), //tau_p
dof_handler (triangulation),
mapping(2)
{}
template <int dim>
void StokesProblem<dim>::setup_bc_constraints ()
{
constraints.clear ();
FEValuesExtractors::Vector velocities(0);
FEValuesExtractors::Scalar xvel(0);
FEValuesExtractors::Scalar yvel(1);
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
// Inlet
VectorTools::interpolate_boundary_values (dof_handler,
4,
Uinlet<dim>(),
constraints,
fe.component_mask(velocities));
// Ends - parallel flow
VectorTools::interpolate_boundary_values (dof_handler,
2,
ZeroFunction<dim>(dim+1+dim*dim),
constraints,
fe.component_mask(yvel));
// Side - Uniform wall motion:
VectorTools::interpolate_boundary_values (dof_handler,
5,
Uplate,
constraints,
fe.component_mask(velocities));
// Fixed Wall --- 1 is flat sections, > 5 are rounded sections with manifolds
VectorTools::interpolate_boundary_values (dof_handler,
1,
ZeroFunction<dim>(dim+1+dim*dim),
constraints,
fe.component_mask(velocities));
VectorTools::interpolate_boundary_values (dof_handler,
6,
ZeroFunction<dim>(dim+1+dim*dim),
constraints,
fe.component_mask(velocities));
VectorTools::interpolate_boundary_values (dof_handler,
7,
ZeroFunction<dim>(dim+1+dim*dim),
constraints,
fe.component_mask(velocities));
VectorTools::interpolate_boundary_values (dof_handler,
8,
ZeroFunction<dim>(dim+1+dim*dim),
constraints,
fe.component_mask(velocities));
constraints.close ();
}
template <int dim>
void StokesProblem<dim>::setup_dofs ()
{
std::cout << " Set_Up_DOF -- Start" <<std::endl;
A_preconditioner.reset ();
system_matrix.clear ();
dof_handler.distribute_dofs (fe);
DoFRenumbering::Cuthill_McKee (dof_handler);
std::vector<unsigned int> block_component (dim+1+dim*dim,0); // initial 0 for u
block_component[dim] = 1; // pressure
for (unsigned int i=0;i<dim*dim;i++)
block_component[dim+1+i] = 2; //block for polymer_T variable...
DoFRenumbering::component_wise (dof_handler, block_component);
std::vector<types::global_dof_index> dofs_per_block (3);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
const unsigned int
n_u = dofs_per_block[0],
n_p = dofs_per_block[1],
n_s = dofs_per_block[2];
std::cout << " n_u: " << n_u
<< " n_p: " << n_p
<< " n_s: " << n_s << std::endl;
// Apply BC Constraints
setup_bc_constraints();
{
BlockDynamicSparsityPattern dsp (3,3);
dsp.block(0,0).reinit (n_u, n_u);
dsp.block(1,0).reinit (n_p, n_u);
dsp.block(0,1).reinit (n_u, n_p);
dsp.block(1,1).reinit (n_p, n_p);
dsp.block(1,2).reinit (n_p, n_s);
dsp.block(2,1).reinit (n_s, n_p);
dsp.block(2,0).reinit (n_s, n_u);
dsp.block(0,2).reinit (n_u, n_s);
dsp.block(2,2).reinit (n_s, n_s);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, true);
sparsity_pattern.copy_from (dsp);
}
system_matrix.reinit (sparsity_pattern);
solution.reinit (3);
solution.block(0).reinit (n_u);
solution.block(1).reinit (n_p);
solution.block(2).reinit (n_s);
solution.collect_sizes ();
interpolated_exact_solution.reinit (3);
interpolated_exact_solution.block(0).reinit (n_u);
interpolated_exact_solution.block(1).reinit (n_p);
interpolated_exact_solution.block(2).reinit (n_s);
interpolated_exact_solution.collect_sizes ();
system_rhs.reinit (3);
system_rhs.block(0).reinit (n_u);
system_rhs.block(1).reinit (n_p);
system_rhs.block(2).reinit (n_s);
system_rhs.collect_sizes ();
std::cout << " Active cells: "<< triangulation.n_active_cells() << std::endl
<< " DoFs: " << dof_handler.n_dofs() << " (" << n_u << '+' << n_p << '+'<< n_s <<')'
<< std::endl;
}
// Assemble Momentum and Continuity Equation
template <int dim>
void StokesProblem<dim>::assemble_stokes ()
{ // Momentum and Continuity Equation
system_matrix=0;
system_rhs=0;
const MappingQ<dim> mapping (degree);
QGauss<dim> quadrature_formula(2*degree+1);
FEValues<dim> fe_values (mapping, 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();
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);
std::vector<Vector<double> > rhs_values (n_q_points,
Vector<double>(dim+1+dim*dim));
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar Txx (dim+1);
const FEValuesExtractors::Scalar Txy (dim+2);
const FEValuesExtractors::Scalar Tyx (dim+3);
const FEValuesExtractors::Scalar Tyy (dim+4);
std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
std::vector<Tensor<2,dim> > grad_phi_u (dofs_per_cell);
std::vector<double> div_phi_u (dofs_per_cell);
std::vector<double> phi_p (dofs_per_cell);
std::vector<double> phi_ur (dofs_per_cell);
std::vector<double> phi_i_s_xx(dofs_per_cell);
std::vector<double> phi_i_s_xy(dofs_per_cell);
std::vector<double> phi_i_s_yx(dofs_per_cell);
std::vector<double> phi_i_s_yy(dofs_per_cell);
// Variables related to previous solution
std::vector<SymmetricTensor<2,dim> > local_symgrad_u (n_q_points);
std::vector<SymmetricTensor<2,dim> > local_T (n_q_points);
std::vector<double> local_Txx (n_q_points);
std::vector<double> local_Txy (n_q_points);
std::vector<double> local_Tyx (n_q_points);
std::vector<double> local_Tyy (n_q_points);
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;
fe_values[velocities].get_function_symmetric_gradients
(previous_solution, local_symgrad_u);
fe_values[Txx].get_function_values(previous_solution,local_Txx);
fe_values[Txy].get_function_values(previous_solution,local_Txy);
fe_values[Tyx].get_function_values(previous_solution,local_Tyx);
fe_values[Tyy].get_function_values(previous_solution,local_Tyy);
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);
grad_phi_u[k] = fe_values[velocities].gradient (k, q);
div_phi_u[k] = fe_values[velocities].divergence (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) +=
(
//Stokes Formulation
//(a) Solvent Contribution
2.*etas*(symgrad_phi_u[i]*symgrad_phi_u[j])
//(b) presure term and continuity term
- div_phi_u[i]*phi_p[j]
- phi_p[i]*div_phi_u[j]
+ phi_p[i]*phi_p[j]
)
* fe_values.JxW(q);
}
local_rhs(i) +=
( grad_phi_u[i][0][0]*(local_Txx[q])
+grad_phi_u[i][1][0]*(local_Txy[q])
+grad_phi_u[i][0][1]*(local_Tyx[q])
+grad_phi_u[i][1][1]*(local_Tyy[q])
)*fe_values.JxW(q);
}
}
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);
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (local_matrix, local_rhs,
local_dof_indices,
system_matrix, system_rhs);
}
// apply Dirichlet boundary condition on velocity
std::map<types::global_dof_index,double> boundary_values;
A_preconditioner
= std::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
A_preconditioner->initialize (system_matrix.block(0,0),
typename InnerPreconditioner<dim>::type::AdditionalData());
}
// Assemble Constitutive Equation
template <int dim>
void StokesProblem<dim>::assemble_polymer_system ()
{
std::cout << " Assemble Polymer Constitutive " << std::endl;
const MappingQ<dim> mapping (degree);
QGauss<dim> quadrature_formula(2*degree+1);
QGauss<dim-1> face_quadrature_formula(2*degree+1);
FEValues<dim> fe_values (mapping, fe, quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
FEFaceValues<dim> fe_face_values (mapping, fe, face_quadrature_formula,
update_values |
update_normal_vectors |
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();
const unsigned int n_face_q_points = face_quadrature_formula.size();
std::vector<Vector<double> > solution_values_face(n_face_q_points, Vector<double>(dim+1+dim*dim));
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs (dofs_per_cell);
// Flow solution;
std::vector<SymmetricTensor<2,dim> > local_symgrad_phi_u (n_q_points);
std::vector<Tensor<2,dim> > local_grad_phi_u (n_q_points);
std::vector<Tensor<1,dim> > local_phi_u (n_q_points);
std::vector<Tensor<1,1> > local_pres (n_q_points);
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar Txx (dim+1);
const FEValuesExtractors::Scalar Txy (dim+2);
const FEValuesExtractors::Scalar Tyx (dim+3);
const FEValuesExtractors::Scalar Tyy (dim+4);
const FEValuesExtractors::Scalar pres (dim);
std::vector<Tensor<1,dim>> face_advection_directions (n_face_q_points);
std::vector<Vector<double>> face_boundary_values (n_face_q_points,Vector<double>(dim+1+dim*dim));
Tensor<1,dim> normal;
double vel_magnitude;
double phi_i_s_xx, phi_i_s_xy, phi_i_s_yy, phi_i_s_yx;
double phi_j_s_xx, phi_j_s_xy, phi_j_s_yy, phi_j_s_yx;
Tensor<1,dim> grad_phi_i_s_xx, grad_phi_i_s_xy, grad_phi_i_s_yy, grad_phi_i_s_yx;
Tensor<1,dim> grad_phi_j_s_xx, grad_phi_j_s_xy, grad_phi_j_s_yy, grad_phi_j_s_yx;
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;
// Use of local solution
fe_values[velocities].get_function_symmetric_gradients (solution, local_symgrad_phi_u);
fe_values[velocities].get_function_gradients (solution, local_grad_phi_u);
fe_values[velocities].get_function_values (solution, local_phi_u);
// fe_values[pres].get_function_values (solution, local_pres);
// Stabilization
double delta = 0.1 * cell->diameter ();
for (unsigned int q=0; q<n_q_points; ++q)
{
vel_magnitude = std::sqrt(local_phi_u[q]*local_phi_u[q]);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
phi_i_s_xx = fe_values[Txx].value(i,q);
phi_i_s_xy = fe_values[Txy].value(i,q);
phi_i_s_yx = fe_values[Tyx].value(i,q);
phi_i_s_yy = fe_values[Tyy].value(i,q);
grad_phi_i_s_xx = fe_values[Txx].gradient (i, q);
grad_phi_i_s_xy = fe_values[Txy].gradient (i, q);
grad_phi_i_s_yx = fe_values[Tyx].gradient (i, q);
grad_phi_i_s_yy = fe_values[Tyy].gradient (i, q);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
phi_j_s_xx = fe_values[Txx].value(j,q);
phi_j_s_xy = fe_values[Txy].value(j,q);
phi_j_s_yx = fe_values[Tyx].value(j,q);
phi_j_s_yy = fe_values[Tyy].value(j,q);
grad_phi_j_s_xx = fe_values[Txx].gradient (j, q);
grad_phi_j_s_xy = fe_values[Txy].gradient (j, q);
grad_phi_j_s_yx = fe_values[Tyx].gradient (j, q);
grad_phi_j_s_yy = fe_values[Tyy].gradient (j, q);
local_matrix(i,j) +=
(
//self term
+ phi_i_s_xx * phi_j_s_xx
+ phi_i_s_xy * phi_j_s_xy
+ phi_i_s_yx * phi_j_s_yx
+ phi_i_s_yy * phi_j_s_yy
+ lam1*(
//advection
+ phi_i_s_xx * (local_phi_u[q]*grad_phi_j_s_xx )
+ phi_i_s_xy * (local_phi_u[q]*grad_phi_j_s_xy )
+ phi_i_s_yx * (local_phi_u[q]*grad_phi_j_s_yx )
+ phi_i_s_yy * (local_phi_u[q]*grad_phi_j_s_yy )
// rotation term //
// Gradient Tensor convention (dux/dy)=local_grad_phi_u[q][0][1])
- phi_i_s_xx * ( 2.*local_grad_phi_u[q][0][0]*phi_j_s_xx +
local_grad_phi_u[q][0][1]*phi_j_s_yx +
local_grad_phi_u[q][0][1]*phi_j_s_xy )
- phi_i_s_xy * ( local_grad_phi_u[q][1][0]*phi_j_s_xx +
local_grad_phi_u[q][0][0]*phi_j_s_xy + //(a)
local_grad_phi_u[q][1][1]*phi_j_s_xy + //(b)
local_grad_phi_u[q][0][1]*phi_j_s_yy )
- phi_i_s_yx * (local_grad_phi_u[q][1][0]*phi_j_s_xx +
local_grad_phi_u[q][0][0]*phi_j_s_yx + //(c)
local_grad_phi_u[q][1][1]*phi_j_s_yx + //(d)
local_grad_phi_u[q][0][1]*phi_j_s_yy )
//actually (a+b+c+d) will be zero effectively
- phi_i_s_yy * ( 2.*local_grad_phi_u[q][1][1]*phi_j_s_yy +
local_grad_phi_u[q][1][0]*phi_j_s_yx +
local_grad_phi_u[q][1][0]*phi_j_s_xy )
//std::cout << "local_"
// stabilization //need to check
+ delta/vel_magnitude*
(
local_phi_u[q]*grad_phi_i_s_xx //(u \cdot \nabla w)
*local_phi_u[q]*grad_phi_j_s_xx
+
local_phi_u[q]*grad_phi_i_s_xy
*local_phi_u[q]*grad_phi_j_s_xy
+
local_phi_u[q]*grad_phi_i_s_yy
*local_phi_u[q]*grad_phi_j_s_yy
+
local_phi_u[q]*grad_phi_i_s_yx
*local_phi_u[q]*grad_phi_j_s_yx
)
)
)*fe_values.JxW(q);
} //close 'j' cycle
local_rhs(i) +=
etap*(
+ phi_i_s_xx * 2.*local_grad_phi_u[q][0][0]
+ phi_i_s_xy * (local_grad_phi_u[q][1][0] + local_grad_phi_u[q][0][1])
+ phi_i_s_yy * 2.*local_grad_phi_u[q][1][1]
+ phi_i_s_yx * (local_grad_phi_u[q][1][0] + local_grad_phi_u[q][0][1])
)*fe_values.JxW(q);
} // i
} // q
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
if (cell->face(face_no)->boundary_id()==1) // strong inflow boundary condition
{
// inflow face
fe_face_values.reinit (cell, face_no);
fe_face_values.get_function_values (solution,solution_values_face);
advection_boundary_values.vector_value_list
(fe_face_values.get_quadrature_points(),face_boundary_values);
for (unsigned int q=0; q<n_face_q_points; ++q)
{
Tensor<1,dim> present_u_face;
for (unsigned int d=0; d<dim; ++d)
present_u_face[d] = solution_values_face[q](d);
//if-inflow condition
if (fe_face_values.normal_vector(q) * present_u_face < 0)
{
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
phi_i_s_xx = fe_face_values[Txx].value(i,q);
phi_i_s_xy = fe_face_values[Txy].value(i,q);
phi_i_s_yy = fe_face_values[Tyy].value(i,q);
phi_i_s_yx = fe_face_values[Tyx].value(i,q);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
phi_j_s_xx = fe_face_values[Txx].value(j,q);
phi_j_s_xy = fe_face_values[Txy].value(j,q);
phi_j_s_yy = fe_face_values[Tyy].value(j,q);
phi_j_s_yx = fe_face_values[Tyx].value(j,q);
local_matrix(i,j) -= (present_u_face
*fe_face_values.normal_vector(q)
*(+phi_i_s_xx*phi_j_s_xx
+phi_i_s_xy*phi_j_s_xy
+phi_i_s_yy*phi_j_s_yy
+phi_i_s_yx*phi_j_s_yx
)
)*fe_face_values.JxW(q);
} //close cycle j
local_rhs(i) -= present_u_face * fe_face_values.normal_vector(q)*
(
+ phi_i_s_xx * face_boundary_values[q][3]
+ phi_i_s_xy * face_boundary_values[q][4]
+ phi_i_s_yx * face_boundary_values[q][5]
+ phi_i_s_yy * face_boundary_values[q][6]
)*fe_face_values.JxW(q);
} // i
} // if-inflow
} // q
} // cell
} // face
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>::solve_stokes ()
{
const InverseMatrix<SparseMatrix<double>,
typename InnerPreconditioner<dim>::type>
A_inverse (system_matrix.block(0,0), *A_preconditioner);
Vector<double> tmp (solution.block(0).size());
{
Vector<double> schur_rhs (solution.block(1).size());
A_inverse.vmult (tmp, system_rhs.block(0));
system_matrix.block(1,0).vmult (schur_rhs, tmp);
schur_rhs -= system_rhs.block(1);
SchurComplement<typename InnerPreconditioner<dim>::type>
schur_complement (system_matrix, A_inverse);
SolverControl solver_control ( 10 * solution.block(1).size(),
1e-5*schur_rhs.l2_norm()); // solver_control here
SolverCG<> cg (solver_control);
SparseILU<double> preconditioner;
preconditioner.initialize (system_matrix.block(1,1),
SparseILU<double>::AdditionalData());
InverseMatrix<SparseMatrix<double>,SparseILU<double> >
m_inverse (system_matrix.block(1,1), preconditioner);
cg.solve (schur_complement, solution.block(1), schur_rhs,
m_inverse);
constraints.distribute (solution);
std::cout << "Flow Solved:" << solver_control.last_step()
<< " outer CG Schur complement iterations for pressure"
<< std::endl;
}
{
system_matrix.block(0,1).vmult (tmp, solution.block(1));
tmp *= -1;
tmp += system_rhs.block(0);
A_inverse.vmult (solution.block(0), tmp);
constraints.distribute (solution);
}
}
template <int dim>
void StokesProblem<dim>::solve_polymer_system ()
{
std::cout << " Solve Advection" << std::endl;
SolverControl solver_control (std::pow(10,8), std::pow(10,-4));
unsigned int restart = 500;
SolverGMRES< Vector<double> >::AdditionalData gmres_additional_data(restart+2);
SolverGMRES< Vector<double> > solver(solver_control, gmres_additional_data);
//'gmres_additional_data' means how much temporary vector you will going to use for grmres solver
//the more the faster, but the more memory will be consummed.
//make preconditioner
SparseILU<double>::AdditionalData additional_data(0,1000); // (0 , additional diagonal terms)
SparseILU<double> preconditioner;
preconditioner.initialize (system_matrix.block(2,2), additional_data);
solver.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), preconditioner);
constraints.distribute (solution);
}
template <int dim>
void
StokesProblem<dim>::output_results (const unsigned int refinement_cycle)
{
std::vector<std::string> solution_names(dim,"Velocity");
solution_names.push_back ("Pressure");
solution_names.push_back ("Txx");
solution_names.push_back ("Txy");
solution_names.push_back ("Tyx");
solution_names.push_back ("Tyy");
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation
(dim, DataComponentInterpretation::component_is_part_of_vector);
data_component_interpretation
.push_back (DataComponentInterpretation::component_is_scalar);
for (unsigned int i=0; i<dim*dim; i++)
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 filenameeps;
filenameeps << "solution-"<< Utilities::int_to_string (refinement_cycle, 2)<< ".vtk";
std::ofstream output (filenameeps.str().c_str());
data_out.write_vtk (output);
//data_out.write_tecplot (output);
// Compute error
const ComponentSelectFunction<dim>
pressure_mask (dim, dim+1);
const ComponentSelectFunction<dim>
velocity_mask(std::make_pair(0, dim), dim+1);
ExactSolution<dim> exact_solution;
Vector<double> cellwise_errors (triangulation.n_active_cells());
// QTrapez<1> q_trapez;
// QIterated<dim> quadrature (q_trapez, degree+2);
QGauss<dim> quadrature (degree+3);
const QTrapezoid<1> q_trapez;
const QIterated<dim> q_iterated (q_trapez, 5);
// L2 error (velocity)
VectorTools::integrate_difference (mapping, dof_handler,
solution,
exact_solution,
cellwise_errors,
quadrature,
VectorTools::L2_norm,
&velocity_mask);
const double u_L2_error = cellwise_errors.l2_norm();
// H1 semi norm (velocity)
VectorTools::integrate_difference (mapping, dof_handler,
solution,
exact_solution,
cellwise_errors,
quadrature,
VectorTools::H1_seminorm,
&velocity_mask);
const double u_H1_error = cellwise_errors.l2_norm();
// sup norm (velocity)
VectorTools::integrate_difference (mapping, dof_handler,
solution,
exact_solution,
cellwise_errors,
q_iterated,
VectorTools::Linfty_norm,
&velocity_mask);
const double u_Linfty_error = cellwise_errors.linfty_norm();
// L2 error (pressure)
VectorTools::integrate_difference (mapping, dof_handler,
solution,
exact_solution,
cellwise_errors,
quadrature,
VectorTools::L2_norm,
&pressure_mask);
const double p_L2_error = cellwise_errors.l2_norm();
// SUP error (pressure)
VectorTools::integrate_difference (mapping, dof_handler,
solution,
exact_solution,
cellwise_errors,
q_iterated,
VectorTools::Linfty_norm,
&pressure_mask);
const double p_Linfty_error = cellwise_errors.linfty_norm();
const unsigned int n_active_cells=triangulation.n_active_cells();
const unsigned int n_dofs=dof_handler.n_dofs();
convergence_table.add_value("cycle", refinement_cycle);
convergence_table.add_value("cells", n_active_cells);
convergence_table.add_value("dofs", dof_handler.n_dofs());
convergence_table.add_value("u_L2", u_L2_error);
convergence_table.add_value("u_H1", u_H1_error);
convergence_table.add_value("u_Linfty", u_Linfty_error);
convergence_table.add_value("p_L2", p_L2_error);
convergence_table.add_value("p_Linfty", p_Linfty_error);
// output numerical solution
DataOut<2> sol_out;
sol_out.attach_dof_handler (dof_handler);
sol_out.add_data_vector (solution, "numerical solution");
sol_out.build_patches ();
std::ofstream sol_output ("num_sol.gpl");
sol_out.write_gnuplot (sol_output);
// output nodal error
if ( refinement_cycle == 0 ) {
std::cout << "Output nodal error" << std::endl;
VectorTools::interpolate (dof_handler,
ExactSolution<dim>(),
interpolated_exact_solution);
interpolated_exact_solution -= solution;
DataOut<2> err_out;
err_out.attach_dof_handler (dof_handler);
err_out.add_data_vector (interpolated_exact_solution, "pointwise_error");
err_out.build_patches ();
std::ofstream err_output ("pointwise_err.gpl");
err_out.write_gnuplot (err_output);
}
Vector<double> interpolated_exact_solution (dof_handler.n_dofs());
VectorTools::interpolate (dof_handler,
ExactSolution<dim>(),
interpolated_exact_solution);
for ( unsigned int i=0; i< solution.block(0).size(); i++ )
interpolated_exact_solution(i) -= solution.block(0)(i);
DataOut<2> err_out;
err_out.attach_dof_handler (dof_handler);
err_out.add_data_vector (interpolated_exact_solution, "pointwise_error");
err_out.build_patches ();
std::ofstream err_output ("pointwise_err.gpl");
err_out.write_gnuplot (err_output);
// output result
std::cout << " Number of active cells: "
<< n_active_cells
<< std::endl
<< " Number of degrees of freedom: "
<< n_dofs
<< std::endl;
std::cout << "p L2 = " << p_L2_error << std::endl;
std::cout << "p SUP= " << p_Linfty_error << std::endl;
std::cout << "u L2 = " << u_L2_error << std::endl;
std::cout << "u SUP= " << u_Linfty_error << std::endl;
std::cout << "u H1 = " << u_H1_error << std::endl;
}
template <int dim>
void StokesProblem<dim>::read_mesh ()
{
{ //Generate mesh and designate boundary
/* Parallel plate
const Point<2> end (5,1);
const Point<2> start (0.0,0.0);
GridGenerator::hyper_rectangle (triangulation, start, end, false);
triangulation.refine_global (3);
*/
GridIn<dim> grid_in;
grid_in.attach_triangulation (triangulation);
//std::ifstream input_file("mesh/slant-rounded-mesh.inp");
std::ifstream input_file("mesh/new_mesh.inp");
std::cout <<"here" << std::endl;
grid_in.read_ucd (input_file);
std::cout <<"here" << std::endl;
//Point<2> center1 (-0.2,-1.7);
Point<2> center1 (-1.7,-0.2);
static const SphericalManifold<dim> manifold_description1 (center1);
triangulation.set_manifold (6, manifold_description1);
//Point<2> center2 (-1,-1.3);
Point<2> center2 (-1.3,-1.0);
static const SphericalManifold<dim> manifold_description2 (center2);
triangulation.set_manifold (7, manifold_description2);
//Point<2> center3 (-0.6,1.5);
Point<2> center3 (1.5,-0.6);
static const SphericalManifold<dim> manifold_description3 (center3);
triangulation.set_manifold (8, manifold_description3);
triangulation.refine_global (3);
}
}
template <int dim>
void StokesProblem<dim>::run ()
{
read_mesh(); //To read mesh file from outside
for (unsigned int refinement_cycle = 0; refinement_cycle<3; ++refinement_cycle)
{
std::cout << "Refinement cycle " << refinement_cycle << std::endl;
setup_dofs ();
if (refinement_cycle < 1 )
{
previous_solution = solution;
previous_solution = 0.;
}
assemble_stokes ();
solve_stokes ();
assemble_polymer_system ();
solve_polymer_system ();
BlockVector<double> difference;
int iteration_number=0 ;
previous_solution =solution ;
do{
iteration_number +=1;
assemble_stokes ();
solve_stokes ();
assemble_polymer_system ();
solve_polymer_system ();
difference = solution;
difference -= previous_solution;
previous_solution=solution;
std::cout << " Iteration Number: " << iteration_number
<< " Difference Norm: " << difference.l2_norm()
<< std::endl << std::flush;
} while (difference.l2_norm()>pow(10,-9)* dof_handler.n_dofs());
output_results (refinement_cycle);
// refine_mesh (refinement_cycle) ;
}
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Viscoelastic;
StokesProblem<2> flow_problem1(1);
flow_problem1.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;
}
