fix bound_fluid.cc

Author: danieldouglas92

cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_darcy_with_fugacity.prm

@@ -51,7 +51,7 @@ subsection Solver parameters
     set GMRES solver restart length                     = 1000
     set Number of cheap Stokes solver steps             = 20000
     set Use full A block as preconditioner              = true
-    set Linear solver tolerance                         = 1e-4
+    set Linear solver tolerance                         = 1e-5
     set Maximum number of expensive Stokes solver steps = 0
   end
   subsection Operator splitting parameters
@@ -99,7 +99,7 @@ subsection Initial composition model
   set List of model names = world builder disequilibrium
   subsection World builder
     set List of relevant compositions = porosity, bound_fluid, peridotite, gabbro, MORB, sediment
-    set Disequilibrium percentage     = 10
+    set Disequilibrium percentage     = 5
   end
 end
 

cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_mckenzie_with_fugacity.prm

@@ -30,7 +30,7 @@ subsection Material model
   end
 end
 
-# Remove the darcy velocity from the list of output variables.
+# Remove the darcy velocity from the list of output variables, output the melt material properties.
 subsection Postprocess
   subsection Visualization
     set List of output variables      = material properties, strain rate, named additional outputs, melt material properties

cookbooks/2D_subduction_with_two_phase_flow/plugin/bound_fluid_source.cc

@@ -45,13 +45,17 @@ namespace aspect
     {
       ReactiveFluidTransport<dim>::evaluate(in, out);
       const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
+      const std::shared_ptr<ReactionRateOutputs<dim>> reaction_rate_out
+        = out.template get_additional_output_object<ReactionRateOutputs<dim>>();
 
-      // Do not allow the bound fluid content to decrease.
-      for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
-        if (out.reaction_terms[q][bound_fluid_idx] <= 0.0)
-          out.reaction_terms[q][bound_fluid_idx] = 0.0;
+      if (this->get_parameters().use_operator_splitting && reaction_rate_out != nullptr
+          && in.requests_property(MaterialProperties::reaction_rates))
+        {
+          for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
+            if (reaction_rate_out->reaction_rates[q][bound_fluid_idx] <= 0.0)
+              reaction_rate_out->reaction_rates[q][bound_fluid_idx] = 0.0;
+        }
     }
-
   }
 }
 
@@ -62,8 +66,8 @@ namespace aspect
   {
     ASPECT_REGISTER_MATERIAL_MODEL(BoundFluidSource,
                                    "reactive fluid transport bound fluid source",
-                                   "A simple material model that is like the "
-                                   "'reactive fluid transport' model, but prevents "
-                                   "the bound fluid content from decreasing.")
+                                   "A material model that is like the 'reactive "
+                                   "fluid transport' model, but modifies the reactions "
+                                   "to prevent the bound fluid content from decreasing.")
   }
 }

cookbooks/2D_subduction_with_two_phase_flow/plugin/world_builder_disequilibrium_fluid.cc

@@ -157,7 +157,7 @@ namespace aspect
                              "for all compositions that are not selected in the list. "
                              "By default the list is empty and the world builder is "
                              "evaluated for all compositional fields.");
-          prm.declare_entry ("Disequilibrium percentage", "10",
+          prm.declare_entry ("Disequilibrium percentage", "5",
                              Patterns::Double(0),
                              "A percentage value that defines the increase in the "
                              "bound_fluid compositional field from its equilibrium value. "
@@ -216,12 +216,10 @@ namespace aspect
     ASPECT_REGISTER_INITIAL_COMPOSITION_MODEL(WorldBuilderDisequilibrium,
                                               "world builder disequilibrium",
                                               "Specify the initial composition through the World Builder. "
-                                              "More information on the World Builder can be found at "
-                                              "\\url{https://geodynamicworldbuilder.github.io}. "
-                                              "Make sure to specify the location of the World Builder file "
-                                              "in the parameter 'World builder file'. It is possible to use "
-                                              "the World Builder only for selected compositional fields by "
-                                              "specifying the parameter 'List of relevant compositions'.")
+                                              "If the composition is not named bound_fluid, the composition "
+                                              "is set to the value returned by the World Builder. If the "
+                                              "composition is named bound_fluid, the value returned by the "
+                                              "World Builder is increased by a user defined percentage. ")
   }
 }
 #endif