Additional Stokes RHS #1547
Additional Stokes RHS #1547
Conversation
|
Can you update the test? |
include/aspect/assembly.h
Outdated
| @@ -711,6 +711,9 @@ namespace aspect | |||
| public SimulatorAccess<dim> | |||
| { | |||
| public: | |||
|
|
|||
| virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const; | |||
There was a problem hiding this comment.
what's this function gonna do?
| * dim components for velocity, 1 component for pressure, 1 component | ||
| * for compaction pressure (if melt transport is enabled). | ||
| */ | ||
| std::vector<std::vector<double> > rhs; |
There was a problem hiding this comment.
You call this class StokesRHS, so I suspect that the force terms should be Tensor<1,dim> objects, assuming you're talking about the force balance equation. If you want additional terms for the mass balance equation, maybe have another std::vector<double>.
If you really wanted to be generic and add terms for other equations as well, I would suggest to not document the order here, but use the same ordering as used by the introspection object (to which both the material models and assemblers have access).
| create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const | ||
| { | ||
| if (this->get_parameters().enable_additional_stokes_rhs | ||
| && !outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >()) |
There was a problem hiding this comment.
I'll keep fighting against using operator ! to compare pointers for non-Null-ness till I die! :-)
Please write the second condition as
&& (outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >() != NULL)
or similar.
| { | ||
| public: | ||
| AdditionalMaterialOutputsStokesRHS(const unsigned int n_points, | ||
| const unsigned int n_stokes_components) |
There was a problem hiding this comment.
how could n_stokes_components ever be different than dim+1?
There was a problem hiding this comment.
Normally dim+1 but dim+2 if you have melt transport enabled. That is also the reason I did not want to use
Tensor<1,dim> u;
double p;
double p_c;
or something of that matter.
There was a problem hiding this comment.
But then it's not just Stokes but everything.
| @@ -113,6 +130,9 @@ namespace aspect | |||
| const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size(); | |||
| const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points; | |||
|
|
|||
| MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> | |||
| *force = scratch.material_model_outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >(); | |||
source/simulator/melt.cc
Outdated
| @@ -88,6 +89,17 @@ namespace aspect | |||
| MeltEquations<dim>::create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &outputs) const | |||
| { | |||
| MeltHandler<dim>::create_material_model_outputs(outputs); | |||
|
|
|||
| if (this->get_parameters().enable_additional_stokes_rhs | |||
| && !outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >()) | |||
source/simulator/melt.cc
Outdated
| && !outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >()) | ||
| { | ||
| const unsigned int n_points = outputs.viscosities.size(); | ||
| const unsigned int n_stokes_comp = dim + 2; |
source/simulator/melt.cc
Outdated
| @@ -215,6 +227,8 @@ namespace aspect | |||
| const unsigned int p_c_component_index = introspection.variable("compaction pressure").first_component_index; | |||
|
|
|||
| MaterialModel::MeltOutputs<dim> *melt_outputs = scratch.material_model_outputs.template get_additional_output<MaterialModel::MeltOutputs<dim> >(); | |||
| MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> | |||
| *force = scratch.material_model_outputs.template get_additional_output<MaterialModel::AdditionalMaterialOutputsStokesRHS<dim> >(); | |||
source/simulator/melt.cc
Outdated
| : std::numeric_limits<double>::quiet_NaN(); | ||
| const double force_pc = (force!=NULL) ? | ||
| force->rhs[q][dim+1] | ||
| : std::numeric_limits<double>::quiet_NaN(); |
source/simulator/parameters.cc
Outdated
| "Whether to ask the material model for additional terms for the right-hand side " | ||
| "of the Stokes equation. This can be used to implement force vectors for " | ||
| "manufactured solution problems. The assembler for the Stokes system will " | ||
| "create an additional output of type AdditionalMaterialOutputsStokesRHS."); |
There was a problem hiding this comment.
the last sentence seems like a very detailed description of internals. there are many levels of knowledge to reach for a user to understand what this means.
tjhei
force-pushed
the
force_vector_the_return
branch
2 times, most recently
from
05f36fc
to
4d1c692
Compare
tjhei
force-pushed
the
force_vector_the_return
branch
from
9503d60
to
ed9d82f
Compare
|
okay, ready for another round! |
| const FullMatrix<double> &/*projection_matrix*/, | ||
| const FullMatrix<double> &/*expansion_matrix*/) | ||
| { | ||
| // TODO: not implemented |
|
W writing: After further consideration, I'm merging from @tjhei's computer. The assertion cannot be added because the function is actually called. |
Latest idea on how to do force vector for the Stokes system.
supersedes #1300