Dear Prof. Arndt,
I have used SolutionTransfer as you say, but if I set a phi_0 fixed, then
project it into finite element space and get a vector phi_0_h, then I get
phi_0_h/2, phi_0_h/4, phi_0_h/8 by using SolutionTransfer, but the norm
of (phi_0_h- phi_0_h/2), (phi_0_h/2- phi_0_h/4), (phi_0_h/4- phi_0_h/8) are
so big, and I cannot get exact value when I compute the convergence rate.
So I wander how to deal with it?
And the attached is the test code, and the results is:
phi_0_norm_58: 0.335377 phi_0_norm_68: 0.317158 phi_0_norm_78: 0.277344
Thank you very much!
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 dealii+unsubscr...@googlegroups.com.
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 <string>
#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);
void run ();
private:
void setup_dofs_pbc_rf (const int refine_number);
void phi_0_trans ();
double norm_compute_phi_0 (const Vector<double> solu_1,
const Vector<double> solu_2);
void create_mesh();
void refine_mesh ();
const unsigned int degree;
Triangulation<dim> triangulation_active;
Triangulation<dim> triangulation;
FE_Q<dim> fe_phi;
DoFHandler<dim> dof_handler_phi;
ConstraintMatrix constraints_phi;
Vector<double> phi_0, phi_0_5, phi_0_6, phi_0_7, phi_0_8, tmp_phi_0_5, tmp_phi_0_6, tmp_phi_0_7, tmp_phi_0_8;
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;
};
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int phi_degree)
:
degree (phi_degree),
fe_phi (degree),
dof_handler_phi (triangulation),
time_step (1e-4),
timestep_number (1),
A (100.0),
B (1.0),
sita_0 (1*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.1, 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] - 0.28125) * (p[0] - 0.28125) + (p[1] - 0.28125) * (p[1] - 0.28125)) ;
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 InitialValuesConv : public Function<dim>
{
public:
InitialValuesConv () : Function<dim>() {}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
};
template <int dim>
double
InitialValuesConv<dim>::value (const Point<dim> &p,
const unsigned int component) const
{
return ZeroFunction<dim>().value (p, component);
}
template <int dim>
void StokesProblem<dim>::setup_dofs_pbc_rf (const int refine_number)
{
dof_handler_phi.distribute_dofs (fe_phi);
DoFRenumbering::component_wise (dof_handler_phi);
const unsigned int n_phi = dof_handler_phi.n_dofs();
constraints_phi.close ();
switch(refine_number)
{
case 5:
phi_0_5.reinit (n_phi);
break;
case 6:
phi_0_6.reinit (n_phi);
break;
case 7:
phi_0_7.reinit (n_phi);
break;
case 8:
phi_0_8.reinit (n_phi);
break;
default:
std::cout << "-----------------------------------------------------------------------------wrong" << std::endl;
}
tmp_phi_0_5.reinit (n_phi);
tmp_phi_0_6.reinit (n_phi);
tmp_phi_0_7.reinit (n_phi);
tmp_phi_0_8.reinit (n_phi);
phi_0.reinit (n_phi);
}
template <int dim>
void
StokesProblem<dim>::phi_0_trans ()
{
SolutionTransfer<dim, Vector<double> > phi_0_5_trans(dof_handler_phi);
SolutionTransfer<dim, Vector<double> > phi_0_6_trans(dof_handler_phi);
SolutionTransfer<dim, Vector<double> > phi_0_7_trans(dof_handler_phi);
SolutionTransfer<dim, Vector<double> > phi_0_8_trans(dof_handler_phi);
triangulation.prepare_coarsening_and_refinement();
phi_0_5_trans.prepare_for_coarsening_and_refinement(phi_0_5);
phi_0_6_trans.prepare_for_coarsening_and_refinement(phi_0_6);
phi_0_7_trans.prepare_for_coarsening_and_refinement(phi_0_7);
phi_0_8_trans.prepare_for_coarsening_and_refinement(phi_0_8);
triangulation.execute_coarsening_and_refinement ();
phi_0_5_trans.interpolate(phi_0_5, tmp_phi_0_5);
phi_0_6_trans.interpolate(phi_0_6, tmp_phi_0_6);
phi_0_7_trans.interpolate(phi_0_7, tmp_phi_0_7);
phi_0_8_trans.interpolate(phi_0_8, tmp_phi_0_8);
phi_0_5 = tmp_phi_0_5;
phi_0_6 = tmp_phi_0_6;
phi_0_7 = tmp_phi_0_7;
phi_0_8 = tmp_phi_0_8;
constraints_phi.distribute(phi_0_5);
constraints_phi.distribute(phi_0_6);
constraints_phi.distribute(phi_0_7);
constraints_phi.distribute(phi_0_8);
std::cout << " phi_0_5: " << phi_0_5.size() << " phi_0_6: " << phi_0_6.size() << " phi_0_7: " << phi_0_7.size() << " phi_0_8: " << phi_0_8.size() << std::endl;
}
template <int dim>
void StokesProblem<dim>::create_mesh()
{
GridGenerator::hyper_cube (triangulation, 0, 0.5625);
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) - (0.5625)) < 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) - (0.5625)) < 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);
}
}
}
}
triangulation.refine_global (5);
std::ofstream out ("./grid_0.vtk");
GridOut grid_out;
grid_out.write_vtk (triangulation, out);
std::cout << "Grid has been written to grid.vtk" << std::endl;
}
template <int dim>
double StokesProblem<dim>::norm_compute_phi_0 (const Vector<double> solu_1,
const Vector<double> solu_2)
{
QGauss<dim> quadrature_formula(degree+4);
Vector<double> difference_phi;
difference_phi.reinit(solu_1.size());
for (unsigned int i=0; i<solu_1.size(); ++i)
{
difference_phi[i] =solu_1[i] - solu_2[i];
}
FEValues<dim> fe_phi_values (fe_phi, quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
Vector<float> difference_per_cell(triangulation.n_active_cells());
InitialValuesConv<dim> zero_for_conv;
VectorTools::integrate_difference(dof_handler_phi,
difference_phi,
zero_for_conv,
difference_per_cell,
QGauss<dim>(degree + 3),
VectorTools::L2_norm);
const double L2_error = VectorTools::compute_global_error(
triangulation, difference_per_cell, VectorTools::L2_norm);
return L2_error;
}
template <int dim>
void StokesProblem<dim>::run ()
{
const unsigned int n_cycles = 8, cycle=0;
double entropy = 0, energy_volume = 0, energy_rate = 0, entropy_rate = 0, old_entropy = 0;
create_mesh ();
//BlockVector<double> solu_5, solu_6, solu_7, solu_8;
for (unsigned int cycle=5; cycle<=n_cycles; ++cycle)
{
std::cout << "Cycle " << cycle << " time_step: " << time_step << std::endl;
int refine_number = 5;
if (cycle == 5)
{
std::cout << "Cycle " << cycle << " refine_number: " << refine_number << std::endl;
}
else
{
triangulation.refine_global (1);
refine_number= cycle;
std::cout << "Cycle: " << cycle << " refine_number: " << refine_number << std::endl;
}
// setup_dofs ();
setup_dofs_pbc_rf (refine_number);
std::cout << "start" << cycle << std::endl;
VectorTools::project (dof_handler_phi, constraints_phi, QGauss<dim>(degree+2),
InitialValues_phi<dim>(),
phi_0);
switch(cycle)
{
case 5:
phi_0_5 = phi_0;
std::cout << "---------------------------------------------------------------------------solu_5 "<< phi_0_5.size() << std::endl;
break;
case 6:
phi_0_6 = phi_0;
std::cout << "----------------------------------------------------------------------------solu_6 "<< phi_0_6.size() << std::endl;
break;
case 7:
phi_0_7 = phi_0;
std::cout << "----------------------------------------------------------------------------solu_7 "<< phi_0_7.size() << std::endl;
break;
case 8:
phi_0_8 = phi_0;
std::cout << "----------------------------------------------------------------------------solu_8 "<< phi_0_8.size()<< std::endl;
break;
default:
std::cout << "-----------------------------------------------------------------------------wrong" << std::endl;
}
std::cout << "phi" << cycle << std::endl;
}
phi_0_trans ();
double phi_0_norm_58 = 1, phi_0_norm_78 = 1, phi_0_norm_68 = 1;
phi_0_norm_58 = norm_compute_phi_0(phi_0_5, phi_0_8);
phi_0_norm_68 = norm_compute_phi_0(phi_0_6, phi_0_8);
phi_0_norm_78 = norm_compute_phi_0(phi_0_7, phi_0_8);
std::cout << " phi_0_norm_58: " << phi_0_norm_58 << " phi_0_norm_68: " << phi_0_norm_68 << " phi_0_norm_78: " << phi_0_norm_78 << std::endl;
}
}
int main ()
{
try
{
using namespace dealii;
using namespace Step22;
StokesProblem<2> flow_problem(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;
}
