Draft: Extending custom reference states to active and passive scalar fields.

.github/workflows/main.yml

@@ -82,6 +82,10 @@ jobs:
         cd "$GITHUB_WORKSPACE"
         sh ./tests/custom_reference/run_test.sh
 
+        # chi custom reference state test
+        cd "$GITHUB_WORKSPACE"
+        sh ./tests/chi_custom_reference/run_test.sh
+
         # coupled bc test
         cd "$GITHUB_WORKSPACE"
         sh ./tests/coupled_bcs/run_test.sh

doc/source/Namelist_Definitions/reference_namelist.txt

@@ -40,7 +40,7 @@
 **gravity_power**
   Specifies the value of *n* (real number) used to determine the radial variation of gravitational acceleration *g* in reference type 1, where :math:`g\propto\left(\frac{r}{r_\mathrm{max}}\right)^n`.
 **ra_constants**
-  Indicates the desired value of specified constant coefficients when reading the value from main_input instead of from a custom-refernce file.  For use with override_constants or override_constant flags.   Syntax is:
+  Indicates the desired value of specified constant coefficients when reading the value from main_input instead of from a custom-reference file.  For use with override_constants or override_constant flags.   Syntax is:
 
   ::
 
@@ -51,9 +51,9 @@
     ...
    /
 **with_custom_constants**
-  Comma separated list of integers indicating which constant coefficients should be read from a custom-refernce file when with_custom_reference is true.  
+  Comma separated list of integers indicating which constant coefficients should be read from a custom-reference file when with_custom_reference is true.  
 **with_custom_functions**
-  Comma separated list of integers indicating which non-constant coefficients should be read from a custom-refernce file when with_custom_reference is true.
+  Comma separated list of integers indicating which non-constant coefficients should be read from a custom-reference file when with_custom_reference is true.
 **with_custom_reference**  
   Logical flag that indicates some constant and non-constant coefficients should be read from a custom-reference file and used to overwrite those values otherwise assigned for reference_Types 1--3.  Default value is .false.
 **custom_reference_file**

doc/source/User_Guide/arbitrary_scalar_fields.txt

@@ -11,8 +11,8 @@ Both types of scalar fields are evolved using:
 .. math::
    :label: passive_scalar_eval
 
-   \frac{\partial \chi_{p_i}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\chi_{p_i}  + d_{1_i}\,\mathrm{g}_{1_i}(r)v_r =\
-        d_{2_i}\,\boldsymbol{\nabla}\cdot\left[\mathrm{g}_{2_i}(r)\,\boldsymbol{\nabla}\chi_{p_i} \right] + d_{3_i}\,\mathrm{g}_{3_i},\
+   \frac{\partial \chi_{p_i}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\chi_{p_i}  + d_{p1_i}\,\mathrm{g}_{p1_i}(r)v_r =\
+        d_{p2_i}\,\boldsymbol{\nabla}\cdot\left[\mathrm{g}_{p2_i}(r)\,\boldsymbol{\nabla}\chi_{p_i} \right] + d_{p3_i}\,\mathrm{g}_{p3_i},\
         \quad 1 \leq i \leq 50
 
 for the passive scalars, and:
@@ -21,21 +21,21 @@ for the passive scalars, and:
 .. math::
    :label: active_scalar_eval
 
-   \frac{\partial \chi_{a_i}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\chi_{a_i}  + d_{4_i}\,\mathrm{g}_{4_i}(r)v_r =\
-        d_{5_i}\,\boldsymbol{\nabla}\cdot\left[\mathrm{g}_{5_i}(r)\,\boldsymbol{\nabla}\chi_{a_i} \right] + d_{6_i}\,\mathrm{g}_{6_i},\
+   \frac{\partial \chi_{a_i}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\chi_{a_i}  + d_{a1_i}\,\mathrm{g}_{a1_i}(r)v_r =\
+        d_{a2_i}\,\boldsymbol{\nabla}\cdot\left[\mathrm{g}_{a2_i}(r)\,\boldsymbol{\nabla}\chi_{a_i} \right] + d_{a3_i}\,\mathrm{g}_{a3_i},\
         \quad 1 \leq i \leq 50
 
 for the active scalars.  
 
 In addition the momentum equation is modified to include a buoyancy coupling to the active scalars through
-a :math:`\sum_i d_{7_i}\,\mathrm{g}_{7_i}(r)\chi_{a_i}` term:
+a :math:`\sum_i d_{a4_i}\,\mathrm{g}_{a5_i}(r)\chi_{a_i}` term:
 
 .. math::
    :label: momentum_active_scalars
 
        \mathrm{f}_1(r)\left[\frac{\partial \boldsymbol{v}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}  %advection
         + c_1\boldsymbol{\hat{z}}\times\boldsymbol{v} \right]  =\ % Coriolis
-       &\left[c_2\,\mathrm{f}_2(r)\Theta - \sum_i d_{7_i}\,\mathrm{g}_{7_i}(r)\chi_{a_i} \right]\,\boldsymbol{\hat{r}} % buoyancy
+       &\left[c_2\,\mathrm{f}_2(r)\Theta - \sum_i d_{a4_i}\,\mathrm{g}_{a5_i}(r)\chi_{a_i} \right]\,\boldsymbol{\hat{r}} % buoyancy
         - c_3\,\mathrm{f}_1(r)\boldsymbol{\nabla}\left(\frac{P}{\mathrm{f}_1(r)}\right) % pressure
         \\
         &+ c_4\left(\boldsymbol{\nabla}\times\boldsymbol{B}\right)\times\boldsymbol{B} % Lorentz Force
@@ -46,34 +46,23 @@ where all other terms are as defined in :ref:`equations_solved`.
 Limitations
 ^^^^^^^^^^^
 
-While under development several of the functions :math:`\mathrm{g}_{j_i}` and constants :math:`d_{j_i}` are hard-coded to be 0:
-
-.. math::
-
-   \begin{aligned}
-   \mathrm{g}_{1_i}(r) &\rightarrow 0\; &d_{1_i} &\rightarrow 0 \\
-   \mathrm{g}_{3_i}(r) &\rightarrow 0\; &d_{3_i} &\rightarrow 0 \\
-   \mathrm{g}_{4_i}(r) &\rightarrow 0\; &d_{4_i} &\rightarrow 0 \\
-   \mathrm{g}_{6_i}(r) &\rightarrow 0\; &d_{6_i} &\rightarrow 0 \\
-   \end{aligned}
-
-which reduces :eq:`passive_scalar_eval` and :eq:`active_scalar_eval` to simple advection-diffusion equations.  The remaining non-zero coefficients are set according to the reference type selected.  Currently only reference_type=1 is tested but development is continuing on supporting other reference types, including custom forms for the
-functions :math:`\mathrm{g}_{j_i}` and the constants :math:`d_{j_i}`.
+The coefficients are set according to the reference type selected.  Currently only reference_type=1 is well-tested but testing is continuing on the other reference types, including custom forms for the
+functions :math:`\mathrm{g}_{sj_i}` and the constants :math:`d_{sj_i}`.
 
 Nondimensional Boussinesq Formulation
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Rayleigh can be run using a nondimensional, Boussinesq formulation
 of the MHD equations (**reference_type=1**). Adopting this nondimensionalization is equivalent to assigning the following to the
-functions :math:`\mathrm{g}_{j_i}` and the constants :math:`d_{j_i}`:
+functions :math:`\mathrm{g}_{sj_i}` and the constants :math:`d_{sj_i}`:
 
 
 .. math::
 
    \begin{aligned}
-   \mathrm{g}_{2_i}(r) &\rightarrow \tilde{\kappa}(r)\; &d_{2_i} &\rightarrow \frac{1}{Pr_{\chi_{p_i}}} \\
-   \mathrm{g}_{5_i}(r) &\rightarrow \tilde{\kappa}(r)\; &d_{5_i} &\rightarrow \frac{1}{Pr_{\chi_{a_i}}} \\
-   \mathrm{g}_{7_i}(r) &\rightarrow \left(\frac{r}{r_o}\right)^n \; &d_{7_i} &\rightarrow \frac{Ra_{\chi_{a_i}}}{Pr_{\chi_{a_i}}} \\
+   \mathrm{g}_{p2_i}(r) &\rightarrow \tilde{\kappa}_{p_i}(r)\; &d_{p2_i} &\rightarrow \frac{1}{Pr_{\chi_{p_i}}} \\
+   \mathrm{g}_{a2_i}(r) &\rightarrow \tilde{\kappa}_{a_i}(r)\; &d_{a2_i} &\rightarrow \frac{1}{Pr_{\chi_{a_i}}} \\
+   \mathrm{g}_{a5_i}(r) &\rightarrow \left(\frac{r}{r_o}\right)^n \; &d_{a4_i} &\rightarrow \frac{Ra_{\chi_{a_i}}}{Pr_{\chi_{a_i}}} \\
    \end{aligned}
 
 When these substitutions, along with those in :ref:`boussinesq`, are made, :eq:`passive_scalar_eval`-:eq:`momentum_active_scalars` transform to:
@@ -82,10 +71,10 @@ When these substitutions, along with those in :ref:`boussinesq`, are made, :eq:`
    :label: boussinesq_scalars
 
    \frac{\partial \chi_{p_i}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\chi_{p_i}  &=\
-        \frac{1}{Pr_{\chi_{p_i}}}\,\boldsymbol{\nabla}\cdot\left[\tilde{\kappa}(r)\,\boldsymbol{\nabla}\chi_{p_i} \right],\
+        \frac{1}{Pr_{\chi_{p_i}}}\,\boldsymbol{\nabla}\cdot\left[\tilde{\kappa}_{p_i}(r)\,\boldsymbol{\nabla}\chi_{p_i} \right],\
         \quad 1 \leq i \leq 50 \\
    \frac{\partial \chi_{a_i}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\chi_{a_i}  &=\
-        \frac{1}{Pr_{\chi_{a_i}}}\,\boldsymbol{\nabla}\cdot\left[\tilde{\kappa}(r)\,\boldsymbol{\nabla}\chi_{a_i} \right],\
+        \frac{1}{Pr_{\chi_{a_i}}}\,\boldsymbol{\nabla}\cdot\left[\tilde{\kappa}_{a_i}(r)\,\boldsymbol{\nabla}\chi_{a_i} \right],\
         \quad 1 \leq i \leq 50 \\
        \left[\frac{\partial \boldsymbol{v}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\boldsymbol{v}  %advection
         + \frac{2}{E}\boldsymbol{\hat{z}}\times\boldsymbol{v} \right]  &= % Coriolis
@@ -99,24 +88,36 @@ Dimensional Anelastic Formulation
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Dimensional, anelastic mode (cgs units; **reference_type=2**
-), is currently untested with arbitrary scalar fields and requires further development.
+), is currently untested with arbitrary scalar fields.
 
 
 Nondimensional Anelastic Formulation
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Non-dimensional, anelastic mode (cgs units; **reference_type=3**
-), is currently untested with arbitrary scalar fields and requires further development.  In addition to the functions set in :ref:`nondim_anelastic`, the following substitutions are made for :math:`\mathrm{g}_{j_i}`
-and :math:`d_{j_i}`:
+), is currently untested with arbitrary scalar fields.  In addition to the functions set in :ref:`nondim_anelastic`, the following substitutions are made for :math:`\mathrm{g}_{sj_i}`
+and :math:`d_{sj_i}`:
 
 .. math::
 
    \begin{aligned}
-   \mathrm{g}_{2_i}(r) &\rightarrow \tilde{\kappa}(r)\; &d_{2_i} &\rightarrow \frac{E}{Pr_{\chi_{p_i}}} \\
-   \mathrm{g}_{5_i}(r) &\rightarrow \tilde{\kappa}(r)\; &d_{5_i} &\rightarrow \frac{E}{Pr_{\chi_{a_i}}} \\
-   \mathrm{g}_{7_i}(r) &\rightarrow \tilde{\rho}(r)\frac{r_\mathrm{max}^2}{r^2}\; &c_{7_i} &\rightarrow \mathrm{Ra}_{\chi_{a_i}}^* \\
+   \mathrm{g}_{p2_i}(r) &\rightarrow \tilde{\kappa}_{p_i}(r)\; &d_{p2_i} &\rightarrow \frac{E}{Pr_{\chi_{p_i}}} \\
+   \mathrm{g}_{a2_i}(r) &\rightarrow \tilde{\kappa}_{a_i}(r)\; &d_{a2_i} &\rightarrow \frac{E}{Pr_{\chi_{a_i}}} \\
+   \mathrm{g}_{a5_i}(r) &\rightarrow \tilde{\rho}(r)\frac{r_\mathrm{max}^2}{r^2}\; &c_{a4_i} &\rightarrow \mathrm{Ra}_{\chi_{a_i}}^* \\
    \end{aligned}
 
+Custom Reference Formulation
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+With **reference_type=4**, custom reference files for active and passive scalar fields can be supplied using the namelist parameters **chi_a_custom_reference_file** and **chi_p_custom_reference_file** respectively.  These files can be constructed using the `scalar_equation_coefficients` class in `post_processing/reference_tools.py`.  
+
+The custom reference file for active scalars can be used to define :math:`d_{aj_i}` and :math:`f_{aj_i}` coefficients.  The custom reference file for passive scalars can be used to defined :math:`d_{pj_i}` and :math:`f_{pj_i}` coefficients.  If :math:`d_{p4_i}` or :math:`d_{p5_i}` are supplied they are ignored.  Additionally both files can be used to define:
+
+.. math::
+    \mathrm{g}_{t4_i}(r) = \frac{d\ln{\mathrm{g}_{t2_i}}}{dr}
+
+or this can be left undefined and it will be calculated in Rayleigh.
+
 
 .. _scalar_setup:
 
@@ -214,7 +215,7 @@ Physical controls
   Determines the radial profile of the diffusivity for active scalar i.
    * type 1 : no radial variation
    * type 2 : diffusivity profile varies as :math:`\rho^{n}` for some real number *n*.
-   * type 3 : diffusivity profile is read from a custom-reference-state file (under development)
+   * type 3 : diffusivity profile is read from a custom-reference-state file
 **kappa_chi_a_top(i)**
   Specifies the value of the diffusivity coefficient at the upper boundary for active scalar i.  This is primarily used for dimensional models or those employing a custom nondimensionalization via Rayleigh's custom-reference interface.   For Rayleigh's intrinsic nondimensional reference states, the following values are assumed:
    * reference_type 1: :math:`\kappa_\mathrm{top}=1/\mathrm{Pr_{\chi_{a_i}}}`
@@ -225,13 +226,78 @@ Physical controls
   Determines the radial profile of the diffusivity for passive scalar i.
    * type 1 : no radial variation
    * type 2 : diffusivity profile varies as :math:`\rho^{n}` for some real number *n*.
-   * type 3 : diffusivity profile is read from a custom-reference-state file (under development)
+   * type 3 : diffusivity profile is read from a custom-reference-state file
 **kappa_chi_p_top(i)**
   Specifies the value of the diffusivity coefficient at the upper boundary for passive scalar i.  This is primarily used for dimensional models or those employing a custom nondimensionalization via Rayleigh's custom-reference interface.   For Rayleigh's intrinsic nondimensional reference states, the following values are assumed:
    * reference_type 1: :math:`\kappa_\mathrm{top}=1/\mathrm{Pr_{\chi_{p_i}}}`
    * reference_type 3: :math:`\kappa_\mathrm{top}=\mathrm{Ek}/\mathrm{Pr_{\chi_{p_i}}}`
 **kappa_chi_p_power(i)**
-  Denotes the value of the exponent *n* in the :math:`\rho^{n}` variation associated with diffusion type 2 for passive scalar i.
+  Denotes the value of the exponent *n* in the :math:`\rho^{n}` variation associated with diffusion type 2 for passive scalar i
+**chi_a_advect_reference_state**
+  List of booleans that determine whether the reference-state chi is advected for each active scalar field.  The default is .true. for all active scalar fields (which makes no difference if the reference-state gradient is set to 0, the default).  When set to .false., the :math:`v_r\frac{\partial\overline{\chi}_a}{\partial r} = d_{a1_i}g_{a1,i}(r)v_r` term is omitted in the active scalar equation.
+**chi_p_advect_reference_state**
+  List of booleans that determine whether the reference-state chi is advected for each passive scalar field.  The default is .true. for all passive scalar fields (which makes no difference if the reference-state gradient is set to 0, the default).  When set to .false., the :math:`v_r\frac{\partial\overline{\chi}_p}{\partial r} = d_{p1_i}g_{p1,i}(r)v_r` term is omitted in the passive scalar equation.
+
+Reference
+*********
+
+**chi_a_custom_reference_file**
+  Name of file from which to read custom-reference-state information for all the active scalar fields when using reference_type 4 or when augmenting reference types 1--3.
+**override_chi_a_constant**
+  Indicates that particular constant coefficients for the active scalar fields, rather than all, should be overridden using main_input values when using reference_type 4.  Multiple constant overrides can be specified, one per line, with the syntax, e.g. for active scalar field 1 and constant coefficients 2 & 3:
+
+  ::
+
+   &Reference_Namelist
+    ...
+    override_chi_a_constant(1, 2) = T
+    override_chi_a_constant(1, 3) = T
+    ...
+   /
+**chi_a_constants**
+  Indicates the desired value of specified constant coefficients for the active scalar fields when reading the value from main_input instead of from a custom-reference file.  For use with override_constants or override_chi_a_constant flags.   Syntax is, e.g. for active scalar field 1 and constant coefficients 2 & 3:
+
+  ::
+
+   &Reference_Namelist
+    ...
+    chi_a_constants(1, 2) = 1.0
+    chi_a_constants(1, 3) = 14.0
+    ...
+   /
+**with_custom_chi_a_constants**
+  Comma separated list of integers for a particular active scalar field indicating which constant coefficients should be read from a custom-reference file when with_custom_reference is true.
+**with_custom_chi_a_functions**
+  Comma separated list of integers for a particular active scalar field indicating which non-constant coefficients should be read from a custom-reference file when with_custom_reference is true.
+**chi_p_custom_reference_file**
+  Name of file from which to read custom-reference-state information for all the passive scalar fields when using reference_type 4 or when augmenting reference types 1--3.
+**override_chi_p_constant**
+  Indicates that particular constant coefficients for the passive scalar fields, rather than all, should be overridden using main_input values when using reference_type 4.  Multiple constant overrides can be specified, one per line, with the syntax, e.g. for passive scalar field 1 and constant coefficients 2 & 3:
+
+  ::
+
+   &Reference_Namelist
+    ...
+    override_chi_p_constant(1, 2) = T
+    override_chi_p_constant(1, 3) = T
+    ...
+   /
+**chi_p_constants**
+  Indicates the desired value of specified constant coefficients for the passive scalar fields when reading the value from main_input instead of from a custom-reference file.  For use with override_constants or override_chi_p_constant flags.   Syntax is, e.g. for passive scalar field 1 and constant coefficients 2 & 3:
+
+  ::
+
+   &Reference_Namelist
+    ...
+    chi_p_constants(1, 2) = 1.0
+    chi_p_constants(1, 3) = 14.0
+    ...
+   /
+**with_custom_chi_p_constants**
+  Comma separated list of integers for a particular passive scalar field indicating which constant coefficients should be read from a custom-reference file when with_custom_reference is true.
+**with_custom_chi_p_functions**
+  Comma separated list of integers for a particular passive scalar field indicating which non-constant coefficients should be read from a custom-reference file when with_custom_reference is true.
+
 
 Output Quantity Codes
 ^^^^^^^^^^^^^^^^^^^^^

examples/custom_reference_states/Anelastic_Dim_CZ.ipynb

@@ -106,7 +106,7 @@
     "import os, sys\n",
     "sys.path.insert(0, os.path.abspath('../../'))\n",
     "\n",
-    "import post_processing.reference_tools as rt # You will need the refernce_tools.py to run this notebook\n",
+    "import post_processing.reference_tools as rt # You will need the reference_tools.py to run this notebook\n",
     "import post_processing.rayleigh_diagnostics as rdiag"
    ]
   },

examples/custom_reference_states/Anelastic_NonDim_CZ.ipynb

@@ -105,7 +105,7 @@
     "import os, sys\n",
     "sys.path.insert(0, os.path.abspath('../../'))\n",
     "\n",
-    "import post_processing.reference_tools as rt # You will need the refernce_tools.py to run this notebook\n",
+    "import post_processing.reference_tools as rt # You will need the reference_tools.py to run this notebook\n",
     "import post_processing.rayleigh_diagnostics as rdiag"
    ]
   },

examples/custom_reference_states/Boussinesq_CZ.ipynb

@@ -110,7 +110,7 @@
     "import os, sys\n",
     "sys.path.insert(0, os.path.abspath('../../'))\n",
     "\n",
-    "import post_processing.reference_tools as rt # You will need the refernce_tools.py to run this notebook\n"
+    "import post_processing.reference_tools as rt # You will need the reference_tools.py to run this notebook\n"
    ]
   },
   {