Dear Wolfgang,
Sorry for an additional question.
I have tried to implement by myself something similar to what is done in
ASPECT (I guess?, it is a bit harder for me to grasp ASPECT concepts
because of the introspection and etc.).
However, I seem to have failed miserably at one point.
My understanding is that I can get
A x = f- Pf
By:
- First I calculate my projection operator. It is the integral of the
pressure shape function over ALL elements associated with a DOF. I
calculate it by looping over the elements and integrating with a quadrature
(see code below)
- Then I assemble my system and my right-hand side as usual
- Before I solve my linear system, I modify my Right-hand-side so that it
does not include the constant mode. I litteraly calculate f= f- Pf
- I solve my linear system for this new RHS.
However, my GMRES stops very quickly after a certain number of newton
iteration with the following : AztecOO::Iterate error code -4: GMRES
Hessenberg ill-conditioned
Clearly I am doing something wrong / I understood something wrong regarding
this, because I know the rest of my code works.
In my system there is no block matrices and everything is lumped in the
same matrix. I havent re-ordered by components either.
Calculation of the projection:
template <int dim>
void
GLSNavierStokesSolver<dim>::initializePressureRHSCorrection()
{
TimerOutput::Scope t(computing_timer, "pressure_RHS_projector");
pressure_shape_function_integrals=0;
QGauss<dim> quadrature_formula(degreeQuadrature_);
FEValues<dim> fe_values (fe,
quadrature_formula,
update_values |
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 FEValuesExtractors::Scalar pressure (dim);
Vector<double>
local_pressure_shape_function_integrals(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
if (cell->is_locally_owned())
{
fe_values.reinit(cell);
local_pressure_shape_function_integrals=0;
for (unsigned int i=0; i<dofs_per_cell;++i)
{
for (unsigned int q=0; q<n_q_points; ++q)
{
const double phi_p = fe_values[pressure].value(i, q);
local_pressure_shape_function_integrals(i) += phi_p *
fe_values.JxW(q);
}
}
cell->get_dof_indices (local_dof_indices);
zero_constraints.distribute_local_to_global(
local_pressure_shape_function_integrals,
local_dof_indices,
pressure_shape_function_integrals);
}
}
pressure_shape_function_integrals.compress (VectorOperation::add);
}
Projection of the System RHS onto it's space without the constant:
template <int dim>
void GLSNavierStokesSolver<dim>::correctRHSClosedSystem()
{
// calculate projection of system_rhs on pressure_shape_function_integrals
for (unsigned int dof_i=0 ; dof_i < pressure_shape_function_integrals.size
() ; ++dof_i)
pressure_projection[dof_i] = pressure_shape_function_integrals[dof_i]*
system_rhs[dof_i];
// Add it to RHS
system_rhs.add(-1.,pressure_projection);
}
Thank you again for everything,
This is greatly appreciated.
Bruno
On Tuesday, 5 March 2019 16:10:09 UTC-5, Wolfgang Bangerth wrote:
>
>
> > A quick question. I think I understand what is done in
> >
> https://github.com/geodynamics/aspect/blob/master/source/simulator/helper_functions.cc#L1041
>
> > Weirdfully, I found the "pickaxe" version easier to understand (Love the
> > comments by the way, this is awesome)
>
> :-) I suspect that Timo gets credit for that one!
>
>
> > A question I have is related to the _pressure_shape_function_integrals_
> member.
> > My understanding is that this is the integration over the element cellI
> of the
> > shape function associated with pressure.
>
> I believe it's actually the integral over all cells the shape function is
> associated with.
>
>
> > However, I have looked over the ASPECT documentation and I have not been
> able
> > to find where this is implemented. Is it a vector that is manually
> filled or
> > is there a helper function that can be used to automatically generate
> it?
>
> Yes, that required a bit of searching. This is filled here:
>
> https://github.com/geodynamics/aspect/blob/master/source/simulator/assemblers/stokes.cc#L601
>
>
> Best
> W.
>
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> <javascript:>
> www: http://www.math.colostate.edu/~bangerth/
>
>
--
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.
