Adding velocity option to isosurfaces refinement strategy

doc/modules/changes/20260730_hyunseong96_2

@@ -0,0 +1,3 @@
+Added: Isosurfaces mesh refinement strategy now has an option for velocity.
+<br>
+(Hyunseong Kim 2026/07/29)

doc/sphinx/parameters/Mesh_20refinement.md

@@ -173,16 +173,18 @@ On the other hand, for discontinuous finite elements (see the ‘Use discont
 
 ‘density’: A mesh refinement criterion that computes refinement indicators from a field that describes the spatial variability of the density, $\rho$. Because this quantity may not be a continuous function ($\rho$ and $C_p$ may be discontinuous functions along discontinuities in the medium, for example due to phase changes), we approximate the gradient of this quantity to refine the mesh. The error indicator defined here takes the magnitude of the approximate gradient and scales it by $h_K^{1+d/2}$ where $h_K$ is the diameter of each cell and $d$ is the dimension. This scaling ensures that the error indicators converge to zero as $h_K\rightarrow 0$ even if the energy density is discontinuous, since the gradient of a discontinuous function grows like $1/h_K$.
 
-‘isosurfaces’: A mesh refinement criterion that computes coarsening and refinement indicators between two isosurfaces of specific field entries (e.g. temperature, composition).
+‘isosurfaces’: A mesh refinement criterion that computes coarsening and refinement indicators between two isosurfaces of specific field entries (e.g. temperature, norm of velocity, composition).
 
 The way these indicators are derived between pairs of isosurfaces is by checking whether the solutions of specific fields are within the ranges of the isosurface values given. If these conditions hold, then coarsening and refinement indicators are set such that the mesh refinement levels lies within the range of levels given. Usage of this plugin allows the user to put a conditional minimum and maximum refinement function onto fields that they are interested in.
 
-For now, only temperature and compositional fields are allowed as field entries. The key words could be ’Temperature’ or one of the names of the compositional fields which are either specified by user or set up as C\_0, C\_1, etc.
+For now, only temperature, norm of velocity and compositional fields are allowed as field entries. The key words could be ’Temperature’, ’Norm of velocity’ or one of the names of the compositional fields which are either specified by user or set up as C\_0, C\_1, etc.
 
 Usage: A list of isosurfaces separated by semi-colons (;). Each isosurface entry consists of multiple entries separated by a comma. The first two entries indicate the minimum and maximum refinement levels respectively. The entries after the first two describe the fields the isosurface applies to, followed by a colon (:), which again is followed by the minimum and maximum field values separated by a bar (|). An example for two isosurface entries is ’0, 2, Temperature: 300 | 600; 2, 2, C\_1: 0.5 | 1’. If both isoterm entries are triggered at the same location and the current refinement level is 1, it means that the first isoline will not set any flag and the second isoline will set a refinement flag. This means the cell will be refined. If both the coarsening and refinement flags are set, preference is given to refinement.
 
 The minimum and maximum refinement levels per isosurface can be provided in absolute values relative to the global minimum and maximum refinement. This is done with the ’min’ and ’max’ key words. For example: ’set Isosurfaces = max-2,  max,    Temperature: 0 | 600 ; min + 1,min+2, Temperature: 1600 | 3000,   C\_2 : 0.0 | 0.5’.
 
+Norm of velocity key word follows the unit selected from parameter ’Use years instead of seconds’. If ”Use years instead of seconds’ is true, velocity unit is \si{\meter}/yr. If false, the unit is \si{\meter\per\second}.
+
 ‘maximum refinement function’: A mesh refinement criterion that ensures a maximum refinement level described by an explicit formula with the depth or position as argument. Which coordinate representation is used is determined by an input parameter. Whatever the coordinate system chosen, the function you provide in the input file will by default depend on variables ‘x’, ‘y’ and ‘z’ (if in 3d). However, the meaning of these symbols depends on the coordinate system. In the Cartesian coordinate system, they simply refer to their natural meaning. If you have selected ‘depth’ for the coordinate system, then ‘x’ refers to the depth variable and ‘y’ and ‘z’ will simply always be zero. If you have selected a spherical coordinate system, then ‘x’ will refer to the radial distance of the point to the origin, ‘y’ to the azimuth angle and ‘z’ to the polar angle measured positive from the north pole. Note that the order of spherical coordinates is r,phi,theta and not r,theta,phi, since this allows for dimension independent expressions. Each coordinate system also includes a final ‘t’ variable which represents the model time, evaluated in years if the ’Use years instead of seconds’ parameter is set, otherwise evaluated in seconds. After evaluating the function, its values are rounded to the nearest integer.
 
 The format of these functions follows the syntax understood by the muparser library, see {ref}`sec:run-aspect:parameters-overview:muparser-format`.
@@ -342,6 +344,8 @@ If the list of scaling factors given in this parameter is empty, then this indic
 **Documentation:** A list of isosurfaces separated by semi-colons (;). Each isosurface entry consists of multiple entries separated by a comma. The first two entries indicate the minimum and maximum refinement levels respectively. The entries after the first two describe the fields the isosurface applies to, followed by a colon (:), which is again followed by the minimum and maximum property values separated by bar (|). An example for an isosurface is ’0, 2, Temperature: 300 | 600; 2, 2, C\_1: 0.5 | 1’. In this example the mesh refinement is kept between level 0 and level 2 if the temperature is between 300 and 600 and at level 2 when the compositional field C\_1 is between 0.5 and 1. If both happen at the same location and the current refinement level is 1, it means that the first isoline will not set any flag and the second isoline will set a refinement flag. This means the cell will be refined. If both the coarsening and refinement flags are set, preference is given to refinement.
 
 The first two entries for each isosurface, describing the minimum and maximum grid levels, can be two numbers or contain one of the key values ’min’ and ’max’. This indicates the key will be replaced with the global minimum and maximum refinement levels. The ’min’ and ’max’ keys also accept adding values to be added or subtracted from them respectively. This is done by adding a ’+’ or ’-’ and a number behind them (e.g. min+2 or max-1). Note that you can’t subtract a value from a minimum value or add a value to the maximum value. If, for example, ‘max-4‘ drops below the minimum or ‘min+4‘ goes above the maximum, it will simply use the global minimum and maximum values respectively. The same holds for any mesh refinement level below the global minimum or above the global maximum.
+
+Norm of velocity key word follows the unit selected from parameter ’Use years instead of seconds’. If ”Use years instead of seconds’ is true, velocity unit is \si{\meter}/yr. If false, the unit is \si{\meter\per\second}.
 ::::
 
 (parameters:Mesh_20refinement/Maximum_20refinement_20function)=

include/aspect/mesh_refinement/isosurfaces.h

@@ -37,6 +37,7 @@ namespace aspect
       enum class PropertyType
       {
         Temperature,
+        Norm_of_velocity,
         Composition
       };
 
@@ -47,8 +48,8 @@ namespace aspect
            * Constructor. Converts a property name into a structure containing a property type
            * and an index. If the property contains multiple items (e.g. the property compositional field has a
            * field index) the index referring to the particular item of that property is stored in the variable index.
-           * @param property_name The name of a property, which can be Temperature for the temperature field or
-           * the name of a compositional field listed in the parameter available_compositions.
+           * @param property_name The name of a property, which can be Temperature for the temperature field, Norm of velocity
+           * for the velocity field or the name of a compositional field listed in the parameter available_compositions.
            * @param available_compositions A list of names of the available compositional fields.
            */
           Property(const std::string &property_name,
@@ -94,7 +95,7 @@ namespace aspect
      * plugin allows for setting a minimum and a maximum refinement level in
      * a part of the model domain where a variable/property (e.g. Temperature)
      * is between two values (e.g. two isotherms of 274K and 1600K). This is
-     * currently implemented for temperature and compositions.
+     * currently implemented for temperature, norm of velocity and compositions.
      *
      * @ingroup MeshRefinement
      */

source/mesh_refinement/isosurfaces.cc

@@ -62,6 +62,12 @@ namespace aspect
             index = 0;
             found = true;
           }
+        else if (property_name == "Norm of velocity")
+          {
+            type = PropertyType::Norm_of_velocity;
+            index = 0;
+            found = true;
+          }
         else
           {
             auto p = std::find(available_compositions.begin(), available_compositions.end(), property_name);
@@ -76,7 +82,7 @@ namespace aspect
           {
             // The property name has not been found. This could be because the compositional field is not present.
             // Abort and warn the user.
-            std::string key_list = "Temperature";
+            std::string key_list = "Temperature, Norm of velocity";
             for (auto &composition : available_compositions)
               key_list += ", " + composition;
 
@@ -193,6 +199,11 @@ namespace aspect
                             {
                               values[index] = in.temperature[i_quad];
                             }
+                          else if (isosurface.properties[index].type == internal::PropertyType::Norm_of_velocity)
+                            {
+                              const double time_scale = this->get_parameters().convert_to_years ? constants::year_in_seconds : 1.0;
+                              values[index] = in.velocity[i_quad].norm() * time_scale;
+                            }
                           else if (isosurface.properties[index].type == internal::PropertyType::Composition)
                             {
                               values[index] = in.composition[i_quad][isosurface.properties[index].index];
@@ -283,7 +294,11 @@ namespace aspect
                              "number behind them (e.g. min+2 or max-1). Note that you can't subtract a value from a minimum value or "
                              "add a value to the maximum value. If, for example, `max-4` drops below the minimum or `min+4` goes above the "
                              "maximum, it will simply use the global minimum and maximum values respectively. The same holds for any "
-                             "mesh refinement level below the global minimum or above the global maximum.");
+                             "mesh refinement level below the global minimum or above the global maximum."
+                             "\n\n"
+                             "Norm of velocity key word follows the unit selected from parameter 'Use years instead of seconds'. "
+                             "If ''Use years instead of seconds' is true, velocity unit is \\si{\\meter}/yr. "
+                             "If false, the unit is \\si{\\meter\\per\\second}.") ;
 
         }
         prm.leave_subsection();
@@ -366,7 +381,7 @@ namespace aspect
                                               "isosurfaces",
                                               "A mesh refinement criterion that computes "
                                               "coarsening and refinement indicators between two isosurfaces of "
-                                              "specific field entries (e.g. temperature, composition)."
+                                              "specific field entries (e.g. temperature, norm of velocity, composition)."
                                               "\n\n"
                                               "The way these indicators are derived between pairs of isosurfaces is by "
                                               "checking whether the solutions of specific "
@@ -377,8 +392,8 @@ namespace aspect
                                               "minimum and maximum refinement function onto fields that they "
                                               "are interested in."
                                               "\n\n"
-                                              "For now, only temperature and compositional fields are allowed as "
-                                              "field entries. The key words could be 'Temperature' or one of the names "
+                                              "For now, only temperature, norm of velocity and compositional fields are allowed as "
+                                              "field entries. The key words could be 'Temperature', 'Norm of velocity' or one of the names "
                                               "of the compositional fields which are either specified by user or set up "
                                               "as C\\_0, C\\_1, etc."
                                               "\n\n"
@@ -397,6 +412,10 @@ namespace aspect
                                               "The minimum and maximum refinement levels per isosurface can be provided in absolute values relative to "
                                               "the global minimum and maximum refinement. This is done with the 'min' and 'max' key words. For example: "
                                               "'set Isosurfaces = max-2,  max,    Temperature: 0 | 600 ; min + 1,min+2, Temperature: 1600 | 3000,   C\\_2 : 0.0 | 0.5'."
+                                              "\n\n"
+                                              "Norm of velocity key word follows the unit selected from parameter 'Use years instead of seconds'. "
+                                              "If ''Use years instead of seconds' is true, velocity unit is \\si{\\meter}/yr. "
+                                              "If false, the unit is \\si{\\meter\\per\\second}."
                                              )
   }
 }

tests/isosurfaces_simple_box_2D_custom_fieldnames_fail/screen-output

@@ -13,7 +13,7 @@ An error occurred in file <isosurfaces.cc> in function
 Additional information: 
     The key given for the isosurface could not be converted. The provided
     key was: C_2. The following keys are allowed for this model:
-    Temperature, field1, FieldB, fIEld_C.
+    Temperature, Norm of velocity, field1, FieldB, fIEld_C.
 
 Stacktrace:
 (rest of the output replaced by default.sh script)

tests/isosurfaces_simple_box_2D_v.prm

@@ -0,0 +1,91 @@
+# This input file tests the isosurface mesh refinement plugin for norm of velocity.
+# Vertical velocity (v) varies by left and right side of the model. left half has 8 m/yr, and right side has 2 m/yr.
+# By the isosurface criterion, refinement levels are 3 for norm of velocity between 5 and 10 m/yr and 1 for norm of velocity between 0 and 5 m/yr.
+
+set Dimension = 2
+set CFL number                             = 1.0
+set End time                               = 0
+set Start time                             = 0
+set Surface pressure                       = 0
+set Use years instead of seconds           = true
+set Nonlinear solver scheme                = single Advection, no Stokes
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 10
+  end
+end
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 10000
+    set Y extent = 10000
+    set X repetitions = 2
+    set Y repetitions = 2
+  end
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Variable names = x, y
+    set Function constants = d= 10000
+    set Function expression = 100
+  end
+end
+
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = top, bottom
+  set List of model names = box
+
+  subsection Box
+    set Bottom temperature    = 100
+    set Top temperature    = 100
+  end
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density = 3340
+    set Reference specific heat = 1200
+    set Thermal expansion coefficient = 3e-5
+    set Viscosity = 1e21
+  end
+end
+
+subsection Mesh refinement
+  set Initial global refinement          = 2
+  set Initial adaptive refinement        = 1
+  set Time steps between mesh refinement = 0
+  set Strategy = isosurfaces
+
+  subsection Isosurfaces
+    set Isosurfaces = 3,    3,    Norm of velocity: 5 | 10; \
+                      1,    1,    Norm of velocity: 0 | 5
+  end
+end
+
+subsection Prescribed Stokes solution
+  set Model name =  function
+
+  subsection Velocity function
+    set Variable names = x, y, t
+    set Function expression = 0; if (x < 5000, 8, 2)
+  end
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization
+
+  subsection Visualization
+    set Output format            = gnuplot
+  end
+end

tests/isosurfaces_simple_box_2D_v/screen-output

@@ -0,0 +1,20 @@
+
+Number of active cells: 64 (on 3 levels)
+Number of degrees of freedom: 948 (578+81+289)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+
+Number of active cells: 148 (on 4 levels)
+Number of degrees of freedom: 2,153 (1,318+176+659)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+
+   Postprocessing:
+     Writing graphical output: output-isosurfaces_simple_box_2D_v/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
+