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

Account for the Darcy Velocity and/or the Fluid Velocity when Determining Outflow Boundaries #6075

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
account for the darcy and fluid velocities when computing outflow bou…
…ndaries
danieldouglas92 committed Nov 4, 2024
commit 6554d2e2321d1d38888993d3831df3a82d42d63d
4 changes: 3 additions & 1 deletion include/aspect/simulator.h
Original file line number Diff line number Diff line change
@@ -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
56 changes: 29 additions & 27 deletions source/simulator/core.cc
Original file line number Diff line number Diff line change
@@ -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, numbers::invalid_unsigned_int);

// if using continuous temperature FE, do the same for the temperature variable:
// evaluate the current boundary temperature and add these constraints as well
@@ -730,40 +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);
95 changes: 89 additions & 6 deletions source/simulator/helper_functions.cc
Original file line number Diff line number Diff line change
@@ -2277,17 +2277,61 @@ 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_fluid_velocity = false;
unsigned int porosity_idx = numbers::invalid_unsigned_int;

// Check to see if we are attempting to replace the boundary conditions for the temperature, or for the
// composition. If we are attempting to replace the composition boundary conditions, check to see if the current
// compositional field is advected with the darcy velocity (via the darcy field advection method) or with the
// fluid velocity (via the melt field advection method, or via a compositional field named porosity when melt
// transport is included).
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 ||
(parameters.include_melt_transport && introspection.get_indices_for_fields_of_type(CompositionalFieldDescription::porosity)[0] == composition_index) )
consider_fluid_velocity = true;
}

// If the compositional field uses the darcy velocity, we need to update the gradients, otherwise we do not
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();
@@ -2318,8 +2362,21 @@ namespace aspect
fe_face_values[introspection.extractors.velocities].get_function_values(current_linearization_point,
face_current_velocity_values);

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

if (consider_fluid_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);
}

// ... 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, or the fluid velocity.
double integrated_flow = 0;

for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
@@ -2331,8 +2388,32 @@ 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_fluid_velocity)
{
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.
@@ -2764,7 +2845,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); \