[cig-commits] commit 2013 by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Tue Nov 19 12:43:40 PST 2013
Revision 2013
Only do the correction of the rhs of the Stokes equation's second component if the material is incompressible (as before) and we have no open boundaries or places where we prescribe a boundary velocity.
U trunk/aspect/include/aspect/simulator.h
U trunk/aspect/source/simulator/assembly.cc
U trunk/aspect/source/simulator/core.cc
U trunk/aspect/source/simulator/helper_functions.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2013&peg=2013
Diff:
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h 2013-11-18 18:22:35 UTC (rev 2012)
+++ trunk/aspect/include/aspect/simulator.h 2013-11-19 20:43:14 UTC (rev 2013)
@@ -247,7 +247,7 @@
*/
/**
- * @name Parameters that have to do with compositional field
+ * @name Parameters that have to do with compositional fields
* @{
*/
unsigned int n_compositional_fields;
@@ -1100,7 +1100,7 @@
double global_Omega_diameter;
double global_volume;
- MeshRefinement::Manager<dim> mesh_refinement_manager;
+ MeshRefinement::Manager<dim> mesh_refinement_manager;
const MappingQ<dim> mapping;
@@ -1111,10 +1111,28 @@
ConstraintMatrix constraints;
ConstraintMatrix current_constraints;
+ /**
+ * The latest correction computed by normalize_pressure(). We store this so
+ * we can undo the correction in denormalize_pressure().
+ */
double pressure_adjustment;
+
+ /**
+ * Scaling factor for the pressure as explained in the Kronbichler/Heister/Bangerth
+ * paper to ensure that the linear system that results from the Stokes equations
+ * is well conditioned.
+ */
double pressure_scaling;
/**
+ * A variable that determines whether we need to do the correction of
+ * the Stokes right hand side vector to ensure that the average divergence
+ * is zero. This is necessary for compressible models, but only if there
+ * are no in/outflow boundaries.
+ */
+ bool do_pressure_rhs_compatibility_modification;
+
+ /**
* @}
*/
Modified: trunk/aspect/source/simulator/assembly.cc
===================================================================
--- trunk/aspect/source/simulator/assembly.cc 2013-11-18 18:22:35 UTC (rev 2012)
+++ trunk/aspect/source/simulator/assembly.cc 2013-11-19 20:43:14 UTC (rev 2013)
@@ -375,7 +375,8 @@
template <int dim>
struct StokesSystem : public StokesPreconditioner<dim>
{
- StokesSystem (const FiniteElement<dim> &finite_element);
+ StokesSystem (const FiniteElement<dim> &finite_element,
+ const bool do_pressure_rhs_compatibility_modification);
StokesSystem (const StokesSystem<dim> &data);
Vector<double> local_rhs;
@@ -386,11 +387,15 @@
template <int dim>
StokesSystem<dim>::
- StokesSystem (const FiniteElement<dim> &finite_element)
+ StokesSystem (const FiniteElement<dim> &finite_element,
+ const bool do_pressure_rhs_compatibility_modification)
:
StokesPreconditioner<dim> (finite_element),
local_rhs (finite_element.dofs_per_cell),
- local_pressure_shape_function_integrals (finite_element.dofs_per_cell)
+ local_pressure_shape_function_integrals (do_pressure_rhs_compatibility_modification ?
+ finite_element.dofs_per_cell
+ :
+ 0)
{}
@@ -401,7 +406,7 @@
:
StokesPreconditioner<dim> (data),
local_rhs (data.local_rhs),
- local_pressure_shape_function_integrals (data.local_pressure_shape_function_integrals)
+ local_pressure_shape_function_integrals (data.local_pressure_shape_function_integrals.size())
{}
@@ -1056,7 +1061,7 @@
if (rebuild_stokes_matrix)
data.local_matrix = 0;
data.local_rhs = 0;
- if (is_compressible)
+ if (do_pressure_rhs_compatibility_modification)
data.local_pressure_shape_function_integrals = 0;
// we only need the strain rates for the viscosity,
@@ -1129,7 +1134,7 @@
0)
)
* scratch.finite_element_values.JxW(q);
- if (is_compressible)
+ if (do_pressure_rhs_compatibility_modification)
for (unsigned int i=0; i<dofs_per_cell; ++i)
data.local_pressure_shape_function_integrals(i) += scratch.phi_p[i] * scratch.finite_element_values.JxW(q);
}
@@ -1155,7 +1160,7 @@
data.local_dof_indices,
system_rhs);
- if (material_model->is_compressible())
+ if (do_pressure_rhs_compatibility_modification)
current_constraints.distribute_local_to_global (data.local_pressure_shape_function_integrals,
data.local_dof_indices,
pressure_shape_function_integrals);
@@ -1172,7 +1177,7 @@
system_matrix = 0;
system_rhs = 0;
- if (material_model->is_compressible())
+ if (do_pressure_rhs_compatibility_modification)
pressure_shape_function_integrals = 0;
const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
@@ -1208,12 +1213,13 @@
UpdateFlags(0))),
parameters.n_compositional_fields),
internal::Assembly::CopyData::
- StokesSystem<dim> (finite_element));
+ StokesSystem<dim> (finite_element,
+ do_pressure_rhs_compatibility_modification));
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
- if (material_model->is_compressible())
+ if (do_pressure_rhs_compatibility_modification)
pressure_shape_function_integrals.compress(VectorOperation::add);
Modified: trunk/aspect/source/simulator/core.cc
===================================================================
--- trunk/aspect/source/simulator/core.cc 2013-11-18 18:22:35 UTC (rev 2012)
+++ trunk/aspect/source/simulator/core.cc 2013-11-19 20:43:14 UTC (rev 2013)
@@ -138,6 +138,7 @@
rebuild_stokes_preconditioner (true)
{
computing_timer.enter_section("Initialization");
+
// first do some error checking for the parameters we got
{
// make sure velocity boundary indicators don't appear in multiple lists
@@ -247,8 +248,35 @@
velocity_boundary_conditions[p->first].reset (bv);
}
+ // determine how to treat the pressure. we have to scale it for the solver
+ // to make velocities and pressures of roughly the same (numerical) size,
+ // and we may have to fix up the right hand side vector before solving for
+ // compressible models if there are no in-/outflow boundaries
pressure_scaling = material_model->reference_viscosity() / geometry_model->length_scale();
+ std::set<types::boundary_id> open_velocity_boundary_indicators
+ = geometry_model->get_used_boundary_indicators();
+ for (std::map<types::boundary_id,std::pair<std::string,std::string> >::const_iterator
+ p = parameters.prescribed_velocity_boundary_indicators.begin();
+ p != parameters.prescribed_velocity_boundary_indicators.end();
+ ++p)
+ open_velocity_boundary_indicators.erase (p->first);
+ for (std::set<types::boundary_id>::const_iterator
+ p = parameters.zero_velocity_boundary_indicators.begin();
+ p != parameters.zero_velocity_boundary_indicators.end();
+ ++p)
+ open_velocity_boundary_indicators.erase (*p);
+ for (std::set<types::boundary_id>::const_iterator
+ p = parameters.tangential_velocity_boundary_indicators.begin();
+ p != parameters.tangential_velocity_boundary_indicators.end();
+ ++p)
+ open_velocity_boundary_indicators.erase (*p);
+ do_pressure_rhs_compatibility_modification = (material_model->is_compressible()
+ &&
+ (parameters.prescribed_velocity_boundary_indicators.size() == 0)
+ &&
+ (open_velocity_boundary_indicators.size() == 0));
+
// make sure that we don't have to fill every column of the statistics
// object in each time step.
statistics.set_auto_fill_mode(true);
@@ -746,7 +774,7 @@
old_old_solution.reinit(introspection.index_sets.system_relevant_partitioning, mpi_communicator);
current_linearization_point.reinit (introspection.index_sets.system_relevant_partitioning, mpi_communicator);
- if (material_model->is_compressible())
+ if (do_pressure_rhs_compatibility_modification)
pressure_shape_function_integrals.reinit (introspection.index_sets.system_partitioning, mpi_communicator);
rebuild_stokes_matrix = true;
Modified: trunk/aspect/source/simulator/helper_functions.cc
===================================================================
--- trunk/aspect/source/simulator/helper_functions.cc 2013-11-18 18:22:35 UTC (rev 2012)
+++ trunk/aspect/source/simulator/helper_functions.cc 2013-11-19 20:43:14 UTC (rev 2013)
@@ -664,10 +664,13 @@
if (parameters.use_locally_conservative_discretization)
AssertThrow(false, ExcNotImplemented());
- const double mean = vector.block(1).mean_value();
- const double correction = -mean*vector.block(1).size()/global_volume;
+ if (do_pressure_rhs_compatibility_modification)
+ {
+ const double mean = vector.block(1).mean_value();
+ const double correction = -mean*vector.block(1).size()/global_volume;
- vector.block(1).add(correction, pressure_shape_function_integrals.block(1));
+ vector.block(1).add(correction, pressure_shape_function_integrals.block(1));
+ }
}
More information about the CIG-COMMITS
mailing list