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.

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

Compositing MaterialModel #1602

Merged
merged 14 commits into from May 29, 2017
10 changes: 4 additions & 6 deletions include/aspect/material_model/compositing.h
Expand Up @@ -18,8 +18,8 @@
<http://www.gnu.org/licenses/>.
*/

#ifndef _aspect_material_model_averaging_h
#define _aspect_material_model_averaging_h
#ifndef _aspect_material_model_compositing_h
#define _aspect_material_model_compositing_h

#include <aspect/material_model/interface.h>
#include <aspect/simulator_access.h>
Expand All @@ -28,8 +28,6 @@ namespace aspect
{
namespace MaterialModel
{
using namespace dealii;

namespace Property
{
/**
Expand Down Expand Up @@ -64,7 +62,7 @@ namespace aspect
}

/**
* A material model that selects properties from other (non-averaging) material models.
* A material model that selects properties from other (non-compositing) material models.
* @ingroup MaterialModels
*/

Expand Down Expand Up @@ -119,7 +117,7 @@ namespace aspect
parse_property_name(const std::string &s);

/**
* Copy desired properties from material model,
* Copy desired properties from material model outputs produced by another material model,
*
* @param model_index Internal index for pointer to the material model evaluated
* @param base_output Properties generated by the material model specified
Expand Down
65 changes: 26 additions & 39 deletions source/material_model/compositing.cc
Expand Up @@ -18,12 +18,7 @@
<http://www.gnu.org/licenses/>.
*/

#include <deal.II/base/std_cxx11/array.h>
#include <aspect/material_model/compositing.h>
#include <utility>
#include <limits>

using namespace dealii;

namespace aspect
{
Expand All @@ -34,15 +29,12 @@ namespace aspect
Property::MaterialProperty
Compositing<dim>::parse_property_name(const std::string &s)
{
if (Property::property_map.count(s) == 0)
AssertThrow (false,
ExcMessage ("The value <" + s + "> for a material "
"property is not one of the valid values."));
return Property::viscosity;
AssertThrow (Property::property_map.count(s) == 0,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should it not be Property::property_map.count(s) > 0?

ExcMessage ("The value <" + s + "> for a material "
"property is not one of the valid values."));
return Property::property_map.at(s);
}

//Copy the requested data for one model
template <int dim>
void
Compositing<dim>::copy_required_properties(const unsigned int model_index,
Expand Down Expand Up @@ -95,8 +87,10 @@ namespace aspect
std::map<std::string, Property::MaterialProperty>::const_iterator prop_it = Property::property_map.begin();
for (; prop_it != Property::property_map.end(); ++prop_it)
{
prm.declare_entry(prop_it->first, "compositing",
Patterns::Selection(MaterialModel::get_valid_model_names_pattern<dim>()),
prm.declare_entry(prop_it->first, "unspecified",
Patterns::Selection(
MaterialModel::get_valid_model_names_pattern<dim>()+"|unspecified"
),
"Material model to use for " + prop_it->first +". Valid values for this "
"parameter are the names of models that are also valid for the "
"``Material models/Model name'' parameter. See the documentation for "
Expand Down Expand Up @@ -129,44 +123,37 @@ namespace aspect
AssertThrow(model_name != "compositing",
ExcMessage("You may not use ``compositing'' as the base model for the "
+ prop_it->first +" property of a compositing material model."));
unsigned int model_ind;
for (model_ind=0; model_ind < model_names.size(); ++model_ind)
std::vector<std::string>::iterator model_position = std::find(model_names.begin(), model_names.end(), model_name);
if ( model_position == model_names.end() )
{
if (model_names[model_ind] == model_name)
break;
model_property_map[prop] = model_names.size();
model_names.push_back(model_name);
}
if (model_ind == model_names.size())
model_names.push_back(model_name);

model_property_map[prop] = model_ind;
else
model_property_map[prop] = std::distance(model_names.begin(), model_position);
}
models.resize(model_names.size());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be great if you could write a comment about what the code above does.


// create the model and initialize their SimulatorAccess base
// class; it will get a chance to read its parameters below after we
// leave the current section
for (unsigned int i=0; i<model_names.size(); ++i)
{
models.at(i).reset(create_material_model<dim>(model_names[i]));
if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(models.at(i).get()))
sim->initialize_simulator (this->get_simulator());
}
}
prm.leave_subsection();
}
prm.leave_subsection();

/* After parsing the parameters for averaging, it is essential to parse
parameters related to the base model. */
// create the models and initialize their SimulatorAccess base
// After parsing the parameters for averaging, it is essential to parse
// parameters related to the base models
for (unsigned int i=0; i<model_names.size(); ++i)
{
models.at(i).reset(create_material_model<dim>(model_names[i]));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One more nitpick (I just mention it because I have the above comment anyway): Now you resize models above (probably for speed?), but then use at here instead of [] further down, probably just copied from the earlier loop? Two comments about that:

  1. Please move the resize command in front of this loop, where models is used, so that related commands are close together.
  2. Use consistently at or [] do not mix between them. In this case [] is totally safe, because i will never exceed the bounds of the vector. In general at is significantly slower (see e.g. https://stackoverflow.com/questions/9376049/vectorat-vs-vectoroperator), although of course in this case with just a few entries this does not matter at all.

if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(models.at(i).get()))
sim->initialize_simulator (this->get_simulator());
models[i]->parse_parameters(prm);
// All models will need to compute all quantities, so do so
this -> model_dependence.viscosity |= models[i]->get_model_dependence().viscosity;
this -> model_dependence.density |= models[i]->get_model_dependence().density;
this -> model_dependence.compressibility |= models[i]->get_model_dependence().compressibility;
this -> model_dependence.specific_heat |= models[i]->get_model_dependence().specific_heat;
this -> model_dependence.thermal_conductivity |= models[i]->get_model_dependence().thermal_conductivity;
this->model_dependence.viscosity |= models[i]->get_model_dependence().viscosity;
this->model_dependence.density |= models[i]->get_model_dependence().density;
this->model_dependence.compressibility |= models[i]->get_model_dependence().compressibility;
this->model_dependence.specific_heat |= models[i]->get_model_dependence().specific_heat;
this->model_dependence.thermal_conductivity |= models[i]->get_model_dependence().thermal_conductivity;
}
}

Expand All @@ -175,7 +162,7 @@ namespace aspect
Compositing<dim>::
is_compressible () const
{
unsigned int ind = model_property_map.at(Property::compressibility);
const unsigned int ind = model_property_map.at(Property::compressibility);
return models[ind]->is_compressible();
}

Expand All @@ -184,7 +171,7 @@ namespace aspect
Compositing<dim>::
reference_viscosity() const
{
unsigned int ind = model_property_map.at(Property::viscosity);
const unsigned int ind = model_property_map.at(Property::viscosity);
return models[ind]->reference_viscosity();
}
}
Expand Down
9 changes: 2 additions & 7 deletions tests/simple_composite.prm
@@ -1,10 +1,5 @@
# Duplicate of the simple incompressible model with using the
# composite model instead

# This is a relatively realistic compressible model setup.
# There is no simple analytical solution to it, it is rather
# intended to serve as a comparison for more complicated
# material models, which should be able to reproduce it.
# Duplicate of simple incompressible test using a compositing model to
# select all material properties from the simple incompressible model.

set CFL number = 1.0
set End time = 1e5
Expand Down