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

Fix n_comp_fields #1568

Merged
merged 2 commits into from
May 10, 2017
Merged

Fix n_comp_fields #1568

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
18fd365
Alter references to n_compositional_fields to ref introspection
class4kayaker May 10, 2017
69af917
Change to reference n_compos_fields using SimulatorAccess
class4kayaker May 10, 2017
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
4 changes: 2 additions & 2 deletions source/postprocess/global_statistics.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -241,9 +241,9 @@ namespace aspect
statistics.add_value("Number of Stokes degrees of freedom", n_stokes_dofs);
statistics.add_value("Number of temperature degrees of freedom",
this->introspection().system_dofs_per_block[this->introspection().block_indices.temperature]);
if (this->get_parameters().n_compositional_fields > 0)
if (this->n_compositional_fields() > 0)
statistics.add_value("Number of degrees of freedom for all compositions",
this->get_parameters().n_compositional_fields
this->n_compositional_fields()
* this->introspection().system_dofs_per_block[this->introspection().block_indices.compositional_fields[0]]);
}

Expand Down
18 changes: 9 additions & 9 deletions source/simulator/assembly.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -480,17 +480,17 @@ namespace aspect

// the values of the compositional fields are stored as blockvectors for each field
// we have to extract them in this structure
std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,
std::vector<std::vector<double> > composition_values (introspection.n_compositional_fields,
std::vector<double> (n_q_points));

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
input_finite_element_values[introspection.extractors.compositional_fields[c]].get_function_values(input_solution,
composition_values[c]);

// then we copy these values to exchange the inner and outer vector, because for the material
// model we need a vector with values of all the compositional fields for every quadrature point
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
material_model_inputs.composition[q][c] = composition_values[c][q];

material_model_inputs.cell = &cell;
Expand Down Expand Up @@ -614,8 +614,8 @@ namespace aspect
aspect::Assemblers::AdvectionAssembler<dim> *adv_assembler
= new aspect::Assemblers::AdvectionAssembler<dim>();

assemblers->advection_system_assembler_properties.resize(1+parameters.n_compositional_fields);
assemblers->advection_system_assembler_on_face_properties.resize(1+parameters.n_compositional_fields);
assemblers->advection_system_assembler_properties.resize(1+introspection.n_compositional_fields);
assemblers->advection_system_assembler_on_face_properties.resize(1+introspection.n_compositional_fields);

aspect::Assemblers::MeltEquations<dim> *melt_equation_assembler = NULL;
if (parameters.include_melt_transport)
Expand Down Expand Up @@ -835,7 +835,7 @@ namespace aspect

if (parameters.use_discontinuous_composition_discretization)
{
for (unsigned int i = 1; i<=parameters.n_compositional_fields; ++i)
for (unsigned int i = 1; i<=introspection.n_compositional_fields; ++i)
{
assemblers->advection_system_assembler_on_face_properties[i].need_face_material_model_data = true;
assemblers->advection_system_assembler_on_face_properties[i].need_face_finite_element_evaluation = true;
Expand Down Expand Up @@ -963,7 +963,7 @@ namespace aspect
StokesPreconditioner<dim> (finite_element, quadrature_formula,
*mapping,
cell_update_flags,
parameters.n_compositional_fields,
introspection.n_compositional_fields,
stokes_dofs_per_cell,
parameters.include_melt_transport),
internal::Assembly::CopyData::
Expand Down Expand Up @@ -1267,7 +1267,7 @@ namespace aspect
face_quadrature_formula,
cell_update_flags,
face_update_flags,
parameters.n_compositional_fields,
introspection.n_compositional_fields,
stokes_dofs_per_cell,
parameters.include_melt_transport,
use_reference_density_profile),
Expand Down Expand Up @@ -1633,7 +1633,7 @@ namespace aspect
Quadrature<dim-1> ()),
update_flags,
face_update_flags,
parameters.n_compositional_fields),
introspection.n_compositional_fields),
internal::Assembly::CopyData::
AdvectionSystem<dim> (finite_element.base_element(advection_field.base_element(introspection)),
allocate_neighbor_contributions));
Expand Down
34 changes: 17 additions & 17 deletions source/simulator/core.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -955,7 +955,7 @@ namespace aspect
// obtain the boundary indicators that belong to Dirichlet-type
// composition boundary conditions and interpolate the composition
// there
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
for (std::set<types::boundary_id>::const_iterator
p = parameters.fixed_composition_boundary_indicators.begin();
p != parameters.fixed_composition_boundary_indicators.end(); ++p)
Expand Down Expand Up @@ -1050,7 +1050,7 @@ namespace aspect
}
}
coupling[x.temperature][x.temperature] = DoFTools::always;
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
const typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand Down Expand Up @@ -1095,7 +1095,7 @@ namespace aspect

if (parameters.use_discontinuous_composition_discretization)
{
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
const typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand Down Expand Up @@ -1202,7 +1202,7 @@ namespace aspect
coupling[x.pressure][x.pressure] = DoFTools::always;
coupling[x.temperature][x.temperature] = DoFTools::always;

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
const typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand Down Expand Up @@ -1247,7 +1247,7 @@ namespace aspect

if (parameters.use_discontinuous_composition_discretization)
{
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
const typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand Down Expand Up @@ -1956,7 +1956,7 @@ namespace aspect
current_linearization_point.block(introspection.block_indices.temperature)
= solution.block(introspection.block_indices.temperature);

for (unsigned int c=0; c < parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c < introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
const typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand All @@ -1976,7 +1976,7 @@ namespace aspect
}
}

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
current_linearization_point.block(introspection.block_indices.compositional_fields[c])
= solution.block(introspection.block_indices.compositional_fields[c]);

Expand Down Expand Up @@ -2062,7 +2062,7 @@ namespace aspect
{
double initial_temperature_residual = 0;
double initial_stokes_residual = 0;
std::vector<double> initial_composition_residual (parameters.n_compositional_fields,0);
std::vector<double> initial_composition_residual (introspection.n_compositional_fields,0);

do
{
Expand All @@ -2081,9 +2081,9 @@ namespace aspect
current_linearization_point.block(introspection.block_indices.temperature)
= solution.block(introspection.block_indices.temperature);
rebuild_stokes_matrix = true;
std::vector<double> relative_composition_residual (parameters.n_compositional_fields,0);
std::vector<double> relative_composition_residual (introspection.n_compositional_fields,0);

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand Down Expand Up @@ -2118,7 +2118,7 @@ namespace aspect

// for consistency we update the current linearization point only after we have solved
// all fields, so that we use the same point in time for every field when solving
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
current_linearization_point.block(introspection.block_indices.compositional_fields[c])
= solution.block(introspection.block_indices.compositional_fields[c]);

Expand All @@ -2144,13 +2144,13 @@ namespace aspect
// write the residual output in the same order as the output when
// solving the equations, output only relative residuals
pcout << " Relative nonlinear residuals: " << relative_temperature_residual;
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
pcout << ", " << relative_composition_residual[c];
pcout << ", " << relative_stokes_residual;
pcout << std::endl;

double max = 0.0;
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
// in models with melt migration the melt advection equation includes the divergence of the velocity
// and can not be expected to converge to a smaller value than the residual of the Stokes equation.
Expand Down Expand Up @@ -2201,7 +2201,7 @@ namespace aspect
current_linearization_point.block(introspection.block_indices.temperature)
= solution.block(introspection.block_indices.temperature);

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand All @@ -2221,7 +2221,7 @@ namespace aspect
}
}

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
current_linearization_point.block(introspection.block_indices.compositional_fields[c])
= solution.block(introspection.block_indices.compositional_fields[c]);

Expand Down Expand Up @@ -2311,7 +2311,7 @@ namespace aspect
current_linearization_point.block(introspection.block_indices.temperature)
= solution.block(introspection.block_indices.temperature);

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
const AdvectionField adv_field (AdvectionField::composition(c));
typename Parameters<dim>::AdvectionFieldMethod::Kind method = adv_field.advection_method(introspection);
Expand All @@ -2331,7 +2331,7 @@ namespace aspect
}
}

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
current_linearization_point.block(introspection.block_indices.compositional_fields[c])
= solution.block(introspection.block_indices.compositional_fields[c]);

Expand Down
8 changes: 4 additions & 4 deletions source/simulator/entropy_viscosity.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,7 @@ namespace aspect
Assert(viscosity_per_cell.size()==triangulation.n_active_cells(), ExcInternalError());

if (advection_field.field_type == AdvectionField::compositional_field)
Assert(parameters.n_compositional_fields > advection_field.compositional_variable, ExcInternalError());
Assert(introspection.n_compositional_fields > advection_field.compositional_variable, ExcInternalError());

viscosity_per_cell = 0.0;

Expand Down Expand Up @@ -255,7 +255,7 @@ namespace aspect
Quadrature<dim-1> (),
update_flags,
face_update_flags,
parameters.n_compositional_fields);
introspection.n_compositional_fields);

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active();
for (; cell<dof_handler.end(); ++cell)
Expand Down Expand Up @@ -301,7 +301,7 @@ namespace aspect
scratch.finite_element_values[introspection.extractors.pressure].get_function_values (old_old_solution,
scratch.old_old_pressure);

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
{
scratch.finite_element_values[introspection.extractors.compositional_fields[c]].get_function_values(old_solution,
scratch.old_composition_values[c]);
Expand Down Expand Up @@ -363,7 +363,7 @@ namespace aspect
scratch.material_model_inputs.velocity[q] = (scratch.old_velocity_values[q] + scratch.old_old_velocity_values[q]) / 2;
scratch.material_model_inputs.pressure_gradient[q] = (scratch.old_pressure_gradients[q] + scratch.old_old_pressure_gradients[q]) / 2;

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
scratch.material_model_inputs.composition[q][c] = (scratch.old_composition_values[c][q] + scratch.old_old_composition_values[c][q]) / 2;
scratch.material_model_inputs.strain_rate[q] = (scratch.old_strain_rates[q] + scratch.old_old_strain_rates[q]) / 2;
}
Expand Down
10 changes: 5 additions & 5 deletions source/simulator/helper_functions.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -435,7 +435,7 @@ namespace aspect

std::vector<Tensor<1,dim> > velocity_values(n_q_points);
std::vector<Tensor<1,dim> > fluid_velocity_values(n_q_points);
std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
std::vector<std::vector<double> > composition_values (introspection.n_compositional_fields,std::vector<double> (n_q_points));

double max_local_speed_over_meshsize = 0;
double min_local_conduction_timestep = std::numeric_limits<double>::max();
Expand Down Expand Up @@ -472,8 +472,8 @@ namespace aspect

if (parameters.use_conduction_timestep)
{
MaterialModel::MaterialModelInputs<dim> in(n_q_points, parameters.n_compositional_fields);
MaterialModel::MaterialModelOutputs<dim> out(n_q_points, parameters.n_compositional_fields);
MaterialModel::MaterialModelInputs<dim> in(n_q_points, introspection.n_compositional_fields);
MaterialModel::MaterialModelOutputs<dim> out(n_q_points, introspection.n_compositional_fields);

in.strain_rate.resize(0); // we do not need the viscosity
in.position = fe_values.get_quadrature_points();
Expand All @@ -486,13 +486,13 @@ namespace aspect
fe_values[introspection.extractors.pressure].get_function_gradients (solution,
in.pressure_gradient);

for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
fe_values[introspection.extractors.compositional_fields[c]].get_function_values (solution,
composition_values[c]);

for (unsigned int q=0; q<fe_values.n_quadrature_points; ++q)
{
for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
in.composition[q][c] = composition_values[c][q];

in.velocity[q] = velocity_values[q];
Expand Down
4 changes: 2 additions & 2 deletions source/simulator/initial_conditions.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ namespace aspect
//
//TODO: it would be great if we had a cleaner way than iterating to 1+n_fields.
// Additionally, the n==1 logic for normalization at the bottom is not pretty.
for (unsigned int n=0; n<1+parameters.n_compositional_fields; ++n)
for (unsigned int n=0; n<1+introspection.n_compositional_fields; ++n)
{
AdvectionField advf = ((n == 0) ? AdvectionField::temperature()
: AdvectionField::composition(n-1));
Expand Down Expand Up @@ -220,7 +220,7 @@ namespace aspect

// Now copy the temperature and initial composition blocks into the solution variables

for (unsigned int n=0; n<1+parameters.n_compositional_fields; ++n)
for (unsigned int n=0; n<1+introspection.n_compositional_fields; ++n)
{
AdvectionField advf = ((n == 0) ? AdvectionField::temperature()
: AdvectionField::composition(n-1));
Expand Down