New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Remove code duplication in initialization of temperature and composition #1154
Changes from 4 commits
c6366ce
3cc745a
209b56b
06e3c39
a636ee7
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -51,130 +51,8 @@ namespace aspect | |
// we need to track whether we need to normalize the totality of fields | ||
bool normalize_composition = false; | ||
|
||
#if DEAL_II_VERSION_GTE(8,5,0) | ||
initial_solution.reinit(system_rhs, false); | ||
|
||
// Interpolate initial temperature | ||
VectorFunctionFromScalarFunctionObject<dim, double> | ||
temperature_function(std_cxx11::bind(&InitialConditions::Interface<dim>::initial_temperature, | ||
std::ref(*initial_conditions), | ||
std_cxx11::_1), | ||
introspection.component_indices.temperature, | ||
introspection.n_components); | ||
|
||
VectorTools::interpolate(*mapping, | ||
dof_handler, | ||
temperature_function, | ||
initial_solution, | ||
introspection.component_masks.temperature); | ||
|
||
// then apply constraints and copy the | ||
// result into vectors with ghost elements. to do so, | ||
// we need the current constraints to be correct for | ||
// the current time | ||
compute_current_constraints (); | ||
current_constraints.distribute(initial_solution); | ||
|
||
// copy temperature/composition block only | ||
const unsigned int t_blockidx = introspection.block_indices.temperature; | ||
solution.block(t_blockidx) = initial_solution.block(t_blockidx); | ||
old_solution.block(t_blockidx) = initial_solution.block(t_blockidx); | ||
old_old_solution.block(t_blockidx) = initial_solution.block(t_blockidx); | ||
|
||
// Interpolate initial compositions | ||
for (unsigned int n=0; n<parameters.n_compositional_fields; ++n) | ||
{ | ||
AdvectionField advf = AdvectionField::composition(n); | ||
initial_solution.reinit(system_rhs, false); | ||
|
||
VectorFunctionFromScalarFunctionObject<dim, double> | ||
composition_function (std_cxx11::bind(&CompositionalInitialConditions::Interface<dim>::initial_composition, | ||
std::ref(*compositional_initial_conditions), | ||
std_cxx11::_1, | ||
n), | ||
introspection.component_indices.compositional_fields[n], | ||
introspection.n_components); | ||
|
||
|
||
VectorTools::interpolate(*mapping, | ||
dof_handler, | ||
composition_function, | ||
initial_solution, | ||
introspection.component_masks.compositional_fields[n]); | ||
|
||
const unsigned int base_element = advf.base_element(introspection); | ||
|
||
// get the temperature/composition support points | ||
const std::vector<Point<dim> > support_points | ||
= finite_element.base_element(base_element).get_unit_support_points(); | ||
Assert (support_points.size() != 0, | ||
ExcInternalError()); | ||
|
||
// create an FEValues object to check compositions at support points | ||
FEValues<dim> fe_values (*mapping, finite_element, | ||
support_points, | ||
update_quadrature_points); | ||
|
||
std::vector<types::global_dof_index> local_dof_indices (finite_element.dofs_per_cell); | ||
|
||
for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(); | ||
cell != dof_handler.end(); ++cell) | ||
if (cell->is_locally_owned()) | ||
{ | ||
fe_values.reinit (cell); | ||
|
||
// Go through the support points for each dof | ||
for (unsigned int i=0; i<finite_element.base_element(base_element).dofs_per_cell; ++i) | ||
{ | ||
// if it is specified in the parameter file that the sum of all compositional fields | ||
// must not exceed one, this should be checked | ||
if (parameters.normalized_fields.size()>0 && n == 1) | ||
{ | ||
double sum = 0; | ||
for (unsigned int m=0; m<parameters.normalized_fields.size(); ++m) | ||
sum += compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i),parameters.normalized_fields[m]); | ||
if (std::abs(sum) > 1.0+1e-6) | ||
{ | ||
max_sum_comp = std::max(sum, max_sum_comp); | ||
normalize_composition = true; | ||
} | ||
} | ||
|
||
} | ||
} | ||
|
||
// if at least one processor decides that it needs | ||
// to normalize, do the same on all processors. | ||
if (Utilities::MPI::max (normalize_composition ? 1 : 0, | ||
mpi_communicator) | ||
== 1) | ||
{ | ||
const double global_max | ||
= Utilities::MPI::max (max_sum_comp, mpi_communicator); | ||
|
||
if (n==1) | ||
pcout << "Sum of compositional fields is not one, fields will be normalized" | ||
<< std::endl; | ||
|
||
for (unsigned int m=0; m<parameters.normalized_fields.size(); ++m) | ||
if (n==parameters.normalized_fields[m]) | ||
initial_solution /= global_max; | ||
} | ||
|
||
// then apply constraints and copy the | ||
// result into vectors with ghost elements. to do so, | ||
// we need the current constraints to be correct for | ||
// the current time | ||
compute_current_constraints (); | ||
current_constraints.distribute(initial_solution); | ||
|
||
// copy temperature/composition block only | ||
const unsigned int blockidx = introspection.block_indices.compositional_fields[n]; | ||
solution.block(blockidx) = initial_solution.block(blockidx); | ||
old_solution.block(blockidx) = initial_solution.block(blockidx); | ||
old_old_solution.block(blockidx) = initial_solution.block(blockidx); | ||
} | ||
#else | ||
// below, we would want to call VectorTools::interpolate on the | ||
// entire FESystem. there currently is no way to restrict the | ||
// interpolation operations to only a subset of vector | ||
|
@@ -197,7 +75,6 @@ namespace aspect | |
{ | ||
AdvectionField advf = ((n == 0) ? AdvectionField::temperature() | ||
: AdvectionField::composition(n-1)); | ||
initial_solution.reinit(system_rhs, false); | ||
|
||
const unsigned int base_element = advf.base_element(introspection); | ||
|
||
|
@@ -214,6 +91,55 @@ namespace aspect | |
|
||
std::vector<types::global_dof_index> local_dof_indices (finite_element.dofs_per_cell); | ||
|
||
#if DEAL_II_VERSION_GTE(8,5,0) | ||
const VectorFunctionFromScalarFunctionObject<dim, double> &advf_init_function = | ||
(advf.is_temperature()? | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Please put the ? and : operators on separate lines. |
||
VectorFunctionFromScalarFunctionObject<dim, double>(std_cxx11::bind(&InitialConditions::Interface<dim>::initial_temperature, | ||
std::ref(*initial_conditions), | ||
std_cxx11::_1), | ||
introspection.component_indices.temperature, | ||
introspection.n_components): | ||
VectorFunctionFromScalarFunctionObject<dim, double>(std_cxx11::bind(&CompositionalInitialConditions::Interface<dim>::initial_composition, | ||
std::ref(*compositional_initial_conditions), | ||
std_cxx11::_1, | ||
n-1), | ||
introspection.component_indices.compositional_fields[n-1], | ||
introspection.n_components)); | ||
|
||
const ComponentMask advf_mask = | ||
(advf.is_temperature()? | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Same here |
||
introspection.component_masks.temperature: | ||
introspection.component_masks.compositional_fields[n-1]); | ||
|
||
VectorTools::interpolate(*mapping, | ||
dof_handler, | ||
advf_init_function, | ||
initial_solution, | ||
advf_mask); | ||
|
||
if (n == 1 && parameters.normalized_fields.size()>0) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. In the other code path (for deal.II 8.4.0) the condition is written the other way |
||
for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(); | ||
cell != dof_handler.end(); ++cell) | ||
if (cell->is_locally_owned()) | ||
{ | ||
fe_values.reinit (cell); | ||
|
||
// Go through the support points for each dof | ||
for (unsigned int i=0; i<finite_element.base_element(base_element).dofs_per_cell; ++i) | ||
{ | ||
// if it is specified in the parameter file that the sum of all compositional fields | ||
// must not exceed one, this should be checked | ||
double sum = 0; | ||
for (unsigned int m=0; m<parameters.normalized_fields.size(); ++m) | ||
sum += compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i),parameters.normalized_fields[m]); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Please break this line after the first parameter of initial_composition (simpler reading). |
||
if (std::abs(sum) > 1.0+1e-6) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think we could replace 1e-6 by There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. And there is another occurence in the other code path, please replace that as well. I realize it was there before your change, but this is the opportunity to clean up 😄 |
||
{ | ||
max_sum_comp = std::max(sum, max_sum_comp); | ||
normalize_composition = true; | ||
} | ||
} | ||
} | ||
#else | ||
for (typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(); | ||
cell != dof_handler.end(); ++cell) | ||
if (cell->is_locally_owned()) | ||
|
@@ -253,16 +179,10 @@ namespace aspect | |
|
||
} | ||
} | ||
#endif | ||
|
||
initial_solution.compress(VectorOperation::insert); | ||
|
||
// we should not have written at all into any of the blocks with | ||
// the exception of the current temperature or composition block | ||
for (unsigned int b=0; b<initial_solution.n_blocks(); ++b) | ||
if (b != advf.block_index(introspection)) | ||
Assert (initial_solution.block(b).l2_norm() == 0, | ||
ExcInternalError()); | ||
|
||
// if at least one processor decides that it needs | ||
// to normalize, do the same on all processors. | ||
if (Utilities::MPI::max (normalize_composition ? 1 : 0, | ||
|
@@ -278,23 +198,30 @@ namespace aspect | |
|
||
for (unsigned int m=0; m<parameters.normalized_fields.size(); ++m) | ||
if (n-1==parameters.normalized_fields[m]) | ||
initial_solution /= global_max; | ||
initial_solution.block(introspection.block_indices.compositional_fields[n-1]) /= global_max; | ||
} | ||
} | ||
|
||
// then apply constraints and copy the | ||
// result into vectors with ghost elements. to do so, | ||
// we need the current constraints to be correct for | ||
// the current time | ||
compute_current_constraints (); | ||
current_constraints.distribute(initial_solution); | ||
// then apply constraints and copy the | ||
// result into vectors with ghost elements. to do so, | ||
// we need the current constraints to be correct for | ||
// the current time | ||
compute_current_constraints (); | ||
current_constraints.distribute(initial_solution); | ||
|
||
// Now copy the temperature and initial composition blocks into the solution variables | ||
|
||
for (unsigned int n=0; n<1+parameters.n_compositional_fields; ++n) | ||
{ | ||
AdvectionField advf = ((n == 0) ? AdvectionField::temperature() | ||
: AdvectionField::composition(n-1)); | ||
|
||
// copy temperature/composition block only | ||
const unsigned int blockidx = advf.block_index(introspection); | ||
|
||
solution.block(blockidx) = initial_solution.block(blockidx); | ||
old_solution.block(blockidx) = initial_solution.block(blockidx); | ||
old_old_solution.block(blockidx) = initial_solution.block(blockidx); | ||
} | ||
#endif | ||
} | ||
|
||
|
||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,61 @@ | ||
# This test ensures that the Simpler material model works | ||
# with compositional fields enabled | ||
|
||
set Dimension = 2 | ||
set End time = 0 | ||
|
||
subsection Boundary temperature model | ||
set Model name = initial temperature | ||
end | ||
|
||
subsection Compositional fields | ||
set Number of fields = 2 | ||
set List of normalized fields = 0,1 | ||
end | ||
|
||
subsection Compositional initial conditions | ||
subsection Function | ||
set Function expression = 1;1 | ||
end | ||
end | ||
|
||
subsection Geometry model | ||
set Model name = box | ||
subsection Box | ||
set X extent = 1 | ||
set X repetitions = 2 | ||
set Y extent = 1 | ||
set Y repetitions = 2 | ||
end | ||
end | ||
|
||
subsection Gravity model | ||
set Model name = vertical | ||
subsection Vertical | ||
set Magnitude = 1 | ||
end | ||
end | ||
|
||
subsection Initial conditions | ||
set Model name = function | ||
subsection Function | ||
set Function expression = 1 | ||
end | ||
end | ||
|
||
subsection Material model | ||
set Model name = simple | ||
end | ||
|
||
subsection Mesh refinement | ||
set Initial adaptive refinement = 0 | ||
set Initial global refinement = 0 | ||
end | ||
|
||
subsection Model settings | ||
set Tangential velocity boundary indicators = left, right, top, bottom | ||
end | ||
|
||
subsection Postprocess | ||
set List of postprocessors = composition statistics | ||
end |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,24 @@ | ||
----------------------------------------------------------------------------- | ||
----------------------------------------------------------------------------- | ||
|
||
Number of active cells: 4 (on 1 levels) | ||
Number of degrees of freedom: 134 (50+9+25+25+25) | ||
|
||
Sum of compositional fields is not one, fields will be normalized | ||
*** Timestep 0: t=0 years | ||
Solving temperature system... 0 iterations. | ||
Solving C_1 system ... 0 iterations. | ||
Solving C_2 system ... 0 iterations. | ||
Rebuilding Stokes preconditioner... | ||
Solving Stokes system... 2+0 iterations. | ||
|
||
Postprocessing: | ||
Compositions min/max/mass: 0.5/0.5/0.5 // 0.5/0.5/0.5 | ||
|
||
Termination requested by criterion: end time | ||
|
||
|
||
+---------------------------------------------+------------+------------+ | ||
+---------------------------------+-----------+------------+------------+ | ||
+---------------------------------+-----------+------------+------------+ | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The tester currently only tests deal.II 8.4.0. Did you test that this changed code works with 8.5.0.pre? If you do not have a development version of deal.II I could also do that for you. Let me know.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I did check that the changed code worked for 8.5.0-pre. I will check again with the most recent version of deal.II though.