Enable different boundary conditions for different fields on different boundaries

include/aspect/boundary_composition/interface.h

@@ -212,6 +212,27 @@ namespace aspect
         bool
         allows_fixed_composition_on_outflow_boundaries() const;
 
+        /*
+         * Return whether on the given boundary the given field
+         * is fixed.
+         */
+        bool
+        get_component_mask_for_field(const types::boundary_id boundary_id,
+                                     const unsigned int compositional_field) const;
+
+        /*
+         * Return whether there are boundaries where only a subset
+         * of fields is fixed.
+         */
+        bool
+        boundaries_with_fixed_subset_of_fields_exist() const;
+
+        std::vector<unsigned int>
+        get_fixed_fields_on_boundary (const types::boundary_id boundary_id) const;
+
+        std::set<types::boundary_id>
+        get_fixed_boundaries_for_field (const unsigned int compositional_field) const;
+
         /**
          * For the current plugin subsystem, write a connection graph of all of the
          * plugins we know about, in the format that the
@@ -235,8 +256,28 @@ namespace aspect
                         << arg1
                         << "> among the names of registered boundary composition objects.");
       private:
+
         /**
-         * A list of enums of boundary composition operators that have been
+         * A list of boundary indicators that indicate for
+         * each plugin in the list of plugin_objects which boundary id
+         * it is responsible for. By default each plugin
+         * is active for all boundaries, but this list
+         * can be modified by derived classes to limit the application
+         * of plugins to specific boundaries.
+         */
+        std::vector <std::vector<types::boundary_id>> boundary_indicators;
+
+
+        /**
+         * A list of masks that specify for each plugin object
+         * for each boundary it applies to which compositional fields
+         *  are prescribed (true) and which are not (false).
+         */
+        std::vector <std::vector<ComponentMask>> masks_fields;
+
+        /**
+         * A list of enums of boundary composition operators for each
+         * boundary composition plugin that has been
          * requested in the parameter file. Each name is associated
          * with a model_name, and is used to modify the composition
          * boundary with the values from the current plugin.
@@ -254,6 +295,12 @@ namespace aspect
          * where material flows out of the domain.
          */
         bool allow_fixed_composition_on_outflow_boundaries;
+
+        /**
+         * Whether one or more boundaries only have fixed boundary conditions
+         * for a subset of fields.
+         */
+        bool boundaries_with_fixed_subset_of_fields = false;
     };
 
 

source/boundary_composition/box.cc

@@ -143,6 +143,8 @@ namespace aspect
       // Verify that each of the lists for boundary values
       // has the requisite number of elements if it is in the set
       // of prescribed boundary indicators.
+      // TODO Even though not all fields need to be fixed on a given boundary
+      // this plugin still requires that values are given for all fields.
       for (unsigned int f=0; f<2*dim; ++f)
         if (this->get_boundary_composition_manager().get_fixed_composition_boundary_indicators().count(f) != 0)
           AssertThrow (composition_values[f].size() == this->n_compositional_fields(),

source/particle/property/interface.cc

@@ -607,8 +607,8 @@ namespace aspect
                       for (unsigned int property_component = 0; property_component < n_property_components; ++property_component)
                         particle_properties.push_back(interpolated_properties[0][start_index+property_component]);
                     }
-                  // Otherwise use the value of the boundary condition of the corresponding compositional field
-                  else
+                  // Otherwise use the value of the boundary condition of the corresponding compositional field.
+                  else if (!manager.boundaries_with_fixed_subset_of_fields_exist())
                     {
                       const types::boundary_id boundary_id = cell->face(boundary_face)->boundary_id();
                       const unsigned int n_property_components = property_information.get_components_by_plugin_index(property_index);
@@ -634,6 +634,48 @@ namespace aspect
                           particle_properties.push_back(field_boundary_value);
                         }
                     }
+                  // Not all fields need to be fixed on a fixed boundary composition boundary. Those fields that
+                  // are not fixed, still need to be interpolated.
+                  else
+                    {
+                      const std::vector<std::vector<double>> interpolated_properties = interpolator.properties_at_points(particle_handler,
+                                                                                        std::vector<Point<dim>> (1,particle_location),
+                                                                                        ComponentMask(property_information.n_components(),true),
+                                                                                        found_cell);
+
+                      const unsigned int start_index = property_information.get_position_by_plugin_index(property_index);
+                      const types::boundary_id boundary_id = cell->face(boundary_face)->boundary_id();
+                      const unsigned int n_property_components = property_information.get_components_by_plugin_index(property_index);
+
+                      for (unsigned int particle_property_component = 0; particle_property_component < n_property_components; ++particle_property_component)
+                        {
+                          // Ask the particle property which field index to use to evaluate boundary condition
+                          const AdvectionField field_to_use = p->advection_field_for_boundary_initialization(particle_property_component);
+
+                          Assert(field_to_use.is_temperature() == false,
+                                 ExcMessage("Interpolating temperature boundary conditions to particles is not supported."));
+                          Assert(field_to_use.compositional_variable < this->n_compositional_fields(),
+                                 ExcMessage("The composition index specified in the function"
+                                            "advection_field_for_boundary_initialization() by the particle property "
+                                            "<" + property_information.get_field_name_by_index(property_index) +"> is "
+                                            "larger than the number of compositional fields. This is not possible "
+                                            "and likely a bug in the particle property plugin."));
+
+                          // Check whether this field is fixed on this boundary, ...
+                          if (manager.get_component_mask_for_field(boundary_id, field_to_use.compositional_variable))
+                            {
+                              const double field_boundary_value = manager.boundary_composition(boundary_id,
+                                                                                               particle_location,
+                                                                                               field_to_use.compositional_variable);
+
+                              particle_properties.push_back(field_boundary_value);
+                            }
+                          // otherwise interpolate.
+                          else
+                            particle_properties.push_back(interpolated_properties[0][start_index+particle_property_component]);
+
+                        }
+                    }
 
                   break;
                 }

source/simulator/assemblers/advection.cc

@@ -848,7 +848,8 @@ namespace aspect
                face->boundary_id()
              )
              != this->get_fixed_composition_boundary_indicators().end())
-           && (!advection_field.is_temperature())))
+           && (!advection_field.is_temperature())
+           && this->get_boundary_composition_manager().get_component_mask_for_field(face->boundary_id(), advection_field.compositional_variable)))
         {
           /*
            * We are on a face of a Dirichlet temperature or composition boundary.

source/simulator/core.cc

@@ -764,7 +764,7 @@ namespace aspect
     {
       // obtain the boundary indicators that belong to Dirichlet-type
       // composition boundary conditions and interpolate the composition
-      // there
+      // there, but only for those fields that are fixed.
       for (unsigned int c=0; c<introspection.n_compositional_fields; ++c)
         {
           // If we do not want to prescribe Dirichlet boundary conditions on outflow boundaries,
@@ -774,7 +774,7 @@ namespace aspect
             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())
+            for (const auto p : boundary_composition_manager.get_fixed_boundaries_for_field(c))
               {
                 VectorFunctionFromScalarFunctionObject<dim> vector_function_object(
                   [&] (const Point<dim> &x) -> double

source/simulator/entropy_viscosity.cc

@@ -380,7 +380,7 @@ namespace aspect
               ?
               boundary_temperature_manager.get_fixed_temperature_boundary_indicators()
               :
-              boundary_composition_manager.get_fixed_composition_boundary_indicators();
+              boundary_composition_manager.get_fixed_boundaries_for_field(advection_field.compositional_variable);
             const std::set<types::boundary_id> &tangential_velocity_boundaries =
               boundary_velocity_manager.get_tangential_boundary_velocity_indicators();
             const std::set<types::boundary_id> &zero_velocity_boundaries =

source/volume_of_fluid/assembler.cc

@@ -196,7 +196,8 @@ namespace aspect
           if ( this->get_fixed_composition_boundary_indicators().find(
                  face->boundary_id()
                )
-               != this->get_fixed_composition_boundary_indicators().end())
+               != this->get_fixed_composition_boundary_indicators().end()
+               && this->get_boundary_composition_manager().get_component_mask_for_field(face->boundary_id(), field.composition_index))
             {
               // We are in the case of a Dirichlet BC, so calculate the fluid flux if the boundary is upwind.
               for (unsigned int q=0; q<n_f_q_points; ++q)

tests/boundary_composition_plugin_per_boundary.prm

@@ -0,0 +1,113 @@
+# A test for specifying fixed compositions for only
+# a subset of the compositional fields. The subset
+# also use different boundary composition plugins to
+# set the fixed values.
+
+set Dimension                              = 2
+set Use years instead of seconds           = true
+set End time                               = 0
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 660000
+    set Y extent = 660000
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 3
+  set Names of fields = field1, field2, field3
+end
+
+subsection Boundary composition model
+  set List of model names = box, function, initial composition
+  set List of model operators = maximum, subtract, add
+  set Fixed composition boundary indicators   = top;field1:box,\
+                                                top;field2:box|function,\
+                                                left,\
+                                                bottom;field1|field2|field3:box|function|initial composition,\
+                                                right
+
+  subsection Box
+    set Left composition = 15, 15, 15
+    set Right composition = 0, 0, 0
+    set Top composition = 0, 0, 0
+    set Bottom composition = 7, 7, 7
+  end
+
+  subsection Function
+    set Function expression = if (x>500000,3,0); if (x>500000,3,0); if (x>500000,3,0);
+  end
+end
+
+subsection Initial composition model
+  set List of model names = function
+
+  subsection Function
+    set Function expression = if ((x-330000)*(x-330000)+(y-330000)*(y-330000) < 100000*100000,-1,1.5);\
+                              if ((x-330000)*(x-330000)+(y-330000)*(y-330000) < 100000*100000,-1,1.5);\
+                              if ((x-330000)*(x-330000)+(y-330000)*(y-330000) < 100000*100000,-1,1.5);
+  end
+end
+
+subsection Initial temperature model
+  set List of model names = function
+
+  subsection Function
+    set Function expression = if ((x-330000)*(x-330000)+(y-330000)*(y-330000) < 100000*100000,-1,1.5)
+  end
+end
+
+subsection Boundary temperature model
+  set List of model names = box, constant, function, initial temperature
+  set List of model operators = maximum, add, subtract, add
+  set Fixed temperature boundary indicators   = top,left,bottom,right
+
+  subsection Box
+    set Left temperature = 15
+  end
+
+  subsection Constant
+    set Boundary indicator to temperature mappings = bottom:7, left:0, top:0, right:0
+  end
+
+  subsection Function
+    set Function expression = if (x>500000,3,0)
+  end
+end
+
+subsection Boundary velocity model
+  set Prescribed velocity boundary indicators = bottom:function,left:function,right:function,top:function
+
+  subsection Function
+    set Function expression = 1;0
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10
+  end
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Viscosity = 1e21
+  end
+end
+
+subsection Mesh refinement
+  set Initial global refinement                = 4
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+end
+
+subsection Postprocess
+  set List of postprocessors = composition statistics
+end

tests/boundary_composition_plugin_per_boundary/screen-output

@@ -0,0 +1,16 @@
+
+Number of active cells: 256 (on 5 levels)
+Number of degrees of freedom: 4,645 (2,178+289+1,089+1,089)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving C_1 system ... 0 iterations.
+   Solving Stokes system (GMG)... 14+0 iterations.
+
+   Postprocessing:
+     Compositions min/max/mass: -1.5/16.5/6.759e+11
+
+Termination requested by criterion: end time
+
+
+

tests/boundary_composition_plugin_per_boundary/statistics

@@ -0,0 +1,16 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for Stokes solver
+# 11: Velocity iterations in Stokes preconditioner
+# 12: Schur complement iterations in Stokes preconditioner
+# 13: Minimal value for composition C_1
+# 14: Maximal value for composition C_1
+# 15: Global mass for composition C_1
+0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 1089 0 0 14 15 15 -1.50000000e+00 1.65000000e+01 6.75945703e+11 

tests/boundary_composition_plugin_per_boundary_fixed_on_outflow.prm

@@ -0,0 +1,12 @@
+# A test for specifying fixed compositions for only
+# a subset of the compositional fields when the composition
+# is also fixed on outflow boundaries.
+include $ASPECT_SOURCE_DIR/tests/boundary_composition_plugin_per_boundary.prm
+
+set Dimension                              = 2
+
+subsection Boundary composition model
+  # The right boundary experiences outflow
+  set Allow fixed composition on outflow boundaries = true
+end
+

tests/boundary_composition_plugin_per_boundary_fixed_on_outflow/screen-output

@@ -0,0 +1,16 @@
+
+Number of active cells: 256 (on 5 levels)
+Number of degrees of freedom: 4,645 (2,178+289+1,089+1,089)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving C_1 system ... 0 iterations.
+   Solving Stokes system (GMG)... 14+0 iterations.
+
+   Postprocessing:
+     Compositions min/max/mass: -1.5/16.5/6.759e+11
+
+Termination requested by criterion: end time
+
+
+

tests/boundary_composition_plugin_per_boundary_fixed_on_outflow/statistics

@@ -0,0 +1,16 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Iterations for temperature solver
+# 9: Iterations for composition solver 1
+# 10: Iterations for Stokes solver
+# 11: Velocity iterations in Stokes preconditioner
+# 12: Schur complement iterations in Stokes preconditioner
+# 13: Minimal value for composition C_1
+# 14: Maximal value for composition C_1
+# 15: Global mass for composition C_1
+0 0.000000000000e+00 0.000000000000e+00 256 2467 1089 1089 0 0 14 15 15 -1.50000000e+00 1.65000000e+01 6.75945703e+11