Skip to content

Commit

Permalink
account for the darcy and fluid velocities when computing outflow bou…
Browse files Browse the repository at this point in the history 
…ndaries
  • Loading branch information
@danieldouglas92
danieldouglas92 committed Oct 22, 2024
1 parent eac4de4 commit 65e611b
Show file tree
Hide file tree
Showing 3 changed files with 108 additions and 35 deletions.
4 changes: 3 additions & 1 deletion include/aspect/simulator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1545,7 +1545,9 @@ namespace aspect
* This function is implemented in
* <code>source/simulator/helper_functions.cc</code>.
*/
void replace_outflow_boundary_ids(const unsigned int boundary_id_offset);
void replace_outflow_boundary_ids(const unsigned int boundary_id_offset,
const bool is_composition,
const unsigned int composition_index);

/**
* Undo the offset of the boundary ids done in replace_outflow_boundary_ids
Expand Down
55 changes: 28 additions & 27 deletions source/simulator/core.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -695,7 +695,7 @@ namespace aspect
// so we want to offset them by 128 and not allow more than 128 boundary ids.
const unsigned int boundary_id_offset = 128;
if (!boundary_temperature_manager.allows_fixed_temperature_on_outflow_boundaries())
replace_outflow_boundary_ids(boundary_id_offset);
replace_outflow_boundary_ids(boundary_id_offset, false, 0);

// if using continuous temperature FE, do the same for the temperature variable:
// evaluate the current boundary temperature and add these constraints as well
Expand Down Expand Up @@ -730,41 +730,42 @@ namespace aspect
// update the composition boundary condition.
boundary_composition_manager.update();

// If we do not want to prescribe Dirichlet boundary conditions on outflow boundaries,
// use the same trick for marking up outflow boundary conditions for compositional fields
// as we did above already for the temperature.
if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
replace_outflow_boundary_ids(boundary_id_offset);

// now do the same for the composition variables:
{
// obtain the boundary indicators that belong to Dirichlet-type
// composition boundary conditions and interpolate the composition
// there
for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
if (parameters.use_discontinuous_composition_discretization[c] == false)
for (const auto p : boundary_composition_manager.get_fixed_composition_boundary_indicators())
{
VectorFunctionFromScalarFunctionObject<dim> vector_function_object(
[&] (const Point<dim> &x) -> double
{
// If we do not want to prescribe Dirichlet boundary conditions on outflow boundaries,
// use the same trick for marking up outflow boundary conditions for compositional fields
// as we did above already for the temperature.
if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
replace_outflow_boundary_ids(boundary_id_offset, true, c);

if (parameters.use_discontinuous_composition_discretization[c] == false)
for (const auto p : boundary_composition_manager.get_fixed_composition_boundary_indicators())
{
return boundary_composition_manager.boundary_composition(p, x, c);
},
introspection.component_masks.compositional_fields[c].first_selected_component(),
introspection.n_components);

VectorTools::interpolate_boundary_values (*mapping,
dof_handler,
p,
vector_function_object,
new_current_constraints,
introspection.component_masks.compositional_fields[c]);
}
VectorFunctionFromScalarFunctionObject<dim> vector_function_object(
[&] (const Point<dim> &x) -> double
{
return boundary_composition_manager.boundary_composition(p, x, c);
},
introspection.component_masks.compositional_fields[c].first_selected_component(),
introspection.n_components);

VectorTools::interpolate_boundary_values (*mapping,
dof_handler,
p,
vector_function_object,
new_current_constraints,
introspection.component_masks.compositional_fields[c]);
if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
restore_outflow_boundary_ids(boundary_id_offset);
}
}
}

if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
restore_outflow_boundary_ids(boundary_id_offset);

if (parameters.include_melt_transport)
melt_handler->add_current_constraints (new_current_constraints);

Expand Down
84 changes: 77 additions & 7 deletions source/simulator/helper_functions.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2271,17 +2271,53 @@ namespace aspect

template <int dim>
void
Simulator<dim>::replace_outflow_boundary_ids(const unsigned int offset)
Simulator<dim>::replace_outflow_boundary_ids(const unsigned int offset,
const bool is_composition,
const unsigned int composition_index)
{
const Quadrature<dim-1> &quadrature_formula = introspection.face_quadratures.temperature;

bool consider_darcy_velocity = false;
bool consider_melt_velocity = false;
unsigned int porosity_idx = numbers::invalid_unsigned_int;
if (is_composition)
{
if (introspection.compositional_field_methods[composition_index] == Parameters<dim>::AdvectionFieldMethod::fem_darcy_field)
{
consider_darcy_velocity = true;
porosity_idx = introspection.find_composition_type(CompositionalFieldDescription::porosity);
}
if (introspection.compositional_field_methods[composition_index] == Parameters<dim>::AdvectionFieldMethod::fem_melt_field)
consider_melt_velocity = true;
}

const UpdateFlags update_flags
= UpdateFlags(
consider_darcy_velocity
?
update_values |
update_gradients |
update_normal_vectors |
update_quadrature_points |
update_JxW_values
:
update_values |
update_normal_vectors |
update_quadrature_points |
update_JxW_values);

FEFaceValues<dim> fe_face_values (*mapping,
finite_element,
quadrature_formula,
update_values | update_normal_vectors |
update_quadrature_points | update_JxW_values);
update_flags);

std::vector<Tensor<1,dim>> face_current_velocity_values (fe_face_values.n_quadrature_points);
std::vector<Tensor<1,dim>> face_current_fluid_velocity_values (fe_face_values.n_quadrature_points);

MaterialModel::MaterialModelInputs<dim> in(fe_face_values.n_quadrature_points, introspection.n_compositional_fields);
MaterialModel::MaterialModelOutputs<dim> out(fe_face_values.n_quadrature_points, introspection.n_compositional_fields);
MeltHandler<dim>::create_material_model_outputs(out);
MaterialModel::MeltOutputs<dim> *fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();

const auto &tangential_velocity_boundaries =
boundary_velocity_manager.get_tangential_boundary_velocity_indicators();
Expand Down Expand Up @@ -2311,11 +2347,17 @@ namespace aspect
fe_face_values.reinit (cell, face_number);
fe_face_values[introspection.extractors.velocities].get_function_values(current_linearization_point,
face_current_velocity_values);

// ... check if the face is an outflow boundary by integrating the normal velocities
// (flux through the boundary) as: int u*n ds = Sum_q u(x_q)*n(x_q) JxW(x_q)...
// do this for the solid velocity, the darcy velocity, and/or the fluid velocity.
double integrated_flow = 0;

if (consider_darcy_velocity)
{
in.reinit(fe_face_values, cell, introspection, solution);
material_model->evaluate(in, out);
}

for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
{
Tensor<1,dim> boundary_velocity;
Expand All @@ -2325,8 +2367,34 @@ namespace aspect
boundary_velocity = boundary_velocity_manager.boundary_velocity(face->boundary_id(),
fe_face_values.quadrature_point(q));

integrated_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
fe_face_values.JxW(q);
if (consider_darcy_velocity)
{
const double porosity = std::max(in.composition[q][porosity_idx], 1e-10);
const Tensor<1,dim> gravity = gravity_model->gravity_vector(in.position[q]);
const double solid_density = out.densities[q];
const double fluid_viscosity = fluid_out->fluid_viscosities[q];
const double fluid_density = fluid_out->fluid_densities[q];
const double permeability = fluid_out->permeabilities[q];
const Tensor<1,dim> boundary_darcy_velocity = boundary_velocity -
permeability / fluid_viscosity / porosity * gravity *
(solid_density - fluid_density);
integrated_flow += (boundary_darcy_velocity * fe_face_values.normal_vector(q)) *
fe_face_values.JxW(q);
}

else if (consider_melt_velocity)
{
const FEValuesExtractors::Vector ex_u_f = introspection.variable("fluid velocity").extractor_vector();
fe_face_values[ex_u_f].get_function_values(current_linearization_point, face_current_fluid_velocity_values);
integrated_flow += (face_current_fluid_velocity_values[q] * fe_face_values.normal_vector(q)) *
fe_face_values.JxW(q);
}

else
{
integrated_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
fe_face_values.JxW(q);
}
}

// ... and change the boundary id of any outflow boundary faces.
Expand Down Expand Up @@ -2758,7 +2826,9 @@ namespace aspect
template void Simulator<dim>::initialize_current_linearization_point (); \
template void Simulator<dim>::interpolate_material_output_into_advection_field(const AdvectionField &adv_field); \
template void Simulator<dim>::check_consistency_of_formulation(); \
template void Simulator<dim>::replace_outflow_boundary_ids(const unsigned int boundary_id_offset); \
template void Simulator<dim>::replace_outflow_boundary_ids(const unsigned int boundary_id_offset, \
const bool is_composition, \
const unsigned int composition_index); \
template void Simulator<dim>::restore_outflow_boundary_ids(const unsigned int boundary_id_offset); \
template void Simulator<dim>::check_consistency_of_boundary_conditions() const; \
template double Simulator<dim>::compute_initial_newton_residual(const LinearAlgebra::BlockVector &linearized_stokes_initial_guess); \
Expand Down

0 comments on commit 65e611b

Please sign in to comment.