From 65e611b52fcd5b167b83996412316fd2a07453f8 Mon Sep 17 00:00:00 2001 From: danieldouglas92 Date: Thu, 5 Sep 2024 13:16:23 -0600 Subject: [PATCH] account for the darcy and fluid velocities when computing outflow boundaries --- include/aspect/simulator.h | 4 +- source/simulator/core.cc | 55 +++++++++--------- source/simulator/helper_functions.cc | 84 +++++++++++++++++++++++++--- 3 files changed, 108 insertions(+), 35 deletions(-) diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h index 4d5c99386d4..ece36e11799 100644 --- a/include/aspect/simulator.h +++ b/include/aspect/simulator.h @@ -1545,7 +1545,9 @@ namespace aspect * This function is implemented in * source/simulator/helper_functions.cc. */ - 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 diff --git a/source/simulator/core.cc b/source/simulator/core.cc index e8173dae562..e3907a11bc6 100644 --- a/source/simulator/core.cc +++ b/source/simulator/core.cc @@ -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 @@ -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 vector_function_object( - [&] (const Point &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 vector_function_object( + [&] (const Point &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); diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc index 8a5a3742909..4a60b980f80 100644 --- a/source/simulator/helper_functions.cc +++ b/source/simulator/helper_functions.cc @@ -2271,17 +2271,53 @@ namespace aspect template void - Simulator::replace_outflow_boundary_ids(const unsigned int offset) + Simulator::replace_outflow_boundary_ids(const unsigned int offset, + const bool is_composition, + const unsigned int composition_index) { const Quadrature &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::AdvectionFieldMethod::fem_darcy_field) + { + consider_darcy_velocity = true; + porosity_idx = introspection.find_composition_type(CompositionalFieldDescription::porosity); + } + if (introspection.compositional_field_methods[composition_index] == Parameters::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 fe_face_values (*mapping, finite_element, quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | update_JxW_values); + update_flags); std::vector> face_current_velocity_values (fe_face_values.n_quadrature_points); + std::vector> face_current_fluid_velocity_values (fe_face_values.n_quadrature_points); + + MaterialModel::MaterialModelInputs in(fe_face_values.n_quadrature_points, introspection.n_compositional_fields); + MaterialModel::MaterialModelOutputs out(fe_face_values.n_quadrature_points, introspection.n_compositional_fields); + MeltHandler::create_material_model_outputs(out); + MaterialModel::MeltOutputs *fluid_out = out.template get_additional_output>(); const auto &tangential_velocity_boundaries = boundary_velocity_manager.get_tangential_boundary_velocity_indicators(); @@ -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 boundary_velocity; @@ -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. @@ -2758,7 +2826,9 @@ namespace aspect template void Simulator::initialize_current_linearization_point (); \ template void Simulator::interpolate_material_output_into_advection_field(const AdvectionField &adv_field); \ template void Simulator::check_consistency_of_formulation(); \ - template void Simulator::replace_outflow_boundary_ids(const unsigned int boundary_id_offset); \ + template void Simulator::replace_outflow_boundary_ids(const unsigned int boundary_id_offset, \ + const bool is_composition, \ + const unsigned int composition_index); \ template void Simulator::restore_outflow_boundary_ids(const unsigned int boundary_id_offset); \ template void Simulator::check_consistency_of_boundary_conditions() const; \ template double Simulator::compute_initial_newton_residual(const LinearAlgebra::BlockVector &linearized_stokes_initial_guess); \