Dear All,
I have 2 questions when implementing of adaptive mesh refinement in
time-dependent problem, and the attached is my test code: test_for_ask.cc
1. I use the form to update my mesh:
[...]
timestep_number = 1;
setup_dofs ();
[...]
for (timestep_number=2; timestep_number<=100; ++timestep_number) {
[...]
refine_mesh (max_grid_level, min_grid_level);
setup_dofs ();
[...]
}
But starting from timestep_number=2, the setup_dof () cannot work. the
error is:
--------------------------------------------------------
An error occurred in line <128> of file
</home/yyyyy/boost/deal.ii-8.5.1/src/source/base/subscriptor.cc> in function
void dealii::Subscriptor::check_no_subscribers() const
The violated condition was:
counter == 0
Additional information:
Object of class N6dealii15SparsityPatternE is still used by 1 other
objects.
(Additional information:
from Subscriber SparseMatrix)
See the entry in the Frequently Asked Questions of deal.II (linked to from
http://www.dealii.org/) for a lot more information on what this error means
and how to fix programs in which it happens.
Stacktrace:
-----------
#0 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::Subscriptor::check_no_subscribers() const
#1 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::Subscriptor::~Subscriptor()
#2 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::SparsityPattern::~SparsityPattern()
#3 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::SparsityPattern::~SparsityPattern()
#4 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::BlockSparsityPatternBase<dealii::SparsityPattern>::reinit(unsigned
int, unsigned int)
#5 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::BlockSparsityPattern::copy_from(dealii::BlockDynamicSparsityPattern
const&)
#6 ./Problem13-2: Step22::StokesProblem<2>::setup_dofs()
#7 ./Problem13-2: Step22::StokesProblem<2>::run()
#8 ./Problem13-2: main
--------------------------------------------------------
make[3]: *** [CMakeFiles/run] Aborted (core dumped)
make[2]: *** [CMakeFiles/run.dir/all] Error 2
make[1]: *** [CMakeFiles/run.dir/rule] Error 2
make: *** [run] Error 2
And my setup_dofs () is written by learning from step-22, which means it
use BlockSparsityPattern:
[...]
system_matrix.clear ();
system_matrix_preconditioner.clear ();
{
BlockDynamicSparsityPattern dsp (2,2);
dsp.block(0,0).reinit (n_phi, n_phi);
dsp.block(1,0).reinit (n_e, n_phi);
dsp.block(0,1).reinit (n_phi, n_e);
dsp.block(1,1).reinit (n_e, n_e);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);
}
system_matrix.reinit (sparsity_pattern);
system_matrix_preconditioner.reinit (sparsity_pattern);
[...]
I had tried several methods all cannot work:(if I set another
BlockSparsityPattern sparsity_pattern_new, and use it in timestep_number =
2, it can go, but stop at timestep_number = 3; or I use
sparsity_pattern_pointer
= std_cxx11::shared_ptr<typename BSP_solution<dim>::type>(new
typename BSP_solution<dim>::type()); it also doesn't work, where:
template <int dim>
struct BSP_solution;
template <>
struct BSP_solution<2>
{
typedef BlockSparsityPattern type;
};
and now I don't know how to correct it. Can anyone give me a help?
2. My refine_mesh() is:
template <int dim>
void
StokesProblem<dim>::refine_mesh (const unsigned int max_grid_level,
const unsigned int min_grid_level)
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
FEValuesExtractors::Scalar phase(0);
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1>(degree+1),
typename FunctionMap<dim>::type(),
old_solution,
estimated_error_per_cell,
fe.component_mask(phase));
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.2,0.1);
if (triangulation.n_levels() > max_grid_level)
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(max_grid_level);
cell != triangulation.end(); ++cell)
cell->clear_refine_flag ();
/* this part for clear_coarsen_flag() cannot work */
if (triangulation.n_levels() < min_grid_level)
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(min_grid_level);
cell != triangulation.end(); ++cell)
cell->clear_coarsen_flag ();
triangulation.execute_coarsening_and_refinement ();
}
at first, if I set :
triangulation_active.refine_global (8)
then I use
refine_mesh (9,6);
the error is:
--------------------------------------------------------
An error occurred in line <11031> of file
</home/yyyyy/boost/deal.ii-8.5.1/src/source/grid/tria.cc> in function
dealii::Triangulation<dim, spacedim>::raw_quad_iterator
dealii::Triangulation<dim, spacedim>::begin_raw_quad(unsigned int) const
[with int dim = 2; int spacedim = 2; dealii::Triangulation<dim,
spacedim>::raw_quad_iterator =
dealii::TriaRawIterator<dealii::CellAccessor<2, 2> >]
The violated condition was:
level<n_global_levels() || level<levels.size()
Additional information:
The given level 6 is not in the valid range!
Stacktrace:
-----------
#0 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::Triangulation<2, 2>::begin_raw_quad(unsigned int) const
#1 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::Triangulation<2, 2>::begin_quad(unsigned int) const
#2 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::Triangulation<2, 2>::begin_active_quad(unsigned int) const
#3 //home/yyyyy/boost/deal.ii-8.5.1/deal.II/lib/libdeal_II.g.so.8.5.1:
dealii::Triangulation<2, 2>::begin_active(unsigned int) const
#4 ./Problem13-2: Step22::StokesProblem<2>::refine_mesh(unsigned int,
unsigned int)
#5 ./Problem13-2: Step22::StokesProblem<2>::run()
#6 ./Problem13-2: main
--------------------------------------------------------
make[3]: *** [CMakeFiles/run] Aborted (core dumped)
make[2]: *** [CMakeFiles/run.dir/all] Error 2
make[1]: *** [CMakeFiles/run.dir/rule] Error 2
make: *** [run] Error 2
I also don't understand what is the valid range above and how to correct it?
Thanks in advance!
Best,
Chucui
--
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/quadrature_lib.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/function.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/convergence_table.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/tensor_function.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/solver_gmres.h>
#include <deal.II/lac/sparse_ilu.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.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/lac/block_matrix_array.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.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_dgq.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_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/numerics/solution_transfer.h>
#include <math.h>
#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/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/constraint_matrix.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/tria_boundary_lib.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_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/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>
// Then we need to include the header file for the sparse direct solver
// UMFPACK:
#include <deal.II/lac/sparse_direct.h>
// This includes the library for the incomplete LU factorization that will be
// used as a preconditioner in 3D:
#include <deal.II/lac/sparse_ilu.h>
// This is C++:
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
using namespace std;
// As in all programs, the namespace dealii is included:
namespace Step22
{
using namespace dealii;
template <int dim>
class StokesProblem
{
public:
StokesProblem (const unsigned int phi_degree,
const unsigned int e_degree,
const unsigned int U_degree,
const unsigned int q_degree);
void run ();
BlockVector<double> project_gradient(const BlockVector<double> solution_single);
private:
void setup_dofs ();
void assemble_system_1 (const double time,
const Vector<double> phi_0,
const Vector<double> e_0,
const Vector<double> U_n,
const Vector<double> q_n,
const Vector<double> old_De,
const Vector<double> older_De,
const Vector<double> old_g,
const Vector<double> older_g,
const Vector<double> old_h,
const Vector<double> older_h,
const BlockVector<double> old_H,
const BlockVector<double> older_H);
void assemble_system_H0 (const Vector<double> phi_0);
void solve_1 ();
void solve_H ();
void refine_mesh (const unsigned int max_grid_level,
const unsigned int min_grid_level);
//void process_solution(const double time);
const unsigned int degree;
Triangulation<dim> triangulation_active;
Triangulation<dim> triangulation;
FESystem<dim> fe;
FESystem<dim> fe_H;
FE_Q<dim> fe_De;
FE_Q<dim> fe_tau;
FE_Q<dim> fe_T;
FE_Q<dim> fe_g;
FE_Q<dim> fe_h;
FE_Q<dim> fe_phi;
FE_Q<dim> fe_e;
FE_Q<dim> fe_U;
FE_Q<dim> fe_q;
DoFHandler<dim> dof_handler;
DoFHandler<dim> dof_handler_H;
DoFHandler<dim> dof_handler_De;
DoFHandler<dim> dof_handler_tau;
DoFHandler<dim> dof_handler_phi;
DoFHandler<dim> dof_handler_T;
DoFHandler<dim> dof_handler_g;
DoFHandler<dim> dof_handler_h;
DoFHandler<dim> dof_handler_e;
DoFHandler<dim> dof_handler_U;
DoFHandler<dim> dof_handler_q;
ConstraintMatrix constraints;
ConstraintMatrix constraints_H;
ConstraintMatrix constraints_De;
ConstraintMatrix constraints_tau;
ConstraintMatrix constraints_T;
ConstraintMatrix constraints_g;
ConstraintMatrix constraints_h;
ConstraintMatrix constraints_phi;
ConstraintMatrix constraints_e;
ConstraintMatrix constraints_U;
ConstraintMatrix constraints_q;
BlockSparsityPattern sparsity_pattern;
BlockSparsityPattern sparsity_pattern_H;
SparsityPattern sparsity_pattern_U;
SparsityPattern sparsity_pattern_q;
SparsityPattern sparsity_pattern_De;
SparsityPattern sparsity_pattern_tau;
SparsityPattern sparsity_pattern_T;
SparsityPattern sparsity_pattern_g;
SparsityPattern sparsity_pattern_h;
BlockSparseMatrix<double> system_matrix;
BlockSparseMatrix<double> system_matrix_preconditioner;
BlockSparseMatrix<double> entropy_matrix;
BlockSparseMatrix<double> system_matrix_H;
BlockSparseMatrix<double> system_matrix_H_preconditioner;
BlockVector<double> solution_1;
BlockVector<double> solution;
BlockVector<double> old_solution;
BlockVector<double> older_solution;
BlockVector<double> older_H, old_H, new_H, system_rhs_H;
BlockVector<double> system_rhs;
BlockVector<double> energy_rhs;
SparseMatrix<double> system_matrix_U;
SparseMatrix<double> system_matrix_q;
SparseMatrix<double> entropy_matrix_U;
SparseMatrix<double> entropy_matrix_q;
SparseMatrix<double> system_matrix_De;
SparseMatrix<double> system_matrix_tau;
SparseMatrix<double> system_matrix_T;
SparseMatrix<double> system_matrix_g;
SparseMatrix<double> system_matrix_h;
Vector<double> phi_0, e_0, U_0, q_0, T_0, g_0, h_0, tau_0, De_0;
Vector<double> system_rhs_De, older_De, old_De, new_De;
Vector<double> system_rhs_U, new_U, older_U, old_U, U_n;
Vector<double> system_rhs_q, new_q, older_q, old_q, q_n;
Vector<double> system_rhs_T, new_T, old_T, older_T;
Vector<double> system_rhs_g, new_g, old_g, older_g;
Vector<double> system_rhs_h, new_h, old_h, older_h;
Vector<double> system_rhs_tau, new_tau, old_tau, older_tau;
double time_step, time;
const double A, B, sita_0, eps, m, eps_4, tau, k_0, C_L, D_u, L_1, L_2, L_0, T_M, C_0, C_1, eps_F, e_L_T_M;//, s;
unsigned int timestep_number = 1;
ConvergenceTable convergence_table;
};
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int phi_degree,
const unsigned int e_degree,
const unsigned int U_degree,
const unsigned int q_degree)
:
degree (phi_degree),
fe (FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), 1),
fe_H (FE_Q<dim>(degree), 1,
FE_Q<dim>(degree), 1),
fe_De (degree),
fe_tau (degree),
fe_T (degree),
fe_g (degree),
fe_h (degree),
fe_phi (degree),
fe_e (degree),
fe_U (degree),
fe_q (degree),
dof_handler (triangulation),
dof_handler_H (triangulation),
dof_handler_De (triangulation),
dof_handler_tau (triangulation),
dof_handler_T (triangulation),
dof_handler_g (triangulation),
dof_handler_h (triangulation),
dof_handler_phi (triangulation),
dof_handler_e (triangulation),
dof_handler_U (triangulation),
dof_handler_q (triangulation),
time_step (0.01),
timestep_number (1),
A (100.0),
B (1.0),
sita_0 (0*numbers::PI/4),
eps (0.01),
m (4.0),
eps_4 (1./15.0),
tau (3.0),
k_0 (0.0),
C_L (0.6),
D_u (0.0001),
L_1 (0.0),
L_2 (0.0),
L_0 (0.5),
T_M (1.0),
C_0 (1.0),
C_1 (1.0),
eps_F (1.0),
e_L_T_M (0.0)
{}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_phi : public Function<dim>
{
public:
InitialValues_phi () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_phi<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
const double C_0 = 1.0,
C_1 = 1.0,
A = 100.0;
const double delta = 0.01, r0 = std::sqrt(0.018);
double phi_value = 0, e_value = 0, q_value = 0, s_value = 0, r = 0, T_value = 0;
//r = std::sqrt((p[0] - 9.0) * (p[0] - 9.0) + (p[1] - 9.0) * (p[1] - 9.0)) ;
// r = std::sqrt((p[0] - 1.125) * (p[0] - 1.125) + (p[1] - 1.125) * (p[1] - 1.125)) ;
r = std::sqrt((p[0] - 2.25) * (p[0] - 2.25) + (p[1] - 2.25) * (p[1] - 2.25)) ;
phi_value = -0.5 * std::tanh((r0-r)/delta) + 0.5;
//phi_value = 0.5 * std::cos(p[0]) * std::cos(p[1]) +0.5;
return phi_value;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_e : public Function<dim>
{
public:
InitialValues_e () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_e<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
const double e_L_T_M = 0.0,
C_L = 0.6,
C_0 = 1.0,
C_1 = 1.0,
A = 100.0,
T_M = 1.0,
L_0 = 0.5,
L_1 = 0.0,
L_2 = 0.0;
const double delta = 0.01, r0 = std::sqrt(0.018);
double phi_value = 0, e_value = 0, r = 0, T_value = 0;
//r = std::sqrt((p[0] - 9.0) * (p[0] - 9.0) + (p[1] - 9.0) * (p[1] - 9.0)) ;
//r = std::sqrt((p[0] - 1.125) * (p[0] - 1.125) + (p[1] - 1.125) * (p[1] - 1.125)) ;
r = std::sqrt((p[0] - 2.25) * (p[0] - 2.25) + (p[1] - 2.25) * (p[1] - 2.25)) ;
phi_value = -0.5 * std::tanh((r0-r)/delta) + 0.5;
double P_phi = 0;
P_phi = 30 * (0.2 * std::pow(phi_value,5) - 0.5 * std::pow(phi_value,4) + (1.0/3.0) * std::pow(phi_value,3));
T_value = 0.5 ;
e_value = e_L_T_M + C_L * (T_value-T_M) + (P_phi-1)*(L_0 + L_1*(T_value-T_M) + L_2*(T_value-T_M)*(T_value-T_M));
return e_value;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_q : public Function<dim>
{
public:
InitialValues_q () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_q<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
const double e_L_T_M = 0.0,
C_L = 0.6,
C_0 = 1.0,
C_1 = 1.0,
A = 100.0,
T_M = 1.0,
L_0 = 0.5,
L_1 = 0.0,
L_2 = 0.0,
eps_F = 1.0;
const double delta = 0.01, r0 = std::sqrt(0.018);
double phi_value = 0, e_value = 0, q_value = 0, s_value = 0, r = 0, T_value = 0;
//r = std::sqrt((p[0] - 9.0) * (p[0] - 9.0) + (p[1] - 9.0) * (p[1] - 9.0)) ;
//r = std::sqrt((p[0] - 1.125) * (p[0] - 1.125) + (p[1] - 1.125) * (p[1] - 1.125)) ;
r = std::sqrt((p[0] - 2.25) * (p[0] - 2.25) + (p[1] - 2.25) * (p[1] - 2.25)) ;
phi_value = -0.5 * std::tanh((r0-r)/delta) + 0.5;
double P_phi = 0, Q_T= 0, F_phi = 0;
P_phi = 30 * (0.2 * std::pow(phi_value,5) - 0.5 * std::pow(phi_value,4) + (1.0/3.0) * std::pow(phi_value,3));
T_value = 0.5;
e_value = e_L_T_M + C_L * (T_value-T_M) + (P_phi-1)*(L_0 + L_1*(T_value-T_M) + L_2*(T_value-T_M)*(T_value-T_M));
Q_T = L_0 * ((1./T_M) - (1./T_value));
F_phi = 0.25 * phi_value * phi_value * (1-phi_value) * (1-phi_value);
s_value = (e_value/T_value) + C_L * (std::log(T_value) - std::log(T_M)) + (e_L_T_M - C_L * T_M) * ((1./T_M) - (1./T_value)) + (P_phi - 1) * Q_T - F_phi;
q_value = std::sqrt(2 * (-s_value - 0.5 * C_0 * phi_value * phi_value - 0.5 * C_1 * e_value * e_value + A));
return q_value;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_U : public Function<dim>
{
public:
InitialValues_U () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_U<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
const double C_0 = 1.0,
C_1 = 1.0,
A = 100.0,
k0 = 0.0,
B = 1.0,
eps_4 = (1./15.0),
sita_0 = 0 * numbers::PI/4,
m = 4.0,
eps = 0.01;
const double delta = 0.01, r0 = std::sqrt(0.018);
double U_value = 0, e_value = 0, kappa_value = 0, s_value = 0, r = 0, T_value = 0;
//r = std::sqrt((p[0] - 9.0) * (p[0] - 9.0) + (p[1] - 9.0) * (p[1] - 9.0)) ;
//r = std::sqrt((p[0] - 1.125) * (p[0] - 1.125) + (p[1] - 1.125) * (p[1] - 1.125)) ;
r = std::sqrt((p[0] - 2.25) * (p[0] - 2.25) + (p[1] - 2.25) * (p[1] - 2.25)) ;
Tensor<1,dim> grad_phi_value;
grad_phi_value[0] = ((1-std::tanh((r0-r)/delta)*std::tanh((r0-r)/delta)) * (p[0]-2.25)) / (4*delta*r);
grad_phi_value[1] = ((1-std::tanh((r0-r)/delta)*std::tanh((r0-r)/delta)) * (p[1]-2.25)) / (4*delta*r);
kappa_value = 1 + eps_4 * std::cos(m * (std::atan2(grad_phi_value[1], grad_phi_value[0]) - sita_0));
U_value = std::sqrt(eps*eps*(kappa_value*kappa_value-k0)*(grad_phi_value[0]*grad_phi_value[0]+grad_phi_value[1]*grad_phi_value[1])+2*B);
//phi_value = 0.5 * std::cos(p[0]) * std::cos(p[1]) +0.5;
return U_value;
}
double gFunction (const double phi_value,
const double e_value)
{
double g_value = 0;
const double e_L_T_M = 0.0,
C_L = 0.6,
C_0 = 1.0,
C_1 = 1.0,
A = 100.0,
T_M = 1.0,
L_0 = 0.5,
L_1 = 0.0,
L_2 = 0.0;
double s_value = 0, T_value = 0, qq_value = 0;
double P_phi = 0, Q_T = 0, F_phi = 0, p_partial_phi = 0, F_partial_phi = 0;
P_phi = 30 * (0.2 * std::pow(phi_value,5) - 0.5 * std::pow(phi_value,4) + (1.0/3.0) * std::pow(phi_value,3));
T_value = (e_value - e_L_T_M - (P_phi - 1)*L_0)/C_L + T_M ;
Q_T = L_0 * ((1./T_M) - (1./T_value ));
F_phi = 0.25 * phi_value * phi_value * (1-phi_value) * (1-phi_value);
p_partial_phi = 30 * phi_value * phi_value * (1-phi_value) * (1-phi_value);
F_partial_phi = 0.5 * phi_value * (phi_value - 1) * (2 * phi_value - 1);
s_value = (e_value /T_value ) + C_L * (std::log(T_value ) - std::log(T_M)) + (e_L_T_M - C_L * T_M) * ((1./T_M) - (1./T_value )) + (P_phi - 1) * Q_T - F_phi ;
qq_value = std::sqrt(2 * (-s_value - 0.5 * C_0 * phi_value * phi_value - 0.5 * C_1 * e_value * e_value + A));
g_value =(-p_partial_phi * Q_T + F_partial_phi - C_0 * phi_value)/(qq_value);
return g_value;
}
double hFunction (const double phi_value,
const double e_value)
{
double h_value = 0;
const double e_L_T_M = 0.0,
C_L = 0.6,
C_0 = 1.0,
C_1 = 1.0,
A = 100.0,
T_M = 1.0,
L_0 = 0.5,
L_1 = 0.0,
L_2 = 0.0;
double s_value = 0, T_value = 0, qq_value = 0;
double P_phi = 0, Q_T= 0, F_phi = 0;
P_phi = 30 * (0.2 * std::pow(phi_value,5) - 0.5 * std::pow(phi_value,4) + (1.0/3.0) * std::pow(phi_value,3));
T_value = (e_value - e_L_T_M - (P_phi-1)*L_0)/C_L + T_M ;
s_value = (e_value/T_value) + C_L * (std::log(T_value) - std::log(T_M)) + (e_L_T_M - C_L * T_M) * ((1./T_M) - (1./T_value)) + (P_phi - 1) * Q_T - F_phi;
qq_value = std::sqrt(2 * (-s_value - 0.5 * C_0 * phi_value * phi_value - 0.5 * C_1 * e_value * e_value + A));
h_value = ((-1) * (1./T_value) - C_1 * e_value) / (qq_value);
return h_value;
}
template <int dim>
double tauFunction (const Tensor<1,dim> grad_phi)
{
double tau1 = 0, tau2 = 0, n_x_1 = 0, n_x_2 = 0, n_y_1 = 0, n_y_2 = 0;
const double eps_4 = (1./15.0),
m = 4.0,
protect_eps_large = 1e-216;
return 3.0;
}
template <int dim>
Tensor<1,dim> HFunction (const Tensor<1,dim> grad_phi)
{
//Tensor<1,dim> return_value;
double kappa_1 = 0, kappa_partial_1 = 0, UU_1 = 0;
const double eps = 0.01,
eps_4 = (1./15.0),//(1./(m^2-1))
m = 4.0,
B = 1,
k_0 = 0.0,
sita_0 = 0 * numbers::PI/4,
protect_eps = 1e-216;
kappa_1 = 1 + eps_4 * std::cos(m * (std::atan2(grad_phi[1], grad_phi[0]) - sita_0));
kappa_partial_1 = -1 * m * eps_4 * std::sin(m * (std::atan2(grad_phi[1], grad_phi[0]) - sita_0));
UU_1 = std::sqrt(eps * eps * (kappa_1 * kappa_1 - k_0) * (grad_phi[0] * grad_phi[0] + grad_phi[1] * grad_phi[1]) + 2 * B);
Tensor<1,dim> H_value;
H_value[0] = ((eps * eps) / UU_1) * (kappa_1 * kappa_partial_1 * (-1) * grad_phi[1] + (kappa_1 * kappa_1 - k_0) * grad_phi[0]);
H_value[1] = ((eps * eps) / UU_1) * (kappa_1 * kappa_partial_1 * grad_phi[0] + (kappa_1 * kappa_1 - k_0) * grad_phi[1]);
return H_value;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_T : public Function<dim>
{
public:
InitialValues_T () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_T<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
double T_value = 0;
T_value = 0.5;
return T_value;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_tau : public Function<dim>
{
public:
InitialValues_tau () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_tau<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
double tau_value = 0;
tau_value = 3.0;
return tau_value;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_De : public Function<dim>
{
public:
InitialValues_De () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_De<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
double De_value = 0;
const double delta = 0.01, r0 = std::sqrt(0.018), D_u = 0.0001;
double phi_value = 0, r = 0, T_value = 0.5;
//r = std::sqrt((p[0] - 9.0) * (p[0] - 9.0) + (p[1] - 9.0) * (p[1] - 9.0)) ;
// r = std::sqrt((p[0] - 1.125) * (p[0] - 1.125) + (p[1] - 1.125) * (p[1] - 1.125)) ;
r = std::sqrt((p[0] - 2.25) * (p[0] - 2.25) + (p[1] - 2.25) * (p[1] - 2.25)) ;
phi_value = -0.5 * std::tanh((r0-r)/delta) + 0.5;
De_value = D_u * T_value * T_value;
return De_value;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_g : public Function<dim>
{
public:
InitialValues_g () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_g<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
double g_value = 0;
const double delta = 0.01, r0 = std::sqrt(0.018);
double phi_value = 0, T_value = 0, e_value = 0, r = 0;
const double e_L_T_M = 0.0,
C_L = 0.6,
C_0 = 1.0,
C_1 = 1.0,
A = 100.0,
T_M = 1.0,
L_0 = 0.5,
L_1 = 0.0,
L_2 = 0.0;
double s_value = 0, qq_value = 0;
double P_phi = 0, Q_T = 0, F_phi = 0, p_partial_phi = 0, F_partial_phi = 0;
r = std::sqrt((p[0] - 2.25) * (p[0] - 2.25) + (p[1] - 2.25) * (p[1] - 2.25)) ;
//r = std::sqrt((p[0] - 9.0) * (p[0] - 9.0) + (p[1] - 9.0) * (p[1] - 9.0)) ;
// r = std::sqrt((p[0] - 1.125) * (p[0] - 1.125) + (p[1] - 1.125) * (p[1] - 1.125)) ;
T_value = 0.5;
phi_value = -0.5 * std::tanh((r0-r)/delta) + 0.5;
P_phi = 30 * (0.2 * std::pow(phi_value,5) - 0.5 * std::pow(phi_value,4) + (1.0/3.0) * std::pow(phi_value,3));
Q_T = L_0 * ((1./T_M) - (1./T_value ));
F_phi = 0.25 * phi_value * phi_value * (1-phi_value) * (1-phi_value);
p_partial_phi = 30 * phi_value * phi_value * (1-phi_value) * (1-phi_value);
F_partial_phi = 0.5 * phi_value * (phi_value - 1) * (2 * phi_value - 1);
s_value = (e_value /T_value ) + C_L * (std::log(T_value ) - std::log(T_M)) + (e_L_T_M - C_L * T_M) * ((1./T_M) - (1./T_value )) + (P_phi - 1) * Q_T - F_phi ;
qq_value = std::sqrt(2 * (-s_value - 0.5 * C_0 * phi_value * phi_value - 0.5 * C_1 * e_value * e_value + A));
g_value =(-p_partial_phi * Q_T + F_partial_phi - C_0 * phi_value)/(qq_value);
return g_value;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
class InitialValues_h : public Function<dim>
{
public:
InitialValues_h () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double InitialValues_h<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
double h_value = 0;
const double delta = 0.01, r0 = std::sqrt(0.018);
double phi_value = 0, T_value = 0, e_value = 0, r = 0;
const double e_L_T_M = 0.0,
C_L = 0.6,
C_0 = 1.0,
C_1 = 1.0,
A = 100.0,
T_M = 1.0,
L_0 = 0.5,
L_1 = 0.0,
L_2 = 0.0;
double s_value = 0, qq_value = 0;
double P_phi = 0, Q_T= 0, F_phi = 0;
r = std::sqrt((p[0] - 2.25) * (p[0] - 2.25) + (p[1] - 2.25) * (p[1] - 2.25)) ;
//r = std::sqrt((p[0] - 9.0) * (p[0] - 9.0) + (p[1] - 9.0) * (p[1] - 9.0)) ;
// r = std::sqrt((p[0] - 1.125) * (p[0] - 1.125) + (p[1] - 1.125) * (p[1] - 1.125)) ;
T_value = 0.5;
phi_value = -0.5 * std::tanh((r0-r)/delta) + 0.5;
P_phi = 30 * (0.2 * std::pow(phi_value,5) - 0.5 * std::pow(phi_value,4) + (1.0/3.0) * std::pow(phi_value,3));
T_value = (e_value - e_L_T_M - (P_phi-1)*L_0)/C_L + T_M ;
s_value = (e_value/T_value) + C_L * (std::log(T_value) - std::log(T_M)) + (e_L_T_M - C_L * T_M) * ((1./T_M) - (1./T_value)) + (P_phi - 1) * Q_T - F_phi;
qq_value = std::sqrt(2 * (-s_value - 0.5 * C_0 * phi_value * phi_value - 0.5 * C_1 * e_value * e_value + A));
h_value = ((-1) * (1./T_value) - C_1 * e_value) / (qq_value);
return h_value;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <int dim>
void StokesProblem<dim>::setup_dofs ()
{
std::vector<unsigned int> block_component (2);
block_component[0] = 0;
block_component[1] = 1;
dof_handler.distribute_dofs (fe);
//DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler, block_component);
constraints.clear ();
const ComponentMask component_mask = ComponentMask();
DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints, component_mask);
DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints, component_mask);
constraints.close ();
std::vector<types::global_dof_index> dofs_per_block (2);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
const unsigned int n_phi = dofs_per_block[0],
n_e = dofs_per_block[1];
//std::cout << "dof 0-7 "<< std::endl;
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< " (" << n_phi << '+' << n_e << ')'
<< std::endl;
system_matrix.clear ();
system_matrix_preconditioner.clear ();
{
BlockDynamicSparsityPattern dsp (2,2);
dsp.block(0,0).reinit (n_phi, n_phi);
dsp.block(1,0).reinit (n_e, n_phi);
dsp.block(0,1).reinit (n_phi, n_e);
dsp.block(1,1).reinit (n_e, n_e);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
std::cout << "make_sparsity_pattern " << std::endl;
sparsity_pattern.copy_from (dsp);
std::cout << "sparsity_pattern_copy " << std::endl;
}
system_matrix.reinit (sparsity_pattern);
system_matrix_preconditioner.reinit (sparsity_pattern);
entropy_matrix.clear ();
entropy_matrix.reinit (sparsity_pattern);
solution.reinit (2);
solution.block(0).reinit (n_phi);
solution.block(1).reinit (n_e);
solution.collect_sizes ();
solution_1.reinit (2);
solution_1.block(0).reinit (n_phi);
solution_1.block(1).reinit (n_e);
solution_1.collect_sizes ();
old_solution.reinit (2);
old_solution.block(0).reinit (n_phi);
old_solution.block(1).reinit (n_e);
old_solution.collect_sizes ();
older_solution.reinit (2);
older_solution.block(0).reinit (n_phi);
older_solution.block(1).reinit (n_e);
older_solution.collect_sizes ();
system_rhs.reinit (2);
system_rhs.block(0).reinit (n_phi);
system_rhs.block(1).reinit (n_e);
system_rhs.collect_sizes ();
energy_rhs.reinit (2);
energy_rhs.block(0).reinit (n_phi);
energy_rhs.block(1).reinit (n_e);
energy_rhs.collect_sizes ();
///////////////////////////////////////////////////////////////////////////////////////////////////////
std::vector<unsigned int> block_component_H (2);
block_component_H[0] = 0;
block_component_H[1] = 1;
dof_handler_H.distribute_dofs (fe_H);
//DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler_H, block_component_H);
constraints_H.clear ();
//const ComponentMask component_mask_H = ComponentMask();
//DoFTools::make_periodicity_constraints(dof_handler_H,0,1,0, constraints_H, component_mask_H);
//DoFTools::make_periodicity_constraints(dof_handler_H,2,3,1, constraints_H, component_mask_H);
constraints_H.close ();
std::vector<types::global_dof_index> dofs_per_block_H (2);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block_H, block_component_H);
const unsigned int n_H_x = dofs_per_block_H[0],
n_H_y = dofs_per_block_H[1];
system_matrix_H.clear ();
system_matrix_H_preconditioner.clear ();
{
BlockDynamicSparsityPattern dsp_H (2,2);
dsp_H.block(0,0).reinit (n_H_x, n_H_x);
dsp_H.block(1,0).reinit (n_H_y, n_H_x);
dsp_H.block(0,1).reinit (n_H_x, n_H_y);
dsp_H.block(1,1).reinit (n_H_y, n_H_y);
dsp_H.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler_H, dsp_H, constraints_H, false);
sparsity_pattern_H.copy_from (dsp_H);
}
system_matrix_H.reinit (sparsity_pattern_H);
system_matrix_H_preconditioner.reinit (sparsity_pattern_H);
system_rhs_H.reinit (2);
system_rhs_H.block(0).reinit (n_H_x);
system_rhs_H.block(1).reinit (n_H_y);
system_rhs_H.collect_sizes ();
older_H.reinit (2);
older_H.block(0).reinit (n_H_x);
older_H.block(1).reinit (n_H_y);
older_H.collect_sizes ();
old_H.reinit (2);
old_H.block(0).reinit (n_H_x);
old_H.block(1).reinit (n_H_y);
old_H.collect_sizes ();
new_H.reinit (2);
new_H.block(0).reinit (n_H_x);
new_H.block(1).reinit (n_H_y);
new_H.collect_sizes ();
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////
dof_handler_phi.distribute_dofs (fe_phi);
dof_handler_e.distribute_dofs (fe_e);
DoFRenumbering::component_wise (dof_handler_phi);
DoFRenumbering::component_wise (dof_handler_e);
constraints_phi.close ();
constraints_e.close ();
DynamicSparsityPattern dsp_phi (n_phi, n_phi);
DoFTools::make_sparsity_pattern (dof_handler_phi, dsp_phi, constraints_phi, false);
DynamicSparsityPattern dsp_e (n_e, n_e);
DoFTools::make_sparsity_pattern (dof_handler_e, dsp_e, constraints_e, false);
dof_handler_U.distribute_dofs (fe_U);
dof_handler_q.distribute_dofs (fe_q);
DoFRenumbering::component_wise (dof_handler_U);
const unsigned int n_U = dof_handler_U.n_dofs();
const unsigned int n_q = dof_handler_q.n_dofs();
// constraints_U.clear ();
//const ComponentMask component_mask_U = ComponentMask();
//DoFTools::make_periodicity_constraints(dof_handler_U,0,1,0, constraints_U, component_mask_U);
//DoFTools::make_periodicity_constraints(dof_handler_U,2,3,1, constraints_U, component_mask_U);
constraints_U.close ();
DoFRenumbering::component_wise (dof_handler_q);
// constraints_q.clear ();
//const ComponentMask component_mask_q = ComponentMask();
//DoFTools::make_periodicity_constraints(dof_handler_q,0,1,0, constraints_q, component_mask_q);
//DoFTools::make_periodicity_constraints(dof_handler_q,2,3,1, constraints_q, component_mask_q);
constraints_q.close ();
DynamicSparsityPattern dsp_U (n_U, n_U);
DoFTools::make_sparsity_pattern (dof_handler_U, dsp_U, constraints_U, false);
sparsity_pattern_U.copy_from (dsp_U);
system_matrix_U.reinit (sparsity_pattern_U);
entropy_matrix_U.reinit (sparsity_pattern_U);
DynamicSparsityPattern dsp_q (n_q, n_q);
DoFTools::make_sparsity_pattern (dof_handler_q, dsp_q, constraints_q, false);
sparsity_pattern_q.copy_from (dsp_q);
system_matrix_q.reinit (sparsity_pattern_q);
entropy_matrix_q.reinit (sparsity_pattern_q);
U_0.reinit (n_U);
phi_0.reinit (n_phi);
e_0.reinit (n_e);
q_0.reinit (n_q);
new_U.reinit (n_U);
old_U.reinit (n_U);
older_U.reinit (n_U);
U_n.reinit (n_U);
new_q.reinit (n_q);
old_q.reinit (n_q);
older_q.reinit (n_q);
q_n.reinit (n_q);
system_rhs_U.reinit (n_U);
system_rhs_q.reinit (n_q);
dof_handler_T.distribute_dofs (fe_T);
DoFRenumbering::component_wise (dof_handler_T);
const unsigned int n_T = dof_handler_T.n_dofs();
constraints_T.close ();
DynamicSparsityPattern dsp_T (n_T, n_T);
DoFTools::make_sparsity_pattern (dof_handler_T, dsp_T, constraints_T, false);
sparsity_pattern_T.copy_from (dsp_T);
system_matrix_T.reinit (sparsity_pattern_T);
new_T.reinit (n_T);
old_T.reinit (n_T);
T_0.reinit (n_T);
older_T.reinit (n_T);
system_rhs_T.reinit (n_T);
///////////////////////////////////////////////////////////////////////////////////////////////////////////
dof_handler_De.distribute_dofs (fe_De);
DoFRenumbering::component_wise (dof_handler_De);
const unsigned int n_De = dof_handler_De.n_dofs();
constraints_De.close ();
DynamicSparsityPattern dsp_De (n_De, n_De);
DoFTools::make_sparsity_pattern (dof_handler_De, dsp_De, constraints_De, false);
sparsity_pattern_De.copy_from (dsp_De);
system_matrix_De.reinit (sparsity_pattern_De);
De_0.reinit (n_De);
new_De.reinit (n_De);
old_De.reinit (n_De);
older_De.reinit (n_De);
system_rhs_De.reinit (n_De);
///////////////////////////////////////////////////////////////////////////////////////////////////////////
dof_handler_tau.distribute_dofs (fe_tau);
DoFRenumbering::component_wise (dof_handler_tau);
const unsigned int n_tau = dof_handler_tau.n_dofs();
constraints_tau.close ();
DynamicSparsityPattern dsp_tau (n_tau, n_tau);
DoFTools::make_sparsity_pattern (dof_handler_tau, dsp_tau, constraints_tau, false);
sparsity_pattern_tau.copy_from (dsp_tau);
system_matrix_tau.reinit (sparsity_pattern_tau);
tau_0.reinit (n_tau);
new_tau.reinit (n_tau);
old_tau.reinit (n_tau);
older_tau.reinit (n_tau);
system_rhs_tau.reinit (n_tau);
///////////////////////////////////////////////////////////////////////////////////////////////////////////
dof_handler_g.distribute_dofs (fe_g);
DoFRenumbering::component_wise (dof_handler_g);
const unsigned int n_g = dof_handler_g.n_dofs();
constraints_g.close ();
DynamicSparsityPattern dsp_g (n_g, n_g);
DoFTools::make_sparsity_pattern (dof_handler_g, dsp_g, constraints_g, false);
sparsity_pattern_g.copy_from (dsp_g);
system_matrix_g.reinit (sparsity_pattern_g);
g_0.reinit (n_g);
new_g.reinit (n_g);
old_g.reinit (n_g);
older_g.reinit (n_g);
system_rhs_g.reinit (n_g);
///////////////////////////////////////////////////////////////////////////////////////////////////////////
dof_handler_h.distribute_dofs (fe_h);
DoFRenumbering::component_wise (dof_handler_h);
const unsigned int n_h = dof_handler_h.n_dofs();
constraints_h.close ();
DynamicSparsityPattern dsp_h (n_h, n_h);
DoFTools::make_sparsity_pattern (dof_handler_h, dsp_h, constraints_h, false);
sparsity_pattern_h.copy_from (dsp_h);
system_matrix_h.reinit (sparsity_pattern_h);
h_0.reinit (n_h);
new_h.reinit (n_h);
old_h.reinit (n_h);
older_h.reinit (n_h);
system_rhs_h.reinit (n_h);
}
template <int dim>
void StokesProblem<dim>::assemble_system_1 (const double time,
const Vector<double> phi_0,
const Vector<double> e_0,
const Vector<double> U_n,
const Vector<double> q_n,
const Vector<double> old_De,
const Vector<double> older_De,
const Vector<double> old_g,
const Vector<double> older_g,
const Vector<double> old_h,
const Vector<double> older_h,
const BlockVector<double> old_H,
const BlockVector<double> older_H)
{
system_matrix=0;
system_rhs=0;
QGauss<dim> quadrature_formula(degree+2);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_phi_values (fe_phi, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_e_values (fe_e, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_U_values (fe_U, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_q_values (fe_q, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_H_values (fe_H, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_g_values (fe_g, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_h_values (fe_h, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_De_values (fe_De, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell);
FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<Vector<double> > old_H_values(n_q_points, Vector<double>(2));
std::vector<Vector<double> > older_H_values(n_q_points, Vector<double>(2));
std::vector<double> phi_0_values(n_q_points);
std::vector<double> e_0_values(n_q_points);
std::vector<double> old_U_values(n_q_points);
std::vector<double> old_q_values(n_q_points);
std::vector<double> old_g_values(n_q_points);
std::vector<double> older_g_values(n_q_points);
std::vector<double> old_h_values(n_q_points);
std::vector<double> older_h_values(n_q_points);
std::vector<double> old_De_values(n_q_points);
std::vector<double> older_De_values(n_q_points);
std::vector<Tensor<1, dim> > phi_0_gradients(n_q_points);
std::vector<Tensor<1, dim> > e_0_gradients(n_q_points);
std::vector<Tensor<1, dim> > old_q_gradients(n_q_points);
std::vector<Tensor<1, dim> > old_g_gradients(n_q_points);
std::vector<Tensor<1, dim> > older_g_gradients(n_q_points);
std::vector<Tensor<1, dim> > old_h_gradients(n_q_points);
std::vector<Tensor<1, dim> > older_h_gradients(n_q_points);
const FEValuesExtractors::Scalar phase (0);
const FEValuesExtractors::Scalar energy (1);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
typename DoFHandler<dim>::active_cell_iterator
cell_phi = dof_handler_phi.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_e = dof_handler_e.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_U = dof_handler_U.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_q = dof_handler_q.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_De = dof_handler_De.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_g = dof_handler_g.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_h = dof_handler_h.begin_active();
typename DoFHandler<dim>::active_cell_iterator
cell_H = dof_handler_H.begin_active();
for (; cell!=endc; ++cell, ++cell_phi, ++cell_e, ++cell_U, ++cell_q, ++cell_De, ++cell_g, ++cell_h, ++cell_H)
{
local_rhs = 0;
local_matrix=0;
fe_values.reinit (cell);
fe_phi_values.reinit (cell_phi);
fe_e_values.reinit (cell_e);
fe_U_values.reinit (cell_U);
fe_q_values.reinit (cell_q);
fe_De_values.reinit (cell_De);
fe_g_values.reinit (cell_g);
fe_h_values.reinit (cell_h);
fe_H_values.reinit (cell_H);
fe_phi_values.get_function_values (phi_0, phi_0_values);
fe_e_values.get_function_values (e_0, e_0_values);
fe_H_values.get_function_values (old_H, old_H_values);
fe_H_values.get_function_values (older_H, older_H_values);
fe_phi_values.get_function_gradients(phi_0, phi_0_gradients);
fe_e_values.get_function_gradients(e_0, e_0_gradients);
fe_U_values.get_function_values (U_n, old_U_values);
fe_q_values.get_function_values (q_n, old_q_values);
fe_g_values.get_function_values (old_g, old_g_values);
fe_g_values.get_function_values (older_g, older_g_values);
fe_h_values.get_function_values (old_h, old_h_values);
fe_h_values.get_function_values (older_h, older_h_values);
fe_De_values.get_function_values (old_De, old_De_values);
fe_De_values.get_function_values (older_De, older_De_values);
fe_q_values.get_function_gradients(q_n, old_q_gradients);
fe_g_values.get_function_gradients(old_g, old_g_gradients);
fe_g_values.get_function_gradients(older_g, older_g_gradients);
fe_h_values.get_function_gradients(old_h, old_h_gradients);
fe_h_values.get_function_gradients(older_h, older_h_gradients);
for (unsigned int q=0; q<n_q_points; ++q)
{
const double old_phi = phi_0_values[q];
const double old_e = e_0_values[q];
const double old_U = old_U_values[q];
const double old_q = old_q_values[q];
// double tau_Bar = 1.5 * old_tau_values[q] - 0.5 * older_tau_values[q];
double De_Bar = 1.5 * old_De_values[q] - 0.5 * older_De_values[q];
double g_Bar = 1.5 * old_g_values[q] - 0.5 * older_g_values[q];
double h_Bar = 1.5 * old_h_values[q] - 0.5 * older_h_values[q];
Tensor<1,dim> grad_old_phi;
Tensor<1,dim> grad_old_e;
Tensor<1,dim> grad_old_q;
Tensor<1,dim> grad_g_Bar;
Tensor<1,dim> grad_h_Bar;
Tensor<1,dim> H_Bar;
for (unsigned int d=0; d<dim; ++d)
{
grad_old_phi[d] = phi_0_gradients[q][d];
grad_old_e[d] = e_0_gradients[q][d];
grad_old_q[d] = old_q_gradients[q][d];
grad_g_Bar[d] = 1.5 * old_g_gradients[q][d] - 0.5 * older_g_gradients[q][d] ;
grad_h_Bar[d] = 1.5 * old_h_gradients[q][d] - 0.5 * older_h_gradients[q][d] ;
H_Bar[d] = 1.5 * old_H_values[q](d) - 0.5 * older_H_values[q](d) ;
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> grad_psi_i_phi = fe_values[phase].gradient (i, q);
const Tensor<1,dim> grad_psi_i_e = fe_values[energy].gradient (i, q);
const double psi_i_phi = fe_values[phase].value (i, q);
const double psi_i_e = fe_values[energy].value (i, q);
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const Tensor<1,dim> grad_psi_j_phi = fe_values[phase].gradient (j, q);
const Tensor<1,dim> grad_psi_j_e = fe_values[energy].gradient (j, q);
const double psi_j_phi = fe_values[phase].value (j, q);
const double psi_j_e = fe_values[energy].value (j, q);
local_matrix(i,j) += ((C_0 / 2.0) * psi_i_phi * psi_j_phi
+ (1 / time_step) * 3 * psi_i_phi * psi_j_phi
+ (eps * eps * k_0) * 0.5 * grad_psi_i_phi * grad_psi_j_phi
+ 0.5 * g_Bar * h_Bar * psi_i_phi * psi_j_e
+ 0.5 * g_Bar * g_Bar * psi_i_phi * psi_j_phi
+ 0.5 * 1 * (H_Bar * grad_psi_j_phi) * (H_Bar * grad_psi_i_phi)
+ (1./ time_step) * psi_i_e * psi_j_e
+ 0.5 * C_1 * De_Bar * grad_psi_i_e * grad_psi_j_e
+ 0.5 * De_Bar * h_Bar * g_Bar * grad_psi_i_e * grad_psi_j_phi
+ 0.5 * De_Bar * h_Bar * grad_g_Bar * grad_psi_i_e * psi_j_phi
+ 0.5 * De_Bar * grad_h_Bar * g_Bar * grad_psi_i_e * psi_j_phi
+ 0.5 * De_Bar * h_Bar * h_Bar * grad_psi_i_e * grad_psi_j_e
+ 0.5 * De_Bar * 2 * h_Bar * grad_h_Bar * grad_psi_i_e * psi_j_e
)
* fe_values.JxW(q);
}
local_rhs(i) += ((0.5 * g_Bar * g_Bar - 0.5 * C_0) * psi_i_phi * old_phi
+ (1 / time_step) * 3 * psi_i_phi * old_phi
- 0.5 * eps * eps * k_0 * grad_psi_i_phi * grad_old_phi
+ 0.5 * (H_Bar * grad_old_phi) * (H_Bar * grad_psi_i_phi)
+ 0.5 * g_Bar * h_Bar * psi_i_phi * old_e
- g_Bar * psi_i_phi * old_q
- old_U * (H_Bar * grad_psi_i_phi)
- 0.5 * C_1 * De_Bar * grad_psi_i_e * grad_old_e
+ (1./ time_step) * psi_i_e * old_e
- De_Bar * h_Bar * grad_psi_i_e * grad_old_q
- De_Bar * grad_h_Bar * grad_psi_i_e * old_q
+ 0.5 * De_Bar * g_Bar * h_Bar * grad_psi_i_e * grad_old_phi
+ 0.5 * De_Bar * grad_g_Bar * h_Bar * grad_psi_i_e * old_phi
+ 0.5 * De_Bar * g_Bar * grad_h_Bar * grad_psi_i_e * old_phi
+ 0.5 * De_Bar * h_Bar * h_Bar * grad_psi_i_e * grad_old_e
+ 0.5 * De_Bar * 2 * h_Bar * grad_h_Bar * grad_psi_i_e * old_e
) * fe_values.JxW(q);
//std::cout << "dof 12-15-1 "<< "|" << i << "|" << local_rhs(i) << std::endl;
}
//std::cout << "dof 12-15-2 "<< "|" << q << "|" << std::endl;
}
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>::assemble_system_H0 (const Vector<double> phi_0)
{
system_matrix_H=0;
system_rhs_H = 0;
QGauss<dim> quadrature_formula(degree+2);
FEValues<dim> fe_phi_values (fe_phi, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
FEValues<dim> fe_H_values (fe_H, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_hessians | update_JxW_values);
std::cout << "dof 12-0 " << std::endl;
const unsigned int dofs_per_cell_H = fe_H.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
Vector<double> local_rhs (dofs_per_cell_H);
FullMatrix<double> local_matrix (dofs_per_cell_H, dofs_per_cell_H);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell_H);
//std::vector<Vector<double> > solution_values(n_q_points, Vector<double>(2));
std::vector<Tensor<1, dim> > phi_0_gradients(n_q_points);
const FEValuesExtractors::Scalar phase (0);//also use as H_x
const FEValuesExtractors::Scalar energy (1);//also use as H_y
//std::cout << "dof 12-0 " << std::endl;
typename DoFHandler<dim>::active_cell_iterator
cell_H = dof_handler_H.begin_active(),
endc_H = dof_handler_H.end();
typename DoFHandler<dim>::active_cell_iterator
cell_phi = dof_handler_phi.begin_active();
for (; cell_H!=endc_H; ++cell_H, ++cell_phi)
{
local_rhs = 0;
local_matrix=0;
fe_phi_values.reinit (cell_phi);
fe_H_values.reinit (cell_H);
//fe_values.get_function_values(solution, solution_values);
fe_phi_values.get_function_gradients(phi_0, phi_0_gradients);
for (unsigned int q=0; q<n_q_points; ++q)
{
Tensor<1,dim> grad_new_phi;
for (unsigned int d=0; d<dim; ++d)
{
grad_new_phi[d] = phi_0_gradients[q][d];
}
for (unsigned int i=0; i<dofs_per_cell_H; ++i)
{
const double psi_i_H_x = fe_H_values[phase].value (i, q); //In fact intervarU is H_x
const double psi_i_H_y = fe_H_values[energy].value (i, q); //In fact phase is H_y
for (unsigned int j=0; j<dofs_per_cell_H; ++j)
{
const double psi_j_H_x = fe_H_values[phase].value (j, q); //In fact intervarU is H_x
const double psi_j_H_y = fe_H_values[energy].value (j, q); //In fact phase is H_y
local_matrix(i,j) += (psi_j_H_x * psi_i_H_x
+psi_j_H_y * psi_i_H_y )
* fe_H_values.JxW(q);
}
local_rhs(i) += (psi_i_H_x * HFunction(grad_new_phi)[0]
+psi_i_H_y * HFunction(grad_new_phi)[1]) * fe_H_values.JxW(q);
}
}
cell_H->get_dof_indices (local_dof_indices);
constraints_H.distribute_local_to_global(local_matrix, local_rhs, local_dof_indices, system_matrix_H, system_rhs_H);
}
}
template <int dim>
void
StokesProblem<dim>::refine_mesh (const unsigned int max_grid_level,
const unsigned int min_grid_level)
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
FEValuesExtractors::Scalar phase(0);
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1>(degree+1),
typename FunctionMap<dim>::type(),
old_solution,
estimated_error_per_cell,
fe.component_mask(phase));
GridRefinement::refine_and_coarsen_fixed_number (triangulation,
estimated_error_per_cell,
0.2,0.1);
if (triangulation.n_levels() > max_grid_level)
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(max_grid_level);
cell != triangulation.end(); ++cell)
cell->clear_refine_flag ();
if (triangulation.n_levels() < min_grid_level)
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(min_grid_level);
cell != triangulation.end(); ++cell)
cell->clear_coarsen_flag ();
triangulation.execute_coarsening_and_refinement ();
}
template <int dim>
void StokesProblem<dim>::solve_1 ()
{
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix);
A_direct.vmult (solution_1, system_rhs);
}
template <int dim>
void StokesProblem<dim>::solve_H ()
{
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix_H);
A_direct.vmult (new_H, system_rhs_H);
}
template <int dim>
void StokesProblem<dim>::run ()
{
const unsigned int n_cycles = 1;
double entropy = 0, energy_volume = 0, energy_rate = 0, entropy_rate = 0, old_entropy = 0;
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation_active, 0, 4.5);
triangulation_active.refine_global (8);
GridGenerator::flatten_triangulation (triangulation_active, triangulation);
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
if (cell->face(f)->at_boundary())
{
if (std::fabs(cell->face(f)->center()(0) - (4.5)) < 1e-12)
cell->face(f)->set_boundary_id (1);
else if(std::fabs(cell->face(f)->center()(1) - (0)) < 1e-12)
cell->face(f)->set_boundary_id (2);
else if(std::fabs(cell->face(f)->center()(1) - (4.5)) < 1e-12)
cell->face(f)->set_boundary_id (3);
else if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
{
cell->face(f)->set_boundary_id (0);
}
}
}
}
}
else
{
//triangulation.refine_global (1);
}
setup_dofs ();
std::cout << "start" << cycle << std::endl;
VectorTools::project (dof_handler_phi, constraints, QGauss<dim>(degree+2),
InitialValues_phi<dim>(),
phi_0);
VectorTools::project (dof_handler_e, constraints, QGauss<dim>(degree+2),
InitialValues_e<dim>(),
e_0);
VectorTools::project (dof_handler_T, constraints, QGauss<dim>(degree+2),
InitialValues_T<dim>(),
T_0);
VectorTools::project (dof_handler_g, constraints, QGauss<dim>(degree+2),
InitialValues_g<dim>(),
g_0);
VectorTools::project (dof_handler_h, constraints, QGauss<dim>(degree+2),
InitialValues_h<dim>(),
h_0);
VectorTools::project (dof_handler_De, constraints, QGauss<dim>(degree+2),
InitialValues_De<dim>(),
De_0);
VectorTools::project (dof_handler_q, constraints_q, QGauss<dim>(degree+2),
InitialValues_q<dim>(),
q_0);
VectorTools::project (dof_handler_U, constraints_U, QGauss<dim>(degree+2),
InitialValues_U<dim>(),
old_U);
std::cout << "assemble_H_0" << cycle << std::endl;
assemble_system_H0 (phi_0);
std::cout << "solve_H_0" << cycle << std::endl;
solve_H ();
constraints_H.distribute (new_H);
older_H = new_H;
old_H = new_H;
older_g = g_0;
old_g = g_0;
older_h = h_0;
old_h = h_0;
older_De = De_0;
old_De = De_0;
// older_tau = tau_0;
// old_tau = tau_0;
older_T = T_0;
old_T = T_0;
std::cout << "start1" << cycle << std::endl;
assemble_system_1 (time= 1*time_step, phi_0, e_0, old_U, q_0, De_0, De_0, g_0, g_0, h_0, h_0, old_H, older_H);
std::cout << "assemble_1" << cycle << std::endl;
solve_1 ();
std::cout << "solve_1" << cycle << std::endl;
constraints.distribute (solution_1);
std::cout << "constraints_1" << cycle << std::endl;
old_solution = solution_1;
timestep_number=2;
std::cout << "time_to_spot" << timestep_number << std::endl;
refine_mesh (9, 6);
std::cout << "refine" << timestep_number << std::endl;
setup_dofs ();
std::cout << "setup_dofs" << timestep_number << std::endl;
}
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step22;
StokesProblem<2> flow_problem(1,1,1,1);
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;
}
