Allow recursive includes in parameter files

benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_t.prm

@@ -1,131 +1,17 @@
-# Like the poiseuille_2d.prm test and based on the nonlinear channel
+# Like input_v.prm and based on the nonlinear channel
 # flow benchmark. This is used to test the traction boundary conditions
 # of the Newton Stokes solver.
 
-set Output directory = output_t
-set Dimension = 2
-set CFL number                             = 1.0
-set Maximum time step                      = 1
-set End time                               = 0
-set Start time                             = 0
-set Adiabatic surface temperature          = 0
-set Surface pressure                       = 0
-set Use years in output instead of seconds = false
-set Nonlinear solver scheme = single Advection, iterated Newton Stokes
-set Max nonlinear iterations = 150
-set Nonlinear solver tolerance = 1e-14
-set Additional shared libraries = libsimple_nonlinear.so
-
-subsection Solver parameters
-  subsection Stokes solver parameters
-    set Linear solver tolerance = 1e-7
-    set Stokes solver type = block AMG
-  end
-
-  subsection Newton solver parameters
-    set Max pre-Newton nonlinear iterations = 3
-    set Nonlinear Newton solver switch tolerance = 1e-20
-    set Max Newton line search iterations = 5
-    set Maximum linear Stokes solver tolerance = 1e-2
-    set Use Newton residual scaling method = false
-    set Use Newton failsafe = false
-    set Stabilization preconditioner = SPD
-    set Stabilization velocity block = SPD
-  end
-end
-
-subsection Boundary temperature model
-  set List of model names = box
-  set Fixed temperature boundary indicators   = bottom, top
-
-  subsection Box
-    set Left temperature = 0
-  end
-end
-
-subsection Initial temperature model
-  set Model name = function
-
-  subsection Function
-    set Function expression = 0
-  end
-end
-
-subsection Gravity model
-  set Model name = vertical
-
-  subsection Vertical
-    set Magnitude = 0
-  end
-end
-
-subsection Geometry model
-  set Model name = box
+include $ASPECT_SOURCE_DIR/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm
 
-  subsection Box
-    set X extent = 10e3
-    set Y extent = 8e3
-    set Y repetitions = 1
-  end
-end
-
-subsection Material model
-  set Model name = simple nonlinear
-
-  subsection Simple nonlinear
-    set Minimum viscosity = 1e19
-    set Maximum viscosity = 1e24
-    set Stress exponent = 3
-    set Viscosity averaging p = -1
-    set Viscosity prefactor = 1e-37
-    set Use deviator of strain-rate = false
-  end
-end
-
-subsection Mesh refinement
-  set Initial adaptive refinement        = 0
-  set Initial global refinement          = 4
-end
+set Output directory = output_t
 
 subsection Boundary velocity model
-  set Zero velocity boundary indicators       = left, right
   # pressure bc: Prescribe a zero vertical velocity component on the vertical boundaries
   set Prescribed velocity boundary indicators = bottom x: function, top x: function
-
-  subsection Function
-    set Function constants = n = 3
-    set Variable names = x,y
-
-    # For velocity boundary conditions both are used, for pressure boundary conditions only the first (x) component
-    set Function expression = 0;(1e-37/(n+1))*((1e9/8e3)^n)*(((5e3)^(n+1))-(abs(x-(5e3))^(n+1)));
-  end
 end
 
 subsection Boundary traction model
   # pressure bc: Prescribe a horizontal traction on the vertical boundaries
   set Prescribed traction boundary indicators = bottom y: function, top y: function
-  subsection Function
-    set Variable names = x,y
-
-    # We want to prescribe a pressure of 1e9 at the bottom boundary
-    # and 0 at the top boundary.
-    # The traction in this case is defined as:
-    # tau =  - pressure * normal_vector.
-    # On the bottom boundary, the outward pointing normal vector is
-    # (0;-1). On the top (0;1).
-    # Therefore:
-    # Top boundary:    tau = - pressure(top) (0;1) = - (0) (0;1) = (0;0).
-    # Bottom boundary: tau = - pressure(bottom) (0;-1) = - (1e9) (0;-1) = (0;1e9).
-    # Conveniently, the traction is the same on both boundaries.
-    set Function expression = 0; y>4e3 ? 0 : 1e9
-  end
-end
-
-subsection Postprocess
-  set List of postprocessors = velocity statistics, pressure statistics, mass flux statistics, visualization
-
-  subsection Visualization
-    set List of output variables = material properties, strain rate
-    set Time between graphical output = 2
-  end
 end

benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_t_gmg.prm

@@ -3,7 +3,6 @@
 
 include $ASPECT_SOURCE_DIR/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_t.prm
 
-set Additional shared libraries = libsimple_nonlinear.so
 set Output directory = results/output-t-gmg
 
 subsection Material model

benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v_gmg.prm

@@ -3,7 +3,6 @@
 
 include $ASPECT_SOURCE_DIR/benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/input_v.prm
 
-set Additional shared libraries = libsimple_nonlinear.so
 set Output directory = results/output-v-gmg
 
 subsection Material model

doc/modules/changes/20241204_gassmoeller

@@ -0,0 +1,8 @@
+Added: ASPECT can now use the 'include' keyword in parameter
+files to include other parameter files, even if those other
+parameter files include parameter files themselves. ASPECT now
+also respects the parameters 'Dimension' and 'Additional shared
+libraries' in the included parameter files, which was not the
+case before.
+<br>
+(Rene Gassmoeller, 2024/12/04)

doc/sphinx/user/run-aspect/parameters-overview/structure.md

@@ -39,7 +39,7 @@ end in backslashes. The underlying implementation always eats whitespace at the
 each continuing line, but not before the backslash. This means that the parameter file
 ```
 set Some parameter = abc\
-def
+  def
 ```
 is equivalent to
 ```
@@ -48,3 +48,35 @@ set Some parameter = abcdef
 that is, with no space between `abc` and `def` despite the leading whitespace at the beginning of
 the second line. If you do want space between these two parts, you need to add it before the
 backslash in the first of the two lines.
+:::
+
+:::{note}
+If you want to run several models that are all small modifications of a base model you can
+`include` the base model in the modified model parameter files to include its parameters.
+This means that in a parameter file `file_a.prm` that contains
+```
+set Some parameter = abc
+include file_b.prm
+set Some other parameter = def
+```
+the content of file `file_b.prm` will be inserted at the position of the
+include statement. If `file_b.prm` includes further include statements, these will
+be recursively substituted until no statements remain or you reach the maximum number
+of supported include statements, in which case ASPECT will stop and output an error message.
+
+Note, that if the same parameter is set twice in parameter files, the last
+occurrence will overwrite the earlier occurrence(s). Thus, in the example above if `file_b.prm` contains both
+`Some parameter` and `Some other parameter`, then the final file will use the value of
+`Some parameter` from `file_b.prm`, but the value of `Some other parameter` from `file_a.prm`.
+
+Also note, that the include statement can include the file path as relative or absolute path,
+and you can also reference the original ASPECT source directory using the string `$ASPECT_SOURCE_DIR`.
+Thus the three include statements
+```
+include file_a.prm
+include /home/user/aspect/file_a.prm
+include $ASPECT_SOURCE_DIR/file_a.prm
+```
+could all include the same file, but the second and third statement are independent from
+your current working directory, while the first one depends on where you execute ASPECT.
+:::

source/main.cc

@@ -382,6 +382,43 @@ read_parameter_file(const std::string &parameter_file_name,
         }
 
       input_as_string = aspect::Utilities::read_and_distribute_file_content(parameter_file_name, comm);
+
+      // Search and replace include directives in the input file.
+      std::match_results<std::string::const_iterator> matches;
+      const std::string search_regex = "(?:^|\n)[ \t]*include[ \t]+(.*?)[ \t]*(?:#|\n|$)";
+      const std::string replace_regex = "(?:^|\n)[ \t]*include[ \t]+.*";
+
+      unsigned int n_included_files = 0;
+
+      while (std::regex_search(input_as_string, matches, std::regex(search_regex)))
+        {
+          // Make sure we are not circularly including files. This is not easily possible
+          // by making sure included files are unique, because we may have multiple
+          // files including the same file in a non-circular way. So we just limit
+          // the number of included files to a reasonable number.
+          Assert(n_included_files < 15,
+                 dealii::ExcMessage("Too many included files in parameter file. You likely have a circular include."));
+
+          // Since the line as a whole matched, the 'matches' variable needs to
+          // contain two entries: [0] denotes the whole line, and [1] the
+          // part that was matched by the '(.*?)' expression.
+          Assert (matches.size() == 2, dealii::ExcInternalError());
+
+          // Expand ASPECT_SOURCE_DIR, but only in the included file name (to keep the formatting
+          // of all parameter files intact, which we will copy into the output directory).
+          const std::string included_filename = aspect::Utilities::expand_ASPECT_SOURCE_DIR(matches.str(1));
+
+          // Prepend a newline character to the included file content, because if the include directive
+          // is not in the first line, we will replace one newlince character
+          // from the original input file in the regex_replace below.
+          const std::string included_file_content = '\n' + aspect::Utilities::read_and_distribute_file_content(included_filename, comm);
+
+          // Replace the include line with the content of the included file. Note that we only replace the first
+          // include line we find (there may be several, which we will replace in subsequent iterations).
+          input_as_string = std::regex_replace(input_as_string, std::regex(replace_regex), included_file_content, std::regex_constants::format_first_only);
+
+          ++n_included_files;
+        }
     }
   else
     {

tests/include_prm_file.prm

@@ -0,0 +1,3 @@
+# Test for including prm files. This is simply including the simpler_box.prm test.
+
+include $ASPECT_SOURCE_DIR/tests/simpler_box.prm # This is an include, and this comments tests if this line is correctly replaced.

tests/include_prm_file/original.prm

@@ -0,0 +1,73 @@
+# Test for including prm files. This is simply including the simpler_box.prm test.
+
+set Dimension = 2
+set CFL number                             = 1.0
+set End time                               = 0
+set Start time                             = 0
+set Adiabatic surface temperature          = 0
+set Surface pressure                       = 0
+set Use years in output instead of seconds = false
+set Nonlinear solver scheme                = single Advection, single Stokes
+
+subsection Boundary temperature model
+  set List of model names = box
+  set Fixed temperature boundary indicators   = 0, 1
+end
+
+subsection Gravity model
+  set Model name = vertical
+end
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1.2
+    set Y extent = 1
+    set Z extent = 1
+  end
+end
+
+subsection Initial temperature model
+  set Model name = perturbed box
+end
+
+subsection Material model
+  set Model name = simpler
+
+  subsection Simpler model
+    set Reference density             = 1
+    set Reference specific heat       = 1250
+    set Reference temperature         = 1
+    set Thermal conductivity          = 1e-6
+    set Thermal expansion coefficient = 2e-5
+    set Viscosity                     = 1
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement        = 0
+  set Initial global refinement          = 5
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = 1
+  set Zero velocity boundary indicators       = 0, 2, 3
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization, velocity statistics, basic statistics,  temperature statistics, heat flux statistics
+
+  subsection Visualization
+    set Interpolate output = false
+    set List of output variables      = material properties
+
+    subsection Material properties
+      set List of material properties = density, viscosity
+    end
+  end
+end
+
+
+set Output directory = output-include_prm_file
+

tests/include_prm_file/screen-output

@@ -0,0 +1,18 @@
+
+Number of active cells: 1,024 (on 6 levels)
+Number of degrees of freedom: 13,764 (8,450+1,089+4,225)
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving temperature system... 0 iterations.
+   Solving Stokes system (GMG)... 17+0 iterations.
+
+   Postprocessing:
+     Writing graphical output:           output-include_prm_file/solution/solution-00000
+     RMS, max velocity:                  9e-09 m/s, 3.23e-08 m/s
+     Temperature min/avg/max:            0 K, 1.035 K, 1.1 K
+     Heat fluxes through boundary parts: 1.667e-07 W, 6.239e-05 W, 0 W, 0 W
+
+Termination requested by criterion: end time
+
+
+

tests/include_prm_file/statistics

@@ -0,0 +1,22 @@
+# 1: Time step number
+# 2: Time (seconds)
+# 3: Time step size (seconds)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Iterations for temperature solver
+# 8: Iterations for Stokes solver
+# 9: Velocity iterations in Stokes preconditioner
+# 10: Schur complement iterations in Stokes preconditioner
+# 11: Visualization file name
+# 12: RMS velocity (m/s)
+# 13: Max. velocity (m/s)
+# 14: Minimal temperature (K)
+# 15: Average temperature (K)
+# 16: Maximal temperature (K)
+# 17: Average nondimensional temperature (K)
+# 18: Outward heat flux through boundary with indicator 0 ("left") (W)
+# 19: Outward heat flux through boundary with indicator 1 ("right") (W)
+# 20: Outward heat flux through boundary with indicator 2 ("bottom") (W)
+# 21: Outward heat flux through boundary with indicator 3 ("top") (W)
+0 0.000000000000e+00 0.000000000000e+00 1024 9539 4225 0 16 18 18 output-include_prm_file/solution/solution-00000 9.00076535e-09 3.23376836e-08 0.00000000e+00 1.03532014e+00 1.10000000e+00 1.03532014e+00 1.66666640e-07 6.23888889e-05 0.00000000e+00 0.00000000e+00 

tests/include_prm_file_recursive.prm

@@ -0,0 +1,5 @@
+# Test for recursively including prm files.
+# This is including the include_prm_file.prm test, which
+# itself is including simpler_box.prm.
+
+include $ASPECT_SOURCE_DIR/tests/include_prm_file.prm