From b52900aa76a01449da3931c2d9db9adef5da7f08 Mon Sep 17 00:00:00 2001
From: danieldouglas92 <daniel.douglas@student.nmt.edu>
Date: Mon, 14 Oct 2024 16:23:59 -0600
Subject: [PATCH] account for the fluid velocity for outflow boundaries

---
 source/simulator/helper_functions.cc | 23 +++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index 55ee60b6c83..2510ee8f6de 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -2291,12 +2291,12 @@ namespace aspect
                                       update_quadrature_points | update_JxW_values);
 
     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);
+    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 = nullptr;
+    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();
@@ -2329,8 +2329,11 @@ namespace aspect
 
                 // ... 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)...
-                double integrated_flow = 0;
+                // do this for the solid velocity, the darcy velocity, and/or the fluid velocity
+                // if applicable.
+                double integrated_solid_flow = 0;
                 double integrated_darcy_flow = 0;
+                double integrated_fluid_flow = 0;
 
                 for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
                   {
@@ -2356,12 +2359,20 @@ namespace aspect
                                                  fe_face_values.JxW(q);
                       }
 
-                    integrated_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
-                                       fe_face_values.JxW(q);
+                    if (parameters.include_melt_transport == true)
+                      {
+                        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_fluid_flow += (face_current_fluid_velocity_values[q] * fe_face_values.normal_vector(q)) *
+                                                 fe_face_values.JxW(q);
+                      }
+
+                    integrated_solid_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
+                                             fe_face_values.JxW(q);
                   }
 
                 // ... and change the boundary id of any outflow boundary faces.
-                if (integrated_flow > 0 || integrated_darcy_flow > 0)
+                if (integrated_solid_flow > 0 || integrated_darcy_flow > 0 || integrated_fluid_flow > 0)
                   face->set_boundary_id(face->boundary_id() + offset);
               }
           }