>
> If the reference solution is 0, you can't see if you are missing a
> constant factor in your equation even when the rate of convergence looks
> correct.
>
I understand thank you !
> Yes, but then you also have to be careful with the boundary conditions. If
> you have an analytical solution with homogeneous Dirichlet boundary
> conditions for \phi,
> I would definitely try that first.
>
It already has homogeneous Dirichlet boundary conditions. \phi is set to 0
at the boundaries as well as all the other variables (u, v, P).
> The ansatz space is the space spanned by the basis functions you are using
> for your solution variables. In case the mapping from the reference cell to
> the
> real cell is affine and you use polynomial base functions (on the
> reference element), the ansatz space is polynomial as well. Then, it is
> easier to manufacture
> a solution that is contained in the ansatz space.
>
I'm not sure to understand how to manufacture a solution for my problem. I
have read step 7 a few times and I understand but for a coupled system like
mine, I don't see how to do it.
Should I choose a solution with u and v that already respect the
icompressibility equation ?
Should I keep the same phi or change it for another one ?
Should I keep the same boundary conditions (everything = 0) on the sides ?
>From the look of it, the problem comes from the stokes equation. In my
simple case described above, the problem in phi is just a poisson equation
and I checked that the solution is the correct one (with a good convergence
rate).
I also completely erased my assembly of the stokes problem and re-wrote it.
It did not change a thing.
By the way, i'm now confident that the problem comes from my dealii code
and not from my mathematica code. I have implemented the same code using
Fenics, the python library, and it gives the exact same results as the
mathematica code. I used the same weak formulation in both the dealii and
fenics codes.
Fenics and Mathematica use triangles for a mesh and Dealii use squares.
Could this cause a difference in the solution, especially knowing that i'm
working on a sphere ? I'm currently comparing the same equations on a
square mesh. We will see how it goes.
I'm joining a cpp file with the parts of the code that only correspond to
the stokes problem. I pretty much used the same format as in the tutorial
so it will look familiar to anyone that read or wrote them.
Will send another e-mail in the afternoon after I studied the problem on a
square mesh.
Thanks again to you Daniel for taking the time to help me, and for having
developed dealii which is by far the best FE solver for complex problems I
have used so far.
--
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/7b2ed826-bf59-42e9-97bd-d523fd820c26%40googlegroups.com.
//Définition de la classe
Elfe::Elfe (
double n,
double s,
double A,
double a,
double X,
unsigned int cycle,
std::string basename):
n(n),
s(s),
A(A),
a(a),
X(X),
basename(basename),
cycle(cycle),
phi_fe (2),
phi_dof_handler (triangulation),
stokes_fe(FE_Q<2>(2), 2,
FE_Q<2>(1), 1),
stokes_dof_handler(triangulation)
{
}
void Elfe::stokes_setup_dofs ()
{ //Fonction qui setup les degree of freedom pour stokes
//On génère le dofhandler
stokes_dof_handler.distribute_dofs (stokes_fe);
//On réorganise les dofs pour mettre u, v, p
std::vector<unsigned int> block_component(3, 0);
block_component[2] = 1;
DoFRenumbering::component_wise (stokes_dof_handler, block_component);
std::vector<types::global_dof_index> dofs_per_block;
dofs_per_block.resize (2);
DoFTools::count_dofs_per_block (stokes_dof_handler, dofs_per_block, block_component);
unsigned int dof_u = dofs_per_block[0];
unsigned int dof_p = dofs_per_block[1];
//On génère les conditions aux bords
stokes_constraints.clear();
DoFTools::make_hanging_node_constraints(stokes_dof_handler, stokes_constraints);
VectorTools::interpolate_boundary_values(stokes_dof_handler,
0,
ZeroFunction<2>(3),
stokes_constraints
);
stokes_constraints.close();
//On crée le sparsity pattern
BlockDynamicSparsityPattern dsp (dofs_per_block, dofs_per_block);
DoFTools::make_sparsity_pattern (stokes_dof_handler, dsp, stokes_constraints, false);
stokes_sparsity_pattern.copy_from (dsp);
//On resize les objets à la bonne taille
stokes_system_matrix.reinit (stokes_sparsity_pattern);
stokes_solution.reinit (dofs_per_block);
stokes_system_rhs.reinit (dofs_per_block);
}
void Elfe::solve_stokes()
{ //Fonction qui résout le problème de Stokes
std::cout << BOLDBLUE << "\n Résolution Stokes " << RESET
<< std::endl;
//On assemble le système
stokes_assemble_system();
std::cout << BLUE << "Inversion" << RESET
<< std::endl;
// On l'inverse
SparseDirectUMFPACK A_direct;
A_direct.initialize(stokes_system_matrix);
A_direct.vmult(stokes_solution, stokes_system_rhs);
stokes_constraints.distribute(stokes_solution);
std::cout << "\x1b[A\x1b[A\x1b[A";
}
void Elfe::stokes_assemble_system ()
{ //Fonction qui assemble le problème de Stokes
std::cout << BLUE << "Assemblage Stokes " << RESET
<< std::endl;
//On réinitialise la matrice et le rhs
stokes_system_matrix = 0;
stokes_system_rhs = 0;
//On charge de quoi extraire les valeurs
QGauss<2> quadrature_formula(4);
const unsigned int dofs_per_cell = stokes_fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
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);
//Stokes values
FEValues<2> stokes_fe_values(stokes_fe, quadrature_formula,
update_values | update_gradients | update_JxW_values);
double u_i, v_i;
double u_i_x, u_i_y, v_i_x, v_i_y;
double u_j, v_j;
double u_j_x, u_j_y, v_j_x, v_j_y;
double P_i, P_j;
const FEValuesExtractors::Vector Velocities(0);
const FEValuesExtractors::Scalar Pressure(2);
//Phi values
FEValues<2> phi_fe_values(phi_fe, quadrature_formula,
update_values | update_gradients | update_JxW_values);
double phi, phi_x, phi_y;
double DX;
double philin = 2*pi*n;
double phiscale = 2*pi*s;
std::vector<Tensor<1, 2> > phi_solution_gradients(n_q_points);
std::vector<double> phi_solution_values(n_q_points);
//Il nous faut 2 cell iterator, un pour phi, un pour stokes
typename DoFHandler<2>::active_cell_iterator
stokes_cell = stokes_dof_handler.begin_active(),
stokes_endc = stokes_dof_handler.end();
typename DoFHandler<2>::active_cell_iterator
phi_cell = phi_dof_handler.begin_active();
for (; stokes_cell!=stokes_endc; ++stokes_cell, ++phi_cell)
{ //Pour chaque cell
stokes_fe_values.reinit(stokes_cell);
phi_fe_values.reinit (phi_cell);
local_matrix = 0;
local_rhs = 0;
//On récupère les valeurs pour les angles phis
phi_fe_values.get_function_values(phi_solution,
phi_solution_values);
phi_fe_values.get_function_gradients(phi_solution,
phi_solution_gradients);
for (unsigned int q=0; q<n_q_points; ++q)
{ //Pour chaque point de quadrature
DX = stokes_fe_values.JxW(q);
phi = phi_solution_values[q];
phi_x = phi_solution_gradients[q][0];
phi_y = phi_solution_gradients[q][1];
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
u_i = stokes_fe_values[Velocities].value (i, q)[0];
v_i = stokes_fe_values[Velocities].value (i, q)[1];
u_i_x = stokes_fe_values[Velocities].gradient(i, q)[0][0];
u_i_y = stokes_fe_values[Velocities].gradient(i, q)[0][1];
v_i_x = stokes_fe_values[Velocities].gradient(i, q)[1][0];
v_i_y = stokes_fe_values[Velocities].gradient(i, q)[1][1];
P_i = stokes_fe_values[Pressure] .value (i, q);
local_rhs(i) += (
u_i*(
8*philin*phiscale*X
*(cos(2*phiscale*phi)*phi_x + sin(2*phiscale*phi)*phi_y)
+ 2*phiscale*phiscale*phi_x*(-4*philin/phiscale)
)
+ v_i*(
8*philin*phiscale*X
*(sin(2*phiscale*phi)*phi_x - cos(2*phiscale*phi)*phi_y)
+ 2*phiscale*phiscale*phi_y*(-4*philin/phiscale)
)
)*DX;
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
u_j = stokes_fe_values[Velocities].value (j, q)[0];
v_j = stokes_fe_values[Velocities].value (j, q)[1];
u_j_x = stokes_fe_values[Velocities].gradient(j, q)[0][0];
u_j_y = stokes_fe_values[Velocities].gradient(j, q)[0][1];
v_j_x = stokes_fe_values[Velocities].gradient(j, q)[1][0];
v_j_y = stokes_fe_values[Velocities].gradient(j, q)[1][1];
P_j = stokes_fe_values[Pressure] .value (j, q);
local_matrix(i, j) += (
//Classical Stokes terms
- u_i_x*u_j_x - u_i_y*u_j_y - v_i_x*v_j_x - v_i_y*v_j_y
+ (u_i_x+v_i_y)*P_j
+ P_i*(u_j_x+v_j_y)
//Leslie Asymetrical terms
+ phiscale/a
*(v_i_x-u_i_y)
*(u_j*phi_x + v_j*phi_y + 1/2/phiscale*(u_j_y-v_j_x))
//Terme provenant du laplacien phi
-2*phiscale*phiscale/a
*(u_i*phi_x + v_i*phi_y)
*(u_j*phi_x + v_j*phi_y + 1/2/phiscale*(u_j_y-v_j_x))
)*DX;
}
}
}
stokes_cell->get_dof_indices (local_dof_indices);
//On applique les boundary conditions
stokes_constraints.distribute_local_to_global(local_matrix,
local_rhs,
local_dof_indices,
stokes_system_matrix,
stokes_system_rhs);
}
}
void Elfe::output_stokes () const
{ //Fonction qui sort u, v, P au format GPL
std::cout << BLUE << " Enregistrement de Stokes\n" << RESET;
std::stringstream stream;
stream << basename;
std::string fileName = stream.str();
DataOut<2> data_out;
data_out.attach_dof_handler (stokes_dof_handler);
data_out.add_data_vector (stokes_solution, "solution");
data_out.build_patches ();
std::ofstream output (fileName);
data_out.write_gnuplot (output);
}
