Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Remove code duplication in initialization of temperature and composition #1154

Merged
merged 5 commits into from Aug 12, 2016
Merged

Remove code duplication in initialization of temperature and composition #1154

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
c6366ce
Remove duplication of constraint computation in initialization
class4kayaker Aug 3, 2016
3cc745a
Add test for composition normalization
class4kayaker Aug 3, 2016
209b56b
Remove constraint/copy duplication
class4kayaker Aug 3, 2016
06e3c39
Collapse redundancies
class4kayaker Aug 3, 2016
a636ee7
Tidy up modification to initial_conditions.cc
class4kayaker Aug 9, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
215 changes: 74 additions & 141 deletions source/simulator/initial_conditions.cc
View file
Open in desktop
Expand Up @@ -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
Expand All @@ -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);

Expand All @@ -214,6 +91,60 @@ namespace aspect

std::vector<types::global_dof_index> local_dof_indices (finite_element.dofs_per_cell);

#if DEAL_II_VERSION_GTE(8,5,0)
Copy link
Member

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.

Copy link
Contributor Author

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.

const VectorFunctionFromScalarFunctionObject<dim, double> &advf_init_function =
(advf.is_temperature()
?
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()
?
introspection.component_masks.temperature
:
introspection.component_masks.compositional_fields[n-1]);

VectorTools::interpolate(*mapping,
dof_handler,
advf_init_function,
initial_solution,
advf_mask);

if (parameters.normalized_fields.size()>0 && n==1)
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]);
if (std::abs(sum) > 1.0+std::numeric_limits<double>::epsilon())
{
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())
Expand Down Expand Up @@ -243,8 +174,9 @@ namespace aspect
{
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)
sum += compositional_initial_conditions->initial_composition(fe_values.quadrature_point(i),
parameters.normalized_fields[m]);
if (std::abs(sum) > 1.0+std::numeric_limits<double>::epsilon())
{
max_sum_comp = std::max(sum, max_sum_comp);
normalize_composition = true;
Expand All @@ -253,16 +185,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,
Expand All @@ -278,23 +204,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
}


Expand Down
61 changes: 61 additions & 0 deletions tests/normalize_initial_composition.prm
View file
Open in desktop
@@ -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
24 changes: 24 additions & 0 deletions tests/normalize_initial_composition/screen-output
View file
Open in desktop
@@ -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


+---------------------------------------------+------------+------------+
+---------------------------------+-----------+------------+------------+
+---------------------------------+-----------+------------+------------+