add test case

source/simulator/helper_functions.cc

@@ -2325,7 +2325,7 @@ namespace aspect
     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>>();
+    MaterialModel::MeltOutputs<dim> *fluid_out = nullptr;
 
     const auto &tangential_velocity_boundaries =
       boundary_velocity_manager.get_tangential_boundary_velocity_indicators();
@@ -2360,12 +2360,13 @@ namespace aspect
                   {
                     in.reinit(fe_face_values, cell, introspection, solution);
                     material_model->evaluate(in, out);
+                    fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();
                   }
 
                 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);
+                    fe_face_values[introspection.variable("fluid velocity").extractor_vector()].get_function_values(current_linearization_point,
+                        face_current_fluid_velocity_values);
                   }
 
                 // ... check if the face is an outflow boundary by integrating the normal velocities

tests/new_outflow_test.prm

@@ -0,0 +1,155 @@
+# This test defines an initial blob of porosity implaced in a solid within a 2D cartesian box.
+# Because the porosity blob is less dense, the blob will rise to the surface of the model and
+# flow out of the top of the box. This tests the implementation of composition dependent outflow
+# boundary conditions.
+set Dimension                              = 2
+set Use years in output instead of seconds = true
+set End time                               = 2.5e3
+set Pressure normalization                 = surface
+set Maximum time step                      = 100
+set Surface pressure                       = 0
+set Use operator splitting                 = true
+
+# Define a function which 
+subsection Boundary velocity model
+  # set Prescribed velocity boundary indicators = top:function, bottom:function, right: function, left: function
+  set Zero velocity boundary indicators = bottom, left, top, right
+
+  subsection Function
+    set Variable names      = x,y,
+    set Function expression = 0;0
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10   # = Ra
+  end
+end
+
+subsection Formulation
+  set Formulation = Boussinesq approximation
+end
+
+subsection Solver parameters
+  set Temperature solver tolerance = 1e-10
+end
+
+# 10 km x 10 km box
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 10e3
+    set Y extent = 10e3
+
+    set X repetitions = 10
+    set Y repetitions = 10
+  end
+end
+
+# Uniform temperature of 293 K
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 293
+  end
+end
+
+subsection Boundary temperature model
+  set Allow fixed temperature on outflow boundaries = false
+  set List of model names = box  
+  set Fixed temperature boundary indicators   = top, bottom, left, right
+
+  subsection Box
+    set Bottom temperature = 293
+    set Top temperature    = 293
+    set Left temperature   = 293
+    set Right temperature  = 293
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 2
+  set Names of fields = porosity, bound_fluid
+  set Compositional field methods = darcy field, field # uncomment for darcy velocity case
+  # set Compositional field methods = field, field # uncomment for fluid velocity case
+end
+
+subsection Melt settings
+  set Include melt transport             = false # uncomment for darcy velocity case
+  # set Include melt transport             = true # uncomment for fluid velocity case
+end
+
+
+# Initialize a porosity of 1% (0.01) everywhere.
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,y,t
+    set Function constants  = radius=2.5e3, center_x=5e3, center_y=5e3, initial_porosity=0.01
+    set Function expression = if( ( (x - center_x)^2 + (y - center_y)^2 ) <= radius^2, initial_porosity, 0.0); 0.0
+  end
+end
+
+subsection Boundary composition model
+  set Fixed composition boundary indicators = right, left
+  set Allow fixed composition on outflow boundaries = false
+  set List of model names = initial composition
+end
+
+# The reactive fluid transport model allows us to set the parameters which
+# influence fluid velocity, such as the fluid viscosity, fluid density, and
+# the reference permeability. We set the fluid weakening factors (shear to
+# bulk viscosity ratio, exponential fluid weakening factor) to 0 for
+# simplicity. We use the zero solubility fluid-solid reaction scheme to prevent
+# the fluid from partitioning into the solid.
+subsection Material model
+  set Model name = reactive fluid transport
+
+  subsection Reactive Fluid Transport Model
+    set Base model                                       = visco plastic
+    set Reference fluid density                          = 1000
+    set Shear to bulk viscosity ratio                    = 0.
+    set Reference fluid viscosity                        = 10
+    set Reference permeability                           = 1e-6
+    set Exponential fluid weakening factor               = 0
+    set Fluid-solid reaction scheme                      = zero solubility
+  end
+
+  # Set the solid density to 3000 kg/m3, and set the minimum/maximum viscosity
+  # to 1e21 Pa s for an isoviscous model.
+  subsection Visco Plastic
+    set Reference temperature                   = 1600
+    set Prefactors for diffusion creep          = 5e-21
+    set Viscous flow law                        = diffusion
+    set Densities                               = 3000
+    set Viscosity averaging scheme              = harmonic
+    set Minimum viscosity                       = 1e21
+    set Maximum viscosity                       = 1e21
+    set Thermal expansivities                   = 0
+  end
+end
+
+# Set the global refinement to 0, 1 km x 1 km mesh.
+subsection Mesh refinement
+  set Initial global refinement          = 0
+  set Time steps between mesh refinement = 0
+end
+
+
+# Output the darcy velocity
+subsection Postprocess
+  set List of postprocessors = visualization
+
+  subsection Visualization
+    set List of output variables      = darcy velocity
+    set Time between graphical output = 250
+    set Output format                 = vtu
+  end
+end
+