Merge pull request #6214 from jdannberg/latent_heat_viscosity

fix latent heat viscosity prefactor

doc/modules/changes/20240124_jdannberg

@@ -0,0 +1,9 @@
+Fixed: The latent heat material model now correctly reads in
+the Viscosity prefactors that change the viscosities of
+individual phases for each compositional field. To make this
+work in a consistent way, the format of this input parameter
+is now the same as for other phase transition inputs (it is
+parsed as a map with keywords rather than a comma-separated
+list), which is an incompatible change in the input file.
+<br>
+(Juliane Dannberg, 2024/01/24)

source/material_model/latent_heat.cc

@@ -73,8 +73,8 @@ namespace aspect
             if ((composition_viscosity_prefactor != 1.0) && (composition.size() > 0))
               {
                 // geometric interpolation
-                out.viscosities[i] = (std::pow(10, ((1-composition[0]) * std::log10(eta*visc_temperature_dependence)
-                                                    + composition[0] * std::log10(eta*composition_viscosity_prefactor*visc_temperature_dependence))));
+                out.viscosities[i] = (std::pow(10., ((1-composition[0]) * std::log10(eta*visc_temperature_dependence)
+                                                     + composition[0] * std::log10(eta*composition_viscosity_prefactor*visc_temperature_dependence))));
               }
             else
               out.viscosities[i] = visc_composition_dependence * visc_temperature_dependence * eta;
@@ -105,10 +105,15 @@ namespace aspect
             // this means, that there are no actual density or viscosity "jumps", but
             // gradual transitions between the materials
             double phase_dependence = 0.0;
-            double viscosity_phase_dependence = 1.0;
+            std::vector<double> viscosity_phase_prefactors (this->n_compositional_fields()+1);
+
+            // Scale viscosity for shallowest phase for each chemical composition
+            viscosity_phase_prefactors[0] = phase_prefactors[0];
+            if (composition.size()>0)
+              viscosity_phase_prefactors[1] = phase_prefactors[phase_function.n_phases_for_each_composition()[0]];
 
             // Loop through phase transitions
-            for (unsigned int phase=0; phase<phase_function.n_phase_transitions(); ++phase)
+            for (unsigned int transition_index=0; transition_index<phase_function.n_phase_transitions(); ++transition_index)
               {
                 const double depth = this->get_geometry_model().depth(position);
                 const double pressure_depth_derivative = (depth > 0.0)
@@ -121,28 +126,37 @@ namespace aspect
                                                                            pressure,
                                                                            depth,
                                                                            pressure_depth_derivative,
-                                                                           phase);
+                                                                           transition_index);
 
                 const double phaseFunction = phase_function.compute_value(phase_in);
 
-                // Note that for the densities we have a list of jumps, so the index used
-                // in the loop corresponds to the index of the phase transition. For the
-                // viscosities we have a list of prefactors, which has one more entry
-                // for the first layer, so we have to use phase+1 as the index.
+                // For the densities we have a list of jumps, so the index used
+                // in the loop corresponds to the index of the phase transition.
                 if (composition.size()==0)
                   {
-                    phase_dependence += phaseFunction * density_jumps[phase];
-                    viscosity_phase_dependence *= 1. + phaseFunction * (phase_prefactors[phase+1]-1.);
+                    phase_dependence += phaseFunction * density_jumps[transition_index];
                   }
                 else if (composition.size()>0)
                   {
-                    if (transition_phases[phase] == 0)     // 1st compositional field
-                      phase_dependence += phaseFunction * density_jumps[phase] * (1.0 - composition[0]);
-                    else if (transition_phases[phase] == 1) // 2nd compositional field
-                      phase_dependence += phaseFunction * density_jumps[phase] * composition[0];
-
-                    viscosity_phase_dependence *= 1. + phaseFunction * (phase_prefactors[phase]-1.);
+                    if (transition_phases[transition_index] == 0)     // 1st compositional field
+                      phase_dependence += phaseFunction * density_jumps[transition_index] * (1.0 - composition[0]);
+                    else if (transition_phases[transition_index] == 1) // 2nd compositional field
+                      phase_dependence += phaseFunction * density_jumps[transition_index] * composition[0];
                   }
+
+                // For the viscosities we have a list of prefactors, which has one more entry
+                // for the first layer for each composition. Consequently, we have to use
+                // transition_index + 1 (+1 for additional compositions) as the index.
+                if (transition_index+1 < phase_function.n_phases_for_each_composition()[0])      // 1st compositional field
+                  viscosity_phase_prefactors[0] = std::pow(10.,
+                                                           phaseFunction * std::log10(viscosity_phase_prefactors[0] * phase_prefactors[transition_index+1])
+                                                           + (1. - phaseFunction) * std::log10(viscosity_phase_prefactors[0]));
+                else if (transition_index+2 < phase_function.n_phases_for_each_composition()[0]
+                         + phase_function.n_phases_for_each_composition()[1])                  // 2nd compositional field
+                  viscosity_phase_prefactors[1] = std::pow(10.,
+                                                           phaseFunction * composition[0] * std::log10(viscosity_phase_prefactors[1] * phase_prefactors[transition_index+2])
+                                                           + (1. - phaseFunction) * std::log10(viscosity_phase_prefactors[1]));
+
               }
 
             // fourth, pressure dependence of density
@@ -164,7 +178,12 @@ namespace aspect
             if (this->get_parameters().formulation == Parameters<dim>::Formulation::boussinesq_approximation)
               out.densities[i] = (reference_rho * density_temperature_dependence + density_composition_dependence + phase_dependence);
 
-            out.viscosities[i] = std::max(minimum_viscosity, std::min(maximum_viscosity, out.viscosities[i] * viscosity_phase_dependence));
+            if (composition.size()==0)
+              out.viscosities[i] = out.viscosities[i]*viscosity_phase_prefactors[0];
+            else
+              out.viscosities[i] = (std::pow(10., ((1-composition[0]) * std::log10(out.viscosities[i]*viscosity_phase_prefactors[0])
+                                                   + composition[0] * std::log10(out.viscosities[i]*viscosity_phase_prefactors[1]))));
+            out.viscosities[i] = std::max(minimum_viscosity, std::min(maximum_viscosity, out.viscosities[i]));
           }
 
           // Calculate entropy derivative
@@ -302,12 +321,15 @@ namespace aspect
                              "and -1 for none of them. "
                              "List must have the same number of entries as Phase transition depths. "
                              "Units: \\si{\\pascal\\per\\kelvin}.");
-          prm.declare_entry ("Viscosity prefactors", "",
-                             Patterns::List (Patterns::Double (0.)),
-                             "A list of prefactors for the viscosity for each phase. The reference "
-                             "viscosity will be multiplied by this factor to get the corresponding "
+          prm.declare_entry ("Viscosity prefactors", "all:1",
+                             Patterns::Anything(),
+                             "A list of prefactors for the viscosity for each phase. The ``Viscosity'' "
+                             "parameter (modified by the ``Composition viscosity prefactor'', depending on "
+                             "composition) will be multiplied by this factor to get the corresponding "
                              "viscosity for each phase. "
-                             "List must have one more entry than Phase transition depths. "
+                             "List must have the same number of entries as there are phases, that is one "
+                             "more than ``Phase transition depths'' for each composition that is used in the "
+                             "model. "
                              "Units: non-dimensional.");
           prm.declare_entry ("Minimum viscosity", "1e19",
                              Patterns::Double (0.),
@@ -355,26 +377,42 @@ namespace aspect
                           (Utilities::split_string_list(prm.get ("Phase transition density jumps")));
           transition_phases = Utilities::string_to_int
                               (Utilities::split_string_list(prm.get ("Corresponding phase for density jump")));
-          phase_prefactors = Utilities::string_to_double
-                             (Utilities::split_string_list(prm.get ("Viscosity prefactors")));
+
+          // Use the same syntax as in the PhaseFunction utilities to read in the viscosity prefactors
+          std::vector<std::string> compositional_field_names = this->introspection().get_composition_names();
+          // Establish that a background field is required here
+          compositional_field_names.insert(compositional_field_names.begin(),"background");
+          Utilities::MapParsing::Options options(compositional_field_names, "Viscosity prefactors");
+          options.allow_missing_keys = true;
+          options.allow_multiple_values_per_key = true;
+          options.check_values_per_key = true;
+          options.n_values_per_key = phase_function.n_phases_for_each_chemical_composition();
+
+          phase_prefactors = Utilities::MapParsing::parse_map_to_double_array (prm.get(options.property_name), options);
 
           const unsigned int n_transitions = phase_function.n_phase_transitions();
           if (density_jumps.size() != n_transitions ||
-              transition_phases.size() != n_transitions ||
-              phase_prefactors.size() != n_transitions+1)
+              transition_phases.size() != n_transitions)
             AssertThrow(false, ExcMessage("Error: At least one list that provides input parameters for phase "
                                           "transitions has the wrong size. The phase function object reports that "
                                           "there are " + std::to_string(n_transitions) + " transitions, "
                                           "therefore the material model expects " + std::to_string(n_transitions) +
-                                          " density jumps and corresponding phases, and "
-                                          + std::to_string(n_transitions+1) + " viscosity prefactors."));
+                                          " density jumps, corresponding phases, and viscosity prefactors."));
 
           // as the phase viscosity prefactors are all applied multiplicatively on top of each other,
           // we have to scale them here so that they are relative factors in comparison to the product
-          // of the prefactors of all phase above the current one
-          for (unsigned int phase=1; phase<phase_prefactors.size(); ++phase)
+          // of the prefactors of all phases above the current one
+          unsigned int index = 0;
+          for (unsigned int c=0; c<options.n_values_per_key.size(); ++c)
             {
-              phase_prefactors[phase] /= phase_prefactors[phase-1];
+              double product = 1.;
+              for (unsigned int phase=0; phase<options.n_values_per_key[c]; ++phase)
+                {
+                  if (phase>0)
+                    phase_prefactors[index] /= product;
+                  product *= phase_prefactors[index];
+                  index += 1;
+                }
             }
 
           if (thermal_viscosity_exponent!=0.0 && reference_T == 0.0)

tests/convection_box_phase_function.prm

@@ -116,7 +116,6 @@ subsection Material model
     set Phase transition Clapeyron slopes = -4050000
     set Phase transition density jumps = 200
     set Corresponding phase for density jump = 0
-    set Viscosity prefactors = 1,1
     set Minimum viscosity = 1e20
     set Maximum viscosity = 1e20
   end

tests/direct_solver_1.prm

@@ -98,7 +98,6 @@ subsection Material model
     # Viscosity is constant.
     set Thermal viscosity exponent                     = 0.0
     set Viscosity                                      = 8.44e21
-    set Viscosity prefactors                           = 1.0, 1.0
     set Composition viscosity prefactor                = 1.0
   end
 end

tests/direct_solver_2.prm

@@ -96,7 +96,6 @@ subsection Material model
     # Viscosity is constant.
     set Thermal viscosity exponent                     = 0.0
     set Viscosity                                      = 8.44e21
-    set Viscosity prefactors                           = 1.0, 1.0
     set Composition viscosity prefactor                = 1.0
   end
 end

tests/latent_heat.prm

@@ -97,7 +97,6 @@ subsection Material model
     # Viscosity is constant.
     set Thermal viscosity exponent                     = 0.0
     set Viscosity                                      = 8.44e21
-    set Viscosity prefactors                           = 1.0, 1.0
     set Composition viscosity prefactor                = 1.0
   end
 end

tests/latent_heat_chemical_viscosity.prm

@@ -0,0 +1,110 @@
+# This is a model setup to test if the viscosity stratification
+# in the latent heat material model works correctly with
+# compositional fields. The test has a different number and depth
+# of phase transitions for the two compositions.
+
+set Dimension = 2
+set Start time                             = 0
+set End time                               = 0
+set Use years in output instead of seconds = false
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 500000
+    set Y extent = 1000000
+    set Y repetitions = 2
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10.0
+  end
+end
+
+############### Boundary conditions
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators   = top, bottom
+  set List of model names = initial temperature
+end
+
+# We prescribe a constant downward flow.
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right, top, bottom
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 1600.0
+    set Variable names      = x,y
+  end
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,y
+    set Function expression = if(x>=250000, 1, 0);
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 1
+  set Names of fields = field1
+end
+
+subsection Material model
+  set Model name = latent heat
+
+  subsection Latent heat
+    set Phase transition density jumps                 = 100, 100, 100, 100, 100
+    set Corresponding phase for density jump           = 0, 0, 1, 1, 1
+    set Density differential for compositional field 1 = 0
+
+    set Phase transition depths                        = background:400000|600000, field1:300000|500000|700000
+    set Phase transition temperatures                  = background:1600|1600, field1:1600|1600|1600
+    set Phase transition Clapeyron slopes              = all:0
+
+    set Phase transition widths                        = all:50000
+    set Reference density                              = 3400
+    set Reference specific heat                        = 1000
+    set Reference temperature                          = 1600
+    set Thermal conductivity                           = 2.38
+
+    set Thermal expansion coefficient                  = 0.0
+    set Compressibility                                = 0.0
+
+    # Viscosity is constant.
+    set Thermal viscosity exponent                     = 0.0
+    set Viscosity                                      = 1e22
+    set Viscosity prefactors                           = background:1|10|1, field1:3.333|1|10|3.333
+    set Composition viscosity prefactor                = 1.0
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 4
+  set Time steps between mesh refinement = 0
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization
+
+    subsection Visualization
+    set List of output variables      = material properties
+    subsection Material properties
+      set List of material properties = viscosity, density
+    end
+
+    set Number of grouped files       = 1
+    set Output format                 = gnuplot
+  end
+end

tests/latent_heat_chemical_viscosity/screen-output

@@ -0,0 +1,26 @@
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+Number of active cells: 512 (on 5 levels)
+Number of degrees of freedom: 9,141 (4,290+561+2,145+2,145)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving temperature system... 0 iterations.
+   Solving field1 system ... 0 iterations.
+   Solving Stokes system (GMG)... 18+0 iterations.
+
+   Postprocessing:
+     Writing graphical output: output-latent_heat_chemical_viscosity/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
++----------------------------------------------+------------+------------+
++----------------------------------+-----------+------------+------------+
++----------------------------------+-----------+------------+------------+
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------