> email to dea...@googlegroups.com
> <mailto:dea...@googlegroups.com>.
I have tried a simple testcase for my problem.I don't know if it's suitable as a testcase but if I use X=0 in my equations above, I know the solution.The solutions is u = v = 0, \phi = (1-x^2 -y^2) and the pressure is then P = 8 \phi_lin^2 (1-x²-y²)It fits with my boundary conditions which are u = v = p = \phi = 0 at the boundaries.Here are the results for dealii with a = 2 and phi_lin = 4*pi.Left is the L2 norm of u and v. As you can see it's going to 0 when i jump to a higher refinement cycle.
I'm a bit at a lose. Why does it work for X = 0 but not when I use a different X ?
I have tried a simple testcase for my problem.I don't know if it's suitable as a testcase but if I use X=0 in my equations above, I know the solution.The solutions is u = v = 0, \phi = (1-x^2 -y^2) and the pressure is then P = 8 \phi_lin^2 (1-x²-y²)It fits with my boundary conditions which are u = v = p = \phi = 0 at the boundaries.Here are the results for dealii with a = 2 and phi_lin = 4*pi.Left is the L2 norm of u and v. As you can see it's going to 0 when i jump to a higher refinement cycle.Is the rate of convergence as expected?
I'm a bit at a lose. Why does it work for X = 0 but not when I use a different X ?Are you sure that your deal.II solution is wrong in fact?
If u = v = 0 in your test case, there might still be a problem with these equations.
I would also try a case where you manufacture a solution and implement a right-hand side term such that this is the solution to your equations.
Since the ansatz space on non-affine mapped geometries is not polynomial anymore, I would also try a cartesian mesh for the case that the solution is containedin your ansatz space.
Best,Daniel
Is the rate of convergence as expected?I believe so. I did a linear regression of the convergence graphs in logscale and got a value of :-7.6 for the speed : I believe it corresponds to the -8 from second order finete element in 2D-3.2 for the pressure : I beliebe it corresponds to the -4 from the first order finite element in 2D
If u = v = 0 in your test case, there might still be a problem with these equations.I'm not sur to understand what you are saying here. Can you elaborate why you think there is a problem in my equations ?
Since the ansatz space on non-affine mapped geometries is not polynomial anymore, I would also try a cartesian mesh for the case that the solution is containedin your ansatz space.I don't really understand what you mean by ansatz space although I have looked it up on google. Do you suggest that I try the X=0 solution on a square mesh ?
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.
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.
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 thereal 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 manufacturea solution that is contained in the ansatz space.
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);
}
}
do not compute to what you think they do. The compiler will read this from left to right, and the first operation it sees is an integer divided by an integer. Integer division in C++ will return the value of the quotient rounded down (towards zero, if positive). There’s a more detailed explanation here:1/2/phiscale*(u_j_y-v_j_x)
On 25 Jul 2019, at 16:10, Félix Bunel <bunel...@gmail.com> wrote:
One final remark for today, if I remove the two asymmetrical terms in my equation so it becomes :
<Auto Generated Inline Image 1.png>
--
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+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a984ed95-fa81-4758-bc46-1357ec5ce2b5%40googlegroups.com.
<Auto Generated Inline Image 1.png>
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 ?
To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.
Okay. When I read you e-mail I sincerely hoped it wasn't the problem. But it was. Now the results match.
I feel so so stupid. Such a rookie mistake. I already made this mistake 3 years ago and forgot about it.
Thank you very much for taking the time to look through this part of the code. I own you big time.Thank you also Daniel for helping me. I will still try a manufactured solution to learn about this process.
Thanks also to everyone that contributed to dealii.
Its the perfect library for people that are not as stupid as me.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/f887ca46-262f-4805-8b94-ef68341cb455%40googlegroups.com.