diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 5164febe..f0fdfa19 100644 --- a/.github/workflows/main.yml +++ b/.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 diff --git a/CHANGELOG.md b/CHANGELOG.md index e83f4187..8fa00c4a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ### Added +- Extend custom references states to active and passive scalar fields by adding a new custom reference state file type (with an accompanying python preprocessing class and namelist options) for scalar fields. This allows the scalar equations to include non-zero (and potentially radially dependent) reference state advection and source terms. \[Cian Wilson; 7-13-2026; [#586](https://github.com/geodynamics/Rayleigh/pull/586)\] + ### Changed - Containers for the TACC systems are no longer supported. \[Philipp Edelmann; 7-12-2026; [#608](https://github.com/geodynamics/Rayleigh/pull/608)\] diff --git a/doc/source/Namelist_Definitions/reference_namelist.txt b/doc/source/Namelist_Definitions/reference_namelist.txt index 07b9888d..9888aa32 100644 --- a/doc/source/Namelist_Definitions/reference_namelist.txt +++ b/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** diff --git a/doc/source/User_Guide/arbitrary_scalar_fields.txt b/doc/source/User_Guide/arbitrary_scalar_fields.txt index de6e56de..73238123 100644 --- a/doc/source/User_Guide/arbitrary_scalar_fields.txt +++ b/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 ^^^^^^^^^^^^^^^^^^^^^ diff --git a/examples/custom_reference_states/Anelastic_Dim_CZ.ipynb b/examples/custom_reference_states/Anelastic_Dim_CZ.ipynb index cc6587e0..4f5ee4f2 100644 --- a/examples/custom_reference_states/Anelastic_Dim_CZ.ipynb +++ b/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" ] }, diff --git a/examples/custom_reference_states/Anelastic_NonDim_CZ.ipynb b/examples/custom_reference_states/Anelastic_NonDim_CZ.ipynb index b47e0188..3e9cf461 100644 --- a/examples/custom_reference_states/Anelastic_NonDim_CZ.ipynb +++ b/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" ] }, diff --git a/examples/custom_reference_states/Boussinesq_CZ.ipynb b/examples/custom_reference_states/Boussinesq_CZ.ipynb index 9c5c3a31..c4352a8f 100644 --- a/examples/custom_reference_states/Boussinesq_CZ.ipynb +++ b/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" ] }, { diff --git a/post_processing/reference_tools.py b/post_processing/reference_tools.py index 1f9e5df4..7e4a3b3a 100644 --- a/post_processing/reference_tools.py +++ b/post_processing/reference_tools.py @@ -130,6 +130,126 @@ def read(self, filename='equation_coefficients', override_nconst=False): print("Version of input file was 1.") print("Converting current equation_coefficients instance's version to %i." %class_version) +class scalar_equation_coefficients: + """ scalar equation coeff class """ + consts = ['dchirefdr_scale', 'kappa_chi_scale', 'source_chi_scale', 'buoy_chi_scale'] + funcs = ['dchirefdr', 'kappa_chi', 'source_chi', 'd_ln_kappa_chi', 'buoy_chi'] + nconst = len(consts) + nfunc = len(funcs) + version = 1 + c_dict = {name:i for i, name in enumerate(consts, start=1)} + f_dict = {name:i for i, name in enumerate(funcs, start=1)} + + def __init__(self, nscalars=1, radius=[], file=None): + """Reads/writes equation_coefficient and custom reference state files. + + Keyword arguments: + nscalars: The number of scalar fields for which to create the reference state. + This is only used if radius is provided and a new state is being created and defaults to 1. + radius: If provided, a new state will be created with these + radii and zero coefficients. + file: If provided, the reference state will be read from the + provided file. + """ + if len(radius) != 0: + if file is not None: + raise RuntimeError("Cannot provide radius and file at the same time.") + self.nscalars = nscalars + nr = len(radius) + self.nr = nr + self.radius = numpy.asarray(radius, dtype='float64') + self.functions = numpy.zeros((self.nscalars, self.nfunc, nr) , dtype='float64' ) + + self.constants = numpy.zeros((self.nscalars, self.nconst) , dtype='float64' ) + self.cset = numpy.zeros((self.nscalars, self.nconst) , dtype='int32' ) + self.fset = numpy.zeros((self.nscalars, self.nfunc) , dtype='int32' ) + elif file is not None: + self.read(filename=file) + + def __getattr__(self, nameindex): + if isinstance(nameindex, (list, tuple)): + name, index = nameindex + else: + name = nameindex + index = 0 + + if name in self.f_dict: + return self.functions[index,self.f_dict[name] - 1,:] + elif name in self.c_dict: + return self.constants[index,self.c_dict[name] - 1] + else: + raise AttributeError("'{}' has no attribute '{}'".format(self.__class__, name)) + + def __setattr__(self, nameindex, value): + if isinstance(nameindex, (list, tuple)): + name, index = nameindex + else: + name = nameindex + index = 0 + + if name in self.f_dict: + self.set_function(value, name, index=index) + elif name in self.c_dict: + self.set_constant(value, name, index=index) + else: + super().__setattr__(name, value) + + def set_function(self,y,f_name,index=0): + + if isinstance(f_name,str): + fi = self.f_dict[f_name] + else: + fi = f_name + + self.functions[index,fi-1,:] = y + self.fset[index,fi-1] = 1 + + def set_constant(self,c,c_name,index=0): + if isinstance(c_name,str): + ci = self.c_dict[c_name] + else: + ci = c_name + self.constants[index,ci-1] = c + self.cset[index,ci-1] = 1 + + def write(self, filename='scoefs.dat'): + pi = numpy.array([314],dtype='int32') + nr = numpy.array([self.nr],dtype='int32') + version = numpy.array([self.version],dtype='int32') + nconst = numpy.array([self.nconst],dtype='int32') + nfunc = numpy.array([self.nfunc],dtype='int32') + nscalars = numpy.array([self.nscalars],dtype='int32') + fd = open(filename,'wb') + pi.tofile(fd) + version.tofile(fd) + nconst.tofile(fd) + nfunc.tofile(fd) + nscalars.tofile(fd) + self.cset.tofile(fd) + self.fset.tofile(fd) + self.constants.tofile(fd) + nr.tofile(fd) + self.radius.tofile(fd) + self.functions.tofile(fd) + fd.close() + + def read(self, filename='equation_coefficients', override_nconst=False): + fd = open(filename,'rb') + picheck = numpy.fromfile(fd,dtype='int32',count=1)[0] + + self.version = numpy.fromfile(fd,dtype='int32', count=1)[0] + self.nconst = numpy.fromfile(fd,dtype='int32', count=1)[0] + self.nfunc = numpy.fromfile(fd,dtype='int32', count=1)[0] + self.nscalars = numpy.fromfile(fd,dtype='int32', count=1)[0] + self.cset = numpy.fromfile(fd,dtype='int32',count=self.nscalars*self.nconst).reshape(self.nscalars,self.nconst) + self.fset = numpy.fromfile(fd,dtype='int32',count=self.nscalars*self.nfunc).reshape(self.nscalars,self.nfunc) + self.constants = numpy.fromfile(fd,dtype='float64',count=self.nscalars*self.nconst).reshape(self.nscalars,self.nconst) + self.nr = numpy.fromfile(fd,dtype='int32',count=1)[0] + self.radius = numpy.fromfile(fd,dtype='float64',count=self.nr) + functions=numpy.fromfile(fd,dtype='float64',count=self.nscalars*self.nr*self.nfunc) + self.functions = numpy.reshape(functions, (self.nscalars,self.nfunc,self.nr)) + fd.close() + class background_state: nr = None radius = None diff --git a/pre_processing/rayleigh_spectral_input.py b/pre_processing/rayleigh_spectral_input.py index 8779dacf..735e8530 100755 --- a/pre_processing/rayleigh_spectral_input.py +++ b/pre_processing/rayleigh_spectral_input.py @@ -64,18 +64,20 @@ def dealias_g2m(g): if m_max < 0: m_max = 0 return m_max -def radial_extents(rmin=None, rmax=None, aspect_ratio=None, shell_depth=None): +def radial_extents(rmin=None, rmax=None, aspect_ratio=None, shell_depth=None, throw=True): """ Return rmin and rmax if they aren't set using aspect_ratio and shell_depth. """ if rmax is None: if shell_depth is None or aspect_ratio is None: - raise Exception("Must supply shell_depth and aspect_ratio if rmax is not set.") - rmax = shell_depth/(1.-aspect_ratio) + if throw: raise Exception("Must supply shell_depth and aspect_ratio if rmax is not set.") + else: + rmax = shell_depth/(1.-aspect_ratio) if rmin is None: - if aspect_ratio is None: - raise Exception("Must supply aspect_ratio if rmin is not set.") - rmin = rmax*aspect_ratio + if rmax is None or aspect_ratio is None: + if throw: raise Exception("Must supply aspect_ratio and rmax or shell_depth if rmin is not set.") + else: + rmin = rmax*aspect_ratio return rmin, rmax def swapwrite(vals,fd,byteswap=False): @@ -349,11 +351,7 @@ def transform_from_rtp_function(self, func, func_kwargs={}, \ raise Exception("'{}' supplied in func_kwargs. You probably didn't mean to do this.".format((k,))) # some sanity checks that we have enough info for the radial domain - if func_radius: - rmin, rmax = radial_extents(rmin, rmax, aspect_ratio, shell_depth) - else: - rmin = None - rmax = None + rmin, rmax = radial_extents(rmin, rmax, aspect_ratio, shell_depth, throw=func_radius) # work out the number of points... if func_theta or func_phi: @@ -685,8 +683,8 @@ def exec_then_eval(theta=None, phi=None, radius=None, \ if expr.find('phi') >= 0: func_args.append("phi") if expr.find('radius') >= 0: func_args.append("radius") - if func_kwargs['rmin'] is None or func_kwargs['rmax'] is None: - func_kwargs['rmin'], func_kwargs['rmax'] = radial_extents(rmin, rmax, aspect_ratio, shell_depth) + if func_kwargs['rmin'] is None or func_kwargs['rmax'] is None: + func_kwargs['rmin'], func_kwargs['rmax'] = radial_extents(rmin, rmax, aspect_ratio, shell_depth, throw=expr.find('radius')>=0) func_argstr, exec_argstr = "", "" if len(func_args) > 0: func_argstr = ", ".join([fa+"=None" for fa in func_args])+", " @@ -758,7 +756,7 @@ def parse_mode(modein): #-------------------------------------------------------------------------------- # )+ end of grouping # - from http://rick.measham.id.au/paste/explain.pl - r = re.compile(r'(?:[^,; '+os.linesep+'([{]|\([^)]*\)|\[[^]]*\]|\{[^}]*\})+') + r = re.compile(r'(?:[^,; '+os.linesep+'([{]|\\([^)]*\\)|\\[[^]]*\\]|\\{[^}]*\\})+') modelist = r.findall(modestr) for index in modelist[:min(3,len(modelist)-1)]: if not index.isdigit(): diff --git a/src/Physics/BoundaryConditions.F90 b/src/Physics/BoundaryConditions.F90 index 7a19ed83..569d91b0 100755 --- a/src/Physics/BoundaryConditions.F90 +++ b/src/Physics/BoundaryConditions.F90 @@ -29,7 +29,7 @@ Module BoundaryConditions Use Generic_Input Use PDE_Coefficients Use Checkpointing, Only : Store_BC_Mask, Load_BC_Mask - Use Controls, Only : full_restart + Use Controls, Only : full_restart, n_scalar_max Implicit None diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 6edaee7a..c8ead3c3 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -69,11 +69,17 @@ Module Controls Logical :: advect_reference_state = .true. ! Set to true to advect the reference state temperature or entropy ! This has no effect for adiabatic reference states. ! Generally only do this if reference state is nonadiabatic + ! Nondimensional variables for the active/passive scalar fields + + ! Allow up to 50 active/passive scalar fields + Integer, Parameter :: n_scalar_max = 50 Logical :: compressible = .false. !run compressible or not Integer :: n_active_scalars = 0 ! number of active scalar fields Integer :: n_passive_scalars = 0 ! number of passive scalar fields + Logical :: chi_a_advect_reference_state(1:n_scalar_max) = .true. + Logical :: chi_p_advect_reference_state(1:n_scalar_max) = .true. ! --- This flag determines if the code is run in benchmark mode ! 0 (default) is no benchmarking. 1-5 are various accuracy benchmarks (see documentation) @@ -99,7 +105,8 @@ Module Controls & momentum_advection, inertia, coriolis, centrifugal, n_active_scalars, n_passive_scalars, & & newtonian_cooling, newtonian_cooling_type, newtonian_cooling_time, & & newtonian_cooling_tvar_amp, newtonian_cooling_profile_file, & - & pseudo_incompressible, compressible, pulse_freq, pulse_sharpness + & pseudo_incompressible, compressible, pulse_freq, pulse_sharpness, & + & chi_a_advect_reference_state, chi_p_advect_reference_state !/////////////////////////////////////////////////////////////////////////// ! Temporal Controls @@ -241,6 +248,8 @@ Subroutine Restore_Physics_Defaults() pseudo_incompressible = .false. compressible = .false. advect_reference_state = .true. + chi_a_advect_reference_state = .true. + chi_p_advect_reference_state = .true. benchmark_mode = 0 benchmark_integration_interval = -1 benchmark_report_interval = -1 diff --git a/src/Physics/Initial_Conditions.F90 b/src/Physics/Initial_Conditions.F90 index 28c8f613..77afca01 100755 --- a/src/Physics/Initial_Conditions.F90 +++ b/src/Physics/Initial_Conditions.F90 @@ -31,7 +31,7 @@ Module Initial_Conditions Use Controls Use Timers Use General_MPI, Only : BCAST2D - Use PDE_Coefficients, Only : n_scalar_max, s_conductive, heating_type,ref, kappa, dlnkappa + Use PDE_Coefficients, Only : s_conductive, heating_type,ref, kappa, dlnkappa Use BoundaryConditions, Only : T_top, T_bottom, fix_tvar_Top, fix_tvar_bottom,& & fix_dtdr_top, fix_dtdr_bottom, dtdr_top, dtdr_bottom, & & C10_bottom, C11_bottom, C1m1_bottom diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 5502cffc..72421235 100644 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -58,8 +58,12 @@ Module PDE_Coefficients Real*8, Allocatable :: dsdr(:) Real*8, Allocatable :: dsdr_over_cp(:) ! (1/c_P) ds/dr Real*8, Allocatable :: d2s_over_cp(:) ! (1/c_P) d2s/dr2 + Real*8, Allocatable :: dchirefadr(:,:) + Real*8, Allocatable :: dchirefpdr(:,:) Real*8, Allocatable :: heating(:) + Real*8, Allocatable :: chi_a_source(:,:) + Real*8, Allocatable :: chi_p_source(:,:) Real*8 :: Centrifugal_Coeff ! Multiplies in the s term in compressible eq Real*8 :: Coriolis_Coeff ! Multiplies z_hat x u in momentum eq. @@ -76,8 +80,6 @@ Module PDE_Coefficients End Type ReferenceInfo Type(ReferenceInfo) :: ref - ! Allow up to 50 active/passive scalar fields - Integer, Parameter :: n_scalar_max = 50 ! Version number for the "equation_coefficients" container that is output to the simulation directory ! (i.e., the human-obtainable version of "ref") @@ -112,23 +114,40 @@ Module PDE_Coefficients Real*8 :: Angular_Velocity = -1.0d0 ! Frame rotation rate (sets Coriolis force) ! Custom reference-state variables (reference_type = 4) - Integer, Parameter :: max_ra_constants = 11 + 2*n_scalar_max - Integer, Parameter :: max_ra_functions = 14 + 2*n_scalar_max - Integer :: n_ra_constants - Integer :: n_ra_functions + Integer, Parameter :: max_ra_constants = 11, max_chi_constants = 4 + Integer, Parameter :: max_ra_functions = 14, max_chi_functions = 5 + Integer :: n_ra_constants, n_chi_constants + Integer :: n_ra_functions, n_chi_functions Logical :: with_custom_reference = .false. Logical :: override_constants = .false. Logical :: override_constant(1:max_ra_constants) = .false. ! in namelist + Logical, target :: override_chi_a_constant(1:n_scalar_max, 1:max_chi_constants) = .false. ! in namelist + Logical, target :: override_chi_p_constant(1:n_scalar_max, 1:max_chi_constants) = .false. ! in namelist Integer :: with_custom_constants(1:max_ra_constants) = 0 ! in namelist Integer :: with_custom_functions(1:max_ra_functions) = 0 ! in namelist + Integer, target :: with_custom_chi_a_constants(1:n_scalar_max, 1:max_chi_constants) = 0 ! in namelist + Integer, target :: with_custom_chi_a_functions(1:n_scalar_max, 1:max_chi_functions) = 0 ! in namelist + Integer, target :: with_custom_chi_p_constants(1:n_scalar_max, 1:max_chi_constants) = 0 ! in namelist + Integer, target :: with_custom_chi_p_functions(1:n_scalar_max, 1:max_chi_functions) = 0 ! in namelist Real*8 :: ra_constants(1:max_ra_constants) = 0.0d0 ! in namelist + Real*8, target :: chi_a_constants(1:n_scalar_max, 1:max_chi_constants) = 0.0d0 ! in namelist + Real*8, target :: chi_p_constants(1:n_scalar_max, 1:max_chi_constants) = 0.0d0 ! in namelist Integer, Allocatable :: ra_constant_set(:) Integer, Allocatable :: ra_function_set(:) + Integer, Allocatable, target :: chi_a_constant_set(:,:), chi_p_constant_set(:,:) + Integer, Allocatable, target :: chi_a_function_set(:,:), chi_p_function_set(:,:) Logical, Allocatable :: use_custom_constant(:) Logical, Allocatable :: use_custom_function(:) + Logical, Allocatable, target :: chi_a_use_custom_constant(:,:), chi_p_use_custom_constant(:,:) + Logical, Allocatable, target :: chi_a_use_custom_function(:,:), chi_p_use_custom_function(:,:) Real*8, Allocatable :: ra_functions(:,:) + Real*8, Allocatable, target :: chi_a_functions(:,:,:), chi_p_functions(:,:,:) Logical :: custom_reference_read = .false. - Character*120 :: custom_reference_file ='nothing' + Logical :: chi_a_custom_reference_read = .false. + Logical :: chi_p_custom_reference_read = .false. + Character*120 :: custom_reference_file ='nothing' + Character*120 :: chi_a_custom_reference_file = 'nothing' + Character*120 :: chi_p_custom_reference_file = 'nothing' ! General nondimensional anelastic polytrope variables (reference_type = 5) Logical :: ND_Volume_Average = .true. @@ -174,8 +193,13 @@ Module PDE_Coefficients & chi_a_convective_rossby_number, chi_p_prandtl_number, & & ND_Volume_Average, ND_Inner_Radius, ND_Outer_Radius, & & Assume_Flux_Ra, Specific_Heat_Ratio, Buoyancy_Number_Visc,& - & Buoyancy_Number_Rot, Sigma_Parameter, Length_Scale, & - & bigz, gas_gamma + & Buoyancy_Number_Rot, Sigma_Parameter, Length_Scale,& + & bigz, gas_gamma,& + & chi_a_custom_reference_file, chi_p_custom_reference_file,& + & override_chi_a_constant, override_chi_p_constant,& + & chi_a_constants, chi_p_constants,& + & with_custom_chi_a_constants, with_custom_chi_p_constants,& + & with_custom_chi_a_functions, with_custom_chi_p_functions !/////////////////////////////////////////////////////////////////////////////////////// ! II. Variables Related to the Transport Coefficients @@ -256,7 +280,11 @@ Subroutine Initialize_Reference() Call stdout%print(" ") Endif - If (with_custom_reference) Call Augment_Reference() + If (with_custom_reference) Then + Call Augment_Reference() + if (n_active_scalars > 0) Call Augment_Scalar_Reference(.true.) + if (n_passive_scalars > 0) Call Augment_Scalar_reference(.false.) + Endif If (rotation) Call Set_Rotation_dt() @@ -267,8 +295,11 @@ End Subroutine Initialize_Reference Subroutine Allocate_Reference_State Implicit None - n_ra_constants = 11 + 2*(n_active_scalars + n_passive_scalars) - n_ra_functions = 14 + 2*(n_active_scalars + n_passive_scalars) + n_ra_constants = 11 + n_ra_functions = 14 + + n_chi_constants = 4 + n_chi_functions = 5 Allocate(ref%density(1:N_R)) Allocate(ref%temperature(1:N_R)) @@ -286,7 +317,11 @@ Subroutine Allocate_Reference_State Allocate(ref%ohmic_amp(1:N_R)) Allocate(ref%viscous_amp(1:N_R)) Allocate(ref%heating(1:N_R)) - Allocate(ref%chi_buoyancy_coeff(n_active_scalars,1:N_R)) + Allocate(ref%chi_buoyancy_coeff(1:N_R,n_active_scalars)) + Allocate(ref%chi_a_source(1:N_R,n_active_scalars)) + Allocate(ref%chi_p_source(1:N_R,n_passive_scalars)) + Allocate(ref%dchirefadr(1:N_R,n_active_scalars)) + Allocate(ref%dchirefpdr(1:N_R,n_passive_scalars)) Allocate(ra_constant_set(1:n_ra_constants)) ra_constant_set = 0 @@ -299,6 +334,28 @@ Subroutine Allocate_Reference_State Allocate(ra_functions(1:N_R, 1:n_ra_functions)) ra_functions(:,:) = Zero + Allocate(chi_a_constant_set(1:n_chi_constants,1:n_active_scalars)) + chi_a_constant_set = 0 + Allocate(chi_a_function_set(1:n_chi_functions,1:n_active_scalars)) + chi_a_function_set = 0 + Allocate(chi_a_use_custom_constant(1:n_chi_constants,1:n_active_scalars)) + chi_a_use_custom_constant = .false. + Allocate(chi_a_use_custom_function(1:n_chi_functions,1:n_active_scalars)) + chi_a_use_custom_function = .false. + Allocate(chi_a_functions(1:n_r, 1:n_chi_functions, 1:n_active_scalars)) + chi_a_functions(:,:,:) = Zero + + Allocate(chi_p_constant_set(1:n_chi_constants,1:n_passive_scalars)) + chi_p_constant_set = 0 + Allocate(chi_p_function_set(1:n_chi_functions,1:n_passive_scalars)) + chi_p_function_set = 0 + Allocate(chi_p_use_custom_constant(1:n_chi_constants,1:n_passive_scalars)) + chi_p_use_custom_constant = .false. + Allocate(chi_p_use_custom_function(1:n_chi_functions,1:n_passive_scalars)) + chi_p_use_custom_function = .false. + Allocate(chi_p_functions(1:n_r, 1:n_chi_functions, 1:n_passive_scalars)) + chi_p_functions(:,:,:) = Zero + ref%density(:) = Zero ref%temperature(:) = Zero ref%dlnrho(:) = Zero @@ -345,6 +402,7 @@ Subroutine Constant_Reference() Integer :: i,j Real*8 :: r_outer, r_inner, prefactor, amp, pscaling Character*12 :: dstring + Character(len=2) :: sind Character*8 :: dofmt = '(ES12.5)' @@ -374,6 +432,10 @@ Subroutine Constant_Reference() ref%dsdr = 0.0d0 ref%dsdr_over_cp = 0.0d0 ref%d2s_over_cp = 0.0d0 + ref%dchirefadr = 0.0d0 + ref%dchirefpdr = 0.0d0 + ref%chi_a_source = 0.0d0 + ref%chi_p_source = 0.0d0 amp = Rayleigh_Number/Prandtl_Number @@ -382,11 +444,18 @@ Subroutine Constant_Reference() Enddo do j = 1, n_active_scalars - amp = -chi_a_Rayleigh_Number(j)/chi_a_Prandtl_Number(j) + If (my_rank .eq. 0) Then + Write(sind,'(I2)') j + Write(dstring,dofmt)chi_a_Rayleigh_Number(j) + Call stdout%print(" ---- Chi Rayleigh Number "//Adjustl(sind)//" : "//trim(dstring)) + Write(dstring,dofmt)chi_a_Prandtl_Number(j) + Call stdout%print(" ---- Chi Prandtl Number "//Adjustl(sind)//" : "//trim(dstring)) + Endif + amp = -chi_a_Rayleigh_Number(j)/chi_a_Prandtl_Number(j) - Do i = 1, N_R - ref%chi_buoyancy_coeff(j,i) = amp*(radius(i)/radius(1))**gravity_power - Enddo + Do i = 1, N_R + ref%chi_buoyancy_coeff(i,j) = amp*(radius(i)/radius(1))**gravity_power + Enddo enddo pressure_specific_heat = 1.0d0 @@ -438,7 +507,8 @@ Subroutine Constant_Reference() ra_functions(:,2) = (radius(:)/radius(1))**gravity_power ra_constants(2) = Rayleigh_Number/Prandtl_Number Do i = 1, n_active_scalars - ra_constants(12+(i-1)*2) = -chi_a_Rayleigh_Number(i)/chi_a_Prandtl_Number(i) + chi_a_constants(i,4) = -chi_a_Rayleigh_Number(i)/chi_a_Prandtl_Number(i) + chi_a_functions(:,5,i) = (radius(:)/radius(1))**gravity_power Enddo End Subroutine Constant_Reference @@ -482,7 +552,7 @@ Subroutine Polytropic_ReferenceND() gravity = (rmax**2)*OneOverRSquared(:) ref%Buoyancy_Coeff = gravity*Modified_Rayleigh_Number*ref%density do i = 1, n_active_scalars - ref%chi_buoyancy_coeff(i,:) = -gravity*chi_a_modified_rayleigh_number(i)*ref%density + ref%chi_buoyancy_coeff(:,i) = -gravity*chi_a_modified_rayleigh_number(i)*ref%density enddo !Compute the background temperature gradient : dTdr = -Dg, d2Tdr2 = 2*D*g/r (for g ~1/r^2) @@ -538,7 +608,8 @@ Subroutine Polytropic_ReferenceND() ra_functions(:,2) = gravity*ref%density ra_constants(2) = Modified_Rayleigh_Number Do i = 1, n_active_scalars - ra_constants(12+(i-1)*2) = -chi_a_modified_rayleigh_number(i) + chi_a_constants(i,4) = -chi_a_modified_rayleigh_number(i) + chi_a_functions(:,5,i) = gravity*ref%density Enddo DeAllocate(dtmparr, gravity) @@ -552,6 +623,7 @@ Subroutine Polytropic_ReferenceND_General() Real*8, Allocatable :: gravity(:), flux_nonrad(:), partial_heating(:), nsquared(:), & & zeta(:), dzeta(:), d2zeta(:), dlnzeta(:), d2lnzeta(:) Character*12 :: dstring + Character(len=2) :: sind Character*8 :: dofmt = '(ES12.5)' Logical :: ND_Length_Shell_Depth = .false. ! Defaults to .true. under following logic Logical :: Adiabatic_Polytrope ! To be determined (to within the tolerance tol) @@ -855,7 +927,14 @@ Subroutine Polytropic_ReferenceND_General() ref%Buoyancy_Coeff(:) = (Rayleigh_Number/Prandtl_Number)*ref%density(:)*gravity(:) Do i = 1, n_active_scalars - ref%Chi_Buoyancy_Coeff(i,:) = -(Chi_A_Rayleigh_Number(i)/Chi_A_Prandtl_Number(i))*& + If (my_rank .eq. 0) Then + Write(sind,'(I2)') i + Write(dstring,dofmt)chi_a_Rayleigh_Number(i) + Call stdout%print(" ---- Chi Rayleigh Number "//Adjustl(sind)//" : "//trim(dstring)) + Write(dstring,dofmt)chi_a_Prandtl_Number(i) + Call stdout%print(" ---- Chi Prandtl Number "//Adjustl(sind)//" : "//trim(dstring)) + Endif + ref%Chi_Buoyancy_Coeff(:,i) = -(Chi_A_Rayleigh_Number(i)/Chi_A_Prandtl_Number(i))*& & ref%density(:)*gravity(:) Enddo @@ -923,9 +1002,9 @@ Subroutine Polytropic_ReferenceND_General() ! Set the buoyancy constants / functions ra_functions(:,2) = gravity*ref%density ra_constants(2) = Rayleigh_Number/Prandtl_Number - Do i = 1, n_active_scalars - ra_constants(12+(i-1)*2) = -Chi_A_Rayleigh_Number(i)/Chi_A_Prandtl_Number(i) + chi_a_constants(i,4) = -(Chi_A_Rayleigh_Number(i)/Chi_A_Prandtl_Number(i)) + chi_a_functions(:,5,i) = ref%density(:)*gravity(:) Enddo DeAllocate(gravity, zeta, dzeta, d2zeta, dlnzeta, d2lnzeta, flux_nonrad, partial_heating) @@ -1030,7 +1109,7 @@ Subroutine Polytropic_Reference() Ref%Buoyancy_Coeff = gravity/Pressure_Specific_Heat*ref%density do i = 1, n_active_scalars - ref%chi_buoyancy_coeff(i,:) = -gravity/pressure_specific_heat*ref%density + ref%chi_buoyancy_coeff(:,i) = -gravity/pressure_specific_heat*ref%density end do Deallocate(zeta, gravity) @@ -1054,8 +1133,9 @@ Subroutine Polytropic_Reference() ra_functions(:,2) = ref%Buoyancy_Coeff ra_constants(2) = 1.0d0 Do i = 1, n_active_scalars - ra_constants(12+(i-1)*2) = -1.0d0 - Enddo + chi_a_constants(i,4) = 1.0d0 + chi_a_functions(:,5,i) = ref%chi_buoyancy_coeff(:,i) + Enddo End Subroutine Polytropic_Reference @@ -1200,6 +1280,146 @@ Subroutine Augment_Reference() End Subroutine Augment_Reference + Subroutine Augment_Scalar_Reference(active) + Implicit None + Logical, intent(in) :: active + Real*8, Allocatable :: temp_functions(:,:,:), temp_constants(:,:) + Character(len=2) :: intstring + Character*12 :: dstring + Character*8 :: dofmt = '(ES12.5)' + Character(len=2) :: sind + ! variables to point at different global parameters depending on whether we are reading in active or passive scalar info + Real*8, pointer :: chi_l_constants(:,:), chi_l_functions(:,:,:) + Logical, pointer :: chi_l_use_custom_constant(:,:), chi_l_use_custom_function(:,:) + Integer, pointer :: chi_l_constant_set(:,:), chi_l_function_set(:,:) + Character*120 :: chi_l_custom_reference_file + Character*7 :: chi_l_type_str + Integer :: i, n_l_scalars + + if (active) then + chi_l_constant_set => chi_a_constant_set + chi_l_function_set => chi_a_function_set + chi_l_constants => chi_a_constants + chi_l_functions => chi_a_functions + chi_l_use_custom_constant => chi_a_use_custom_constant + chi_l_use_custom_function => chi_a_use_custom_function + chi_l_custom_reference_file = chi_a_custom_reference_file + n_l_scalars = n_active_scalars + chi_l_type_str = 'active' + else + chi_l_constant_set => chi_p_constant_set + chi_l_function_set => chi_p_function_set + chi_l_constants => chi_p_constants + chi_l_functions => chi_p_functions + chi_l_use_custom_constant => chi_p_use_custom_constant + chi_l_use_custom_function => chi_p_use_custom_function + chi_l_custom_reference_file = chi_p_custom_reference_file + n_l_scalars = n_passive_scalars + chi_l_type_str = 'passive' + endif + + If (my_rank .eq. 0) Then + if (active) Then + Call stdout%print('Active scalar field reference state will be augmented.') + Call stdout%print('Only source, buoyancy, or background dchidr may be modified.') + else + Call stdout%print('Passive scalar field reference state will be augmented.') + Call stdout%print('Only source or background dchidr may be modified.') + endif + Call stdout%print('Source requires both d_3 and g_3 to be set.') + if (active) Call stdout%print('Buoyancy requires both d_4 and g_5 to be set.') + Call stdout%print('Reference dchidr requires g_1 to be set (sets d_1 to 1 if unspecified).') + Call stdout%print('Reading '//Adjustl(chi_l_type_str)//' scalars from: '//TRIM(ADJUSTL(chi_l_custom_reference_file)) ) + Endif + + ! Before reading the custom reference file, guard against + ! potential overwrites of all chi_functions and chi_constants which, + ! in this case were determined by ref_types 1,2,3. Only certain + ! combinations are allowed to be overwritten. + + Allocate(temp_functions(1:n_r, 1:n_chi_functions, 1:n_l_scalars)) + Allocate(temp_constants(1:n_l_scalars, 1:n_chi_constants)) + + temp_functions(:,:,:) = chi_l_functions(:,:,:) + ! Note that chi_a_constants is allocated up to n_scalar_max x max_chi_constants, + ! which could be more than n_active_scalars xn_chi_constants + temp_constants(:,:) = chi_l_constants(1:n_l_scalars, 1:n_chi_constants) + + Call Read_Scalar_Custom_Reference_File(chi_l_custom_reference_file, active) + + Do i = 1, n_l_scalars + If (my_rank .eq. 0) Write(sind,'(I2)') i + + If (chi_l_use_custom_constant(3,i) .and. chi_l_use_custom_function(3,i)) Then + If (my_rank .eq. 0) Then + Call stdout%print('Chi source for '//Adjustl(chi_l_type_str)// & + ' scalar field '//Adjustl(sind)//' has been set to:') + Call stdout%print('d_3*g_3') + Call stdout%print(' ') + Endif + ! can't declare items in type as a target so have to do this if statement instead + if (active) then + ref%chi_a_source(:,i) = chi_l_constants(i,3)*chi_l_functions(:,3,i) + else + ref%chi_p_source(:,i) = chi_l_constants(i,3)*chi_l_functions(:,3,i) + endif + temp_functions(:,3,i) = chi_l_functions(:,3,i) + temp_constants(i,3) = chi_l_constants(i,3) + Endif + + If (active) Then + If (chi_l_use_custom_constant(4,i) .and. chi_l_use_custom_function(5,i)) Then + If (my_rank .eq. 0) Then + Call stdout%print('Chi buoyancy source for active scalar field '//Adjustl(sind)//' has been set to:') + Call stdout%print('d_4*f_5') + Call stdout%print(' ') + Endif + ref%chi_buoyancy_coeff(:,i) = chi_l_constants(i,4)*chi_l_functions(:,5,i) + temp_functions(:,5,i) = chi_l_functions(:,5,i) + temp_constants(i,4) = chi_l_constants(i,4) + Endif + Endif + + If (chi_l_use_custom_function(1,i)) Then + ! Set d_1 = 1 by default + If (.not. chi_l_use_custom_constant(1,i)) Then + If (my_rank .eq. 0) Then + Call stdout%print("User didn't set d_1 for "//Adjustl(chi_l_type_str)// & + ' scalar field '//Adjustl(sind)//'.') + Call stdout%print('Now setting d_1 to 1.') + Endif + chi_l_constants(i,1) = 1.0d0 + Else + If (my_rank .eq. 0) Then + Write(dstring,dofmt) chi_l_constants(i,1) + Call stdout%print('User set d_1 to '//Adjustl(dstring)//' for '// & + Adjustl(chi_l_type_str)//' scalar field '// & + Adjustl(sind)//'.') + Endif + Endif + If (my_rank .eq. 0) Then + Call stdout%print('Background dchidr for '//Adjustl(chi_l_type_str)// & + ' scalar field '//Adjustl(sind)//' has been set to:') + Call stdout%print('d_1*g_1') + Call stdout%print(' ') + Endif + ! can't declare items in type as a target so have to do this if statement instead + if (active) then + ref%dchirefadr(:,i) = chi_l_constants(i,1)*chi_l_functions(:,1,i) + else + ref%dchirefpdr(:,i) = chi_l_constants(i,1)*chi_l_functions(:,1,i) + endif + temp_functions(:,1,i) = chi_l_functions(:,1,i) + temp_constants(i,1) = chi_l_constants(i,1) + Endif + Enddo + + chi_l_functions(:,:,:) = temp_functions(:,:,:) + chi_l_constants(1:n_l_scalars, 1:n_chi_constants) = temp_constants(:,:) + DeAllocate(temp_functions, temp_constants) + + End Subroutine Augment_Scalar_Reference + Subroutine Constant_Energy_Heating() Implicit None ! rho T dSdt = alpha : energy deposition per unit volume is constant @@ -1262,6 +1482,20 @@ Subroutine Get_Custom_Reference() Call Read_Custom_Reference_File(custom_reference_file) + if (n_active_scalars > 0) then + If (my_rank .eq. 0) Then + Call stdout%print('And active scalars from: '//TRIM(ADJUSTL(chi_a_custom_reference_file)) ) + Endif + Call Read_Scalar_Custom_Reference_File(chi_a_custom_reference_file, .true.) + endif + + if (n_passive_scalars > 0) then + If (my_rank .eq. 0) Then + Call stdout%print('And passive scalars from: '//TRIM(ADJUSTL(chi_p_custom_reference_file)) ) + Endif + Call Read_Scalar_Custom_Reference_File(chi_p_custom_reference_file, .false.) + endif + Do i=1,4 fi = fi_to_check(i) If (ra_function_set(fi) .eq. 0) Then @@ -1277,7 +1511,7 @@ Subroutine Get_Custom_Reference() ref%d2lnrho(:) = ra_functions(:,9) ref%buoyancy_coeff(:) = ra_constants(2)*ra_functions(:,2) Do i = 1, n_active_scalars - ref%chi_buoyancy_coeff(i,:) = ra_constants(12+(i-1)*2)*ra_functions(:,2) + ref%chi_buoyancy_coeff(:,i) = chi_a_constants(i,4)*chi_a_functions(:,5,i) Enddo ref%temperature(:) = ra_functions(:,4) @@ -1292,6 +1526,12 @@ Subroutine Get_Custom_Reference() Endif ref%heating(:) = ra_functions(:,6)/(ref%density*ref%temperature)*ra_constants(10) + do i = 1, n_active_scalars + ref%chi_a_source(:,i) = chi_a_constants(i,3)*chi_a_functions(:,3,i) + enddo + do i = 1, n_passive_scalars + ref%chi_p_source(:,i) = chi_p_constants(i,3)*chi_p_functions(:,3,i) + enddo ref%Coriolis_Coeff = ra_constants(1) If (Angular_Velocity .gt. 0) Then @@ -1315,6 +1555,14 @@ Subroutine Get_Custom_Reference() ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat Call log_deriv(ref%dsdr_over_cp(:), ref%d2s_over_cp(:), no_log=.true.) + do i = 1, n_active_scalars + ref%dchirefadr(:,i) = chi_a_constants(i,1)*chi_a_functions(:,1,i) + enddo + + do i = 1, n_passive_scalars + ref%dchirefpdr(:,i) = chi_p_constants(i,1)*chi_p_functions(:,1,i) + enddo + End Subroutine Get_Custom_Reference Subroutine Write_Equation_Coefficients_File(filename) @@ -1356,7 +1604,7 @@ Subroutine Read_Custom_Reference_File(filename) Character*120, Intent(In) :: filename Character*120 :: ref_file Integer :: pi_integer,nr_ref, eqversion - Integer :: i, k, j, n_scalars, dummy + Integer :: i, k, j, dummy Integer :: cset(1:n_ra_constants), fset(1:n_ra_functions) Real*8 :: input_constants(1:n_ra_constants) Real*8, Allocatable :: ref_arr_old(:,:), rtmp(:), rtmp2(:) @@ -1396,14 +1644,13 @@ Subroutine Read_Custom_Reference_File(filename) Read(15) eqversion If (eqversion .eq. 1) Then !Read(15) cset(1:n_ra_constants-1) ! c_11 didn't exist yet - Read(15) cset(1:10) ! equation_coefficients couldn't write the custom active/passive scalar constants yet, and c_11 didn't exist yet + Read(15) cset(1:10) ! c_11 didn't exist yet Read(15) fset(1:n_ra_functions) Read(15) input_constants(1:10) cset(11) = 1 ! treat this as if c_11 = 1 was specified in the custom reference file input_constants(11) = 1.0d0 Else - Read(15) dummy ! n_ra_constants, but its up to the user to make sure this is the same - ! same number specified in Rayleigh via "n_active_scalars" + Read(15) dummy ! n_ra_constants, which Rayleigh should already know Read(15) dummy ! n_ra_functions, which Rayleigh should already know Read(15) cset(1:n_ra_constants) Read(15) fset(1:n_ra_functions) @@ -1529,12 +1776,6 @@ Subroutine Read_Custom_Reference_File(filename) If ((fset(13) .eq. 0) .and. (fset(7) .eq. 1)) Then Call log_deriv(ra_functions(:,7), ra_functions(:,13)) !dlneta Endif - n_scalars = n_active_scalars + n_passive_scalars - do i = 0, (n_scalars - 1) - If ((fset(16+i*2) .eq. 0) .and. (fset(15+i*2) .eq. 1)) Then - Call log_deriv(ra_functions(:,15+i*2), ra_functions(:,16+i*2)) !dlnkappa_chi - Endif - end do Else Call stdout%print('Error. This file appears to be corrupt (check Endian convention).') Endif @@ -1581,6 +1822,276 @@ Subroutine Read_Custom_Reference_File(filename) End Subroutine Read_Custom_Reference_File + Subroutine Read_Scalar_Custom_Reference_File(filename, active) + Character*120, Intent(In) :: filename + Logical, Intent(In) :: active ! if .true., read in active scalar info; otherwise, read in passive scalar info + Character*120 :: ref_file + Integer :: pi_integer,nr_ref, eqversion + Integer :: i, k, j, n_scalars, dummy, ios + Integer :: cset(1:n_chi_constants,1:n_scalar_max), fset(1:n_chi_functions,1:n_scalar_max) + Real*8 :: input_constants(1:n_chi_constants,1:n_scalar_max) + Real*8, Allocatable :: ref_arr_old(:,:,:), rtmp(:), rtmp2(:) + Real*8, Allocatable :: old_radius(:) + Character*12 :: dstring + Character*8 :: dofmt = '(ES12.5)' + Character(len=2) :: cind, sind + Character*3 :: intstr + ! variables to point at different global parameters depending on whether we are reading in active or passive scalar info + Integer :: n_l_scalars + Real*8, pointer :: chi_l_constants(:,:), chi_l_functions(:,:,:) + Logical, pointer :: override_chi_l_constant(:,:) + Logical, pointer :: chi_l_use_custom_constant(:,:), chi_l_use_custom_function(:,:) + Integer, pointer :: chi_l_constant_set(:,:), chi_l_function_set(:,:) + Integer, pointer :: with_custom_chi_l_constants(:,:), with_custom_chi_l_functions(:,:) + Character*7 :: chi_l_type_str + + + if (active) then + n_l_scalars = n_active_scalars + override_chi_l_constant => override_chi_a_constant + chi_l_constant_set => chi_a_constant_set + chi_l_function_set => chi_a_function_set + chi_l_constants => chi_a_constants + chi_l_functions => chi_a_functions + with_custom_chi_l_constants => with_custom_chi_a_constants + with_custom_chi_l_functions => with_custom_chi_a_functions + chi_l_use_custom_constant => chi_a_use_custom_constant + chi_l_use_custom_function => chi_a_use_custom_function + chi_l_type_str = 'active' + else + n_l_scalars = n_passive_scalars + override_chi_l_constant => override_chi_p_constant + chi_l_constant_set => chi_p_constant_set + chi_l_function_set => chi_p_function_set + chi_l_constants => chi_p_constants + chi_l_functions => chi_p_functions + with_custom_chi_l_constants => with_custom_chi_p_constants + with_custom_chi_l_functions => with_custom_chi_p_functions + chi_l_use_custom_constant => chi_p_use_custom_constant + chi_l_use_custom_function => chi_p_use_custom_function + chi_l_type_str = 'passive' + end if + + cset(:,:) = 0 + fset(:,:) = 0 + input_constants(:,:) = 0.0d0 + + ref_file = Trim(my_path)//filename + + Open(unit=15,file=ref_file,form='unformatted', status='old',access='stream',iostat=ios) + if (ios /= 0) then + Call stdout%print('WARNING: custom reference file for '//Adjustl(chi_l_type_str)//' scalars not found, skipping.') + return + endif + + !Verify Endianness + Read(15)pi_integer + If (pi_integer .ne. 314) Then + close(15) + Open(unit=15,file=ref_file,form='unformatted', status='old', & + CONVERT = 'BIG_ENDIAN' , access='stream') + Read(15)pi_integer + If (pi_integer .ne. 314) Then + Close(15) + Open(unit=15,file=ref_file,form='unformatted', status='old', & + CONVERT = 'LITTLE_ENDIAN' , access='stream') + Read(15)pi_integer + Endif + Endif + + If (pi_integer .eq. 314) Then + + ! Read in constants and their 'set' flags + Read(15) eqversion + Read(15) dummy ! n_chi_constants + Read(15) dummy ! n_chi_functions + Read(15) n_scalars + if (n_scalars .gt. n_l_scalars) then + Call stdout%print('WARNING: number of scalars in the custom reference file') + Call stdout%print('exceeds the number of '//Adjustl(chi_l_type_str)//' scalars') + Call stdout%print('specified in main_input.') + Call stdout%print('Rayleigh will ignore the extra reference states specified in the file') + Call stdout%print('and assume that the first reference states specified correspond to the') + Call stdout%print('fields described in main_input.') + Call stdout%print('Please check your scalar custom reference file and main_input settings.') + end if + Read(15) ((cset(i, j), i = 1, n_chi_constants), j = 1, n_scalars) + Read(15) ((fset(i, j), i = 1, n_chi_functions), j = 1, n_scalars) + Read(15) ((input_constants(i, j), i = 1, n_chi_constants), j = 1, n_scalars) + + ! cset(i,j) is 1 if a constant(i) was set for scalar j; it is 0 otherwise. + ! The logic below deals with a constant set in both the reference + ! file and in main_input. Main_input values take precedence if + ! override_chi_a_constant(s) is set or if the reference file constant + ! was not set. + Do j = 1, n_l_scalars + Do i = 1, n_chi_constants + If ( (.not. override_constants) .and. (.not. override_chi_l_constant(j,i)) ) Then + chi_l_constants(j,i) = chi_l_constants(j,i) + cset(i,j)*(input_constants(i,j)-chi_l_constants(j,i)) + Endif + Enddo + Enddo + + ! determine which functions/constants were set by the user + chi_l_function_set(:,1:n_l_scalars) = fset(:,1:n_l_scalars) + Do j = 1, n_l_scalars + Do i = 1, n_chi_constants + If ((cset(i,j) .eq. 1) .or. override_chi_l_constant(j,i) .or. override_constants) Then + chi_l_constant_set(i,j) = 1 + Endif + Enddo + Enddo + + ! Print the values of the constants + Do j = 1, n_l_scalars + Do i = 1, n_chi_constants + If (my_rank .eq. 0) Then + Write(cind, '(I2)') i + Write(sind,'(I2)') j + Write(dstring,dofmt) chi_l_constants(j,i) + Call stdout%print('d_'//Adjustl(cind)//' for '//Adjustl(chi_l_type_str)// & + ' scalar '//Adjustl(sind)//' = '//Trim(dstring)) + Endif + Enddo + Enddo + + ! Read the reference file's radial grid + Read(15) nr_ref + Allocate(ref_arr_old(1:nr_ref, 1:n_chi_functions, 1:n_scalars)) + Allocate(old_radius(1:nr_ref)) + + Read(15)(old_radius(i),i=1,nr_ref) + Do k = 1, n_scalars + Do j = 1, n_chi_functions + Read(15)(ref_arr_old(i, j, k),i=1,nr_ref) + EndDo + EndDo + + !Check to see if radius is tabulated in ascending or descending order. + !If it is found to be in ascending order, reverse the radius and the + !input array of functions + If (old_radius(1) .lt. old_radius(nr_ref)) Then + + If (my_rank .eq. 0) Call stdout%print('Reversing Radial Indices in Custom Ref File!') + + Allocate(rtmp(1:nr_ref)) + + rtmp = old_radius + Do i = 1, nr_ref + old_radius(i) = rtmp(nr_ref-i+1) + Enddo + + Do k = 1, n_scalars + Do j = 1, n_chi_functions + rtmp(:) = ref_arr_old(:,j,k) + Do i = 1, nr_ref + ref_arr_old(i,j,k) = rtmp(nr_ref-i+1) + Enddo + Enddo + Enddo + + DeAllocate(rtmp) + + Endif + + Close(15) + if (active) then + chi_a_custom_reference_read = .true. + else + chi_p_custom_reference_read = .true. + end if + + If (nr_ref .ne. n_r) Then + !Interpolate onto the current radial grid if necessary + !Note that the underlying assumption here is that same # of grid points + ! means same grid - come back to this later for generality + Allocate( rtmp2(1:n_r)) + Allocate( rtmp(1:nr_ref)) + + Do k = 1, n_l_scalars + Do j = 1, n_chi_functions + rtmp(:) = ref_arr_old(:,j,k) + rtmp2(:) = 0.0d0 + Call Spline_Interpolate(rtmp, old_radius, rtmp2, radius) + + chi_l_functions(1:n_r,j,k) = rtmp2 + Enddo + Enddo + + DeAllocate(rtmp,rtmp2) + Else + + ! Bit redundant here, but may want to do filtering on ref_arr array + chi_l_functions(1:n_r, 1:n_chi_functions, 1:n_l_scalars) = & + ref_arr_old(1:n_r,1:n_chi_functions,1:n_l_scalars) + + If (my_rank .eq. 0) Then + call stdout%print(" WARNING: nr = nr_old. Assuming grids are the same.") + Endif + Endif + DeAllocate(ref_arr_old,old_radius) + + ! Finally, if the logarithmic derivatives of kappa_chi was + ! not specified, then we compute it here. + ! only calculate the log derivative if kappa_chi was set, otherwise there + ! are divide by zero issues + Do k = 1, n_l_scalars + If ((fset(4,k) .eq. 0) .and. (fset(2,k) .eq. 1)) Then + Call log_deriv(chi_l_functions(:,2,k), chi_l_functions(:,4,k)) !dlnkappa_chi + Endif + Enddo + Else + Call stdout%print('Error. This file appears to be corrupt (check Endian convention).') + Endif + + ! only used if user wants to change reference_type=1,2,3 + If (with_custom_reference) Then + Do i = 1, n_l_scalars + Do k = 1, n_chi_constants + j = with_custom_chi_l_constants(i,k) + if ((j .gt. 0) .and. (j .le. n_chi_constants)) Then + If (chi_l_constant_set(j,i) .eq. 1) Then + chi_l_use_custom_constant(j,i) = .true. + Else + If (my_rank .eq. 0) Then + Write(intstr,'(i3)') j + Write(sind, '(I2)') i + Call stdout%print(' ') + Call stdout%print('You set with_custom_constant: '//TRIM(ADJUSTL(intstr))) + Call stdout%print('for '//Adjustl(chi_l_type_str)//' scalar '//Adjustl(sind)) + Call stdout%print('But this constant was not set in either main_input or ') + Call stdout%print('the custom reference file. Selection will be ignored.') + Call stdout%print(' ') + Endif + Endif + Endif + Enddo + Enddo + + Do i = 1, n_l_scalars + Do k = 1, n_chi_functions + j = with_custom_chi_l_functions(i,k) + if ((j .gt. 0) .and. (j .le. n_chi_functions)) Then + If (chi_l_function_set(j,i) .eq. 1) Then + chi_l_use_custom_function(j,i) = .true. + Else + If (my_rank .eq. 0) Then + Write(intstr,'(i3)') j + Write(sind, '(I2)') i + Call stdout%print(' ') + Call stdout%print('You set with_custom_function: '//TRIM(ADJUSTL(intstr))) + Call stdout%print('for '//Adjustl(chi_l_type_str)//' scalar '//Adjustl(sind)) + Call stdout%print('But this function was not set. Selection will be ignored.') + Call stdout%print(' ') + Endif + Endif + Endif + Enddo + Enddo + Endif + + End Subroutine Read_Scalar_Custom_Reference_File + Subroutine Log_Deriv(arr1,arr2, no_log) Implicit None Real*8, Intent(In) :: arr1(:) @@ -1663,7 +2174,8 @@ Subroutine Restore_Reference_Defaults Angular_Velocity = -1.0d0 ! Frame rotation rate (sets Coriolis force) ! Custom reference-state variables (reference_type = 4) - ! NOTE: n_ra_constants / n_ra_functions do not have default values (but maybe do, if you consider the Allocate_Reference_State() routine) + ! NOTE: n_ra_constants / n_ra_functions do not have default values + ! (but maybe do, if you consider the Allocate_Reference_State() routine) custom_reference_file ='nothing' ! Internal heating variables @@ -1719,17 +2231,22 @@ Subroutine Restore_Reference_Defaults If (allocated(ref%Temperature)) DeAllocate(ref%Temperature) If (allocated(ref%dlnT)) DeAllocate(ref%dlnT) - If (allocated(ref%entropy)) DeAllocate(ref%entropy) - If (allocated(ref%exp_entropy)) DeAllocate(ref%exp_entropy) + If (allocated(ref%entropy)) DeAllocate(ref%entropy) + If (allocated(ref%exp_entropy)) DeAllocate(ref%exp_entropy) If (allocated(ref%dsdr)) DeAllocate(ref%dsdr) If (allocated(ref%dsdr_over_cp)) DeAllocate(ref%dsdr_over_cp) If (allocated(ref%d2s_over_cp)) DeAllocate(ref%d2s_over_cp) If (allocated(ref%heating)) DeAllocate(ref%heating) - ! NOTE: ref%Coriolis_Coeff and ref%Lorentz_Coeff have no default values (but maybe do, if you consider the Allocate_Reference_State() routine) + ! NOTE: ref%Coriolis_Coeff and ref%Lorentz_Coeff have no default values + ! (but maybe do, if you consider the Allocate_Reference_State() routine) If (allocated(ref%Buoyancy_Coeff)) DeAllocate(ref%Buoyancy_Coeff) If (allocated(ref%chi_buoyancy_coeff)) DeAllocate(ref%chi_buoyancy_coeff) + If (allocated(ref%chi_a_source)) DeAllocate(ref%chi_a_source) + If (allocated(ref%chi_p_source)) DeAllocate(ref%chi_p_source) + If (allocated(ref%dchirefadr)) DeAllocate(ref%dchirefadr) + If (allocated(ref%dchirefpdr)) DeAllocate(ref%dchirefpdr) If (allocated(ref%dpdr_w_term)) DeAllocate(ref%dpdr_w_term) If (allocated(ref%pressure_dwdr_term)) DeAllocate(ref%pressure_dwdr_term) @@ -1776,9 +2293,13 @@ Subroutine Initialize_Transport_Coefficients() Implicit None Integer :: i Real*8, Allocatable :: temp_functions(:,:), temp_constants(:) - Logical :: restore, need_custom + Real*8, Allocatable :: temp_chi_a_functions(:,:,:), temp_chi_a_constants(:,:) + Real*8, Allocatable :: temp_chi_p_functions(:,:,:), temp_chi_p_constants(:,:) + Logical :: restore, restore_chi_a, restore_chi_p, need_custom restore = .false. + restore_chi_a = .false. + restore_chi_p = .false. Call Allocate_Transport_Coefficients @@ -1789,45 +2310,79 @@ Subroutine Initialize_Transport_Coefficients() need_custom = .true. EndIf + If ((.not. custom_reference_read) .and. need_custom) Then + Allocate(temp_functions(1:n_r, 1:n_ra_functions)) + Allocate(temp_constants(1:n_ra_constants)) + temp_functions(:,:) = ra_functions(:,:) + ! Note that ra_constants is allocated up to max_ra_constants, + ! which could be more than n_ra_constants + temp_constants(:) = ra_constants(1:n_ra_constants) + restore = .true. + ! If we read the custom file, we may overwrite things besides the diffusion coefficients + ! We "back up" the current reference state in temp_constants and temp_functions + ! Below, we modify only the "temp" equation coefficients associated with custom diffusions + ! Then we restore ra_constants and ra_functions from the "temp" arrays + Call Read_Custom_Reference_File(custom_reference_file) + EndIf + + ! reset to false + need_custom = .false. Do i = 1, n_active_scalars If (kappa_chi_a_type(i) .eq. 3) Then need_custom = .true. Endif Enddo + If ((.not. chi_a_custom_reference_read) .and. need_custom) Then + Allocate(temp_chi_a_functions(1:n_r, 1:n_chi_functions, 1:n_active_scalars)) + Allocate(temp_chi_a_constants(1:n_active_scalars, 1:n_chi_constants)) + temp_chi_a_functions(:,:,:) = chi_a_functions(:,:,:) + ! Note that chi_a_constants is allocated up to n_scalar_max x max_chi_constants, + ! which could be more than n_active_scalars xn_chi_constants + temp_chi_a_constants(:,:) = chi_a_constants(1:n_active_scalars, 1:n_chi_constants) + restore_chi_a = .true. + ! If we read the custom file, we may overwrite things besides the diffusion coefficients + ! We "back up" the current reference state in temp_constants and temp_functions + ! Below, we modify only the "temp" equation coefficients associated with custom diffusions + ! Then we restore ra_constants and ra_functions from the "temp" arrays + Call Read_Scalar_Custom_Reference_File(chi_a_custom_reference_file, .true.) + EndIf + + ! reset to false + need_custom = .false. Do i = 1, n_passive_scalars If (kappa_chi_p_type(i) .eq. 3) Then need_custom = .true. Endif Enddo - If ((.not. custom_reference_read) .and. need_custom) Then - Allocate(temp_functions(1:n_r, 1:n_ra_functions)) - Allocate(temp_constants(1:n_ra_constants)) - temp_functions(:,:) = ra_functions(:,:) - ! Note that ra_constants is allocated up to max_ra_constants, - ! which could be more than n_ra_constants - temp_constants(:) = ra_constants(1:n_ra_constants) - restore = .true. + If ((.not. chi_p_custom_reference_read) .and. need_custom) Then + Allocate(temp_chi_p_functions(1:n_r, 1:n_chi_functions, 1:n_passive_scalars)) + Allocate(temp_chi_p_constants(1:n_passive_scalars, 1:n_chi_constants)) + temp_chi_p_functions(:,:,:) = chi_p_functions(:,:,:) + ! Note that chi_p_constants is allocated up to n_scalar_max x max_chi_constants, + ! which could be more than n_passive_scalars x n_chi_constants + temp_chi_p_constants(:,:) = chi_p_constants(1:n_passive_scalars, 1:n_chi_constants) + restore_chi_p = .true. ! If we read the custom file, we may overwrite things besides the diffusion coefficients ! We "back up" the current reference state in temp_constants and temp_functions ! Below, we modify only the "temp" equation coefficients associated with custom diffusions ! Then we restore ra_constants and ra_functions from the "temp" arrays - Call Read_Custom_Reference_File(custom_reference_file) + Call Read_Scalar_Custom_Reference_File(chi_p_custom_reference_file, .false.) EndIf Call Initialize_Diffusivity(nu,dlnu,nu_top,nu_type,nu_power,5,3,11) Call Initialize_Diffusivity(kappa,dlnkappa,kappa_top,kappa_type,kappa_power,6,5,12) do i = 1, n_active_scalars - Call Initialize_Diffusivity(kappa_chi_a(i,:),dlnkappa_chi_a(i,:),& + Call Initialize_Scalar_Diffusivity(kappa_chi_a(i,:),dlnkappa_chi_a(i,:),& kappa_chi_a_top(i),kappa_chi_a_type(i),kappa_chi_a_power(i),& - 12+(i-1)*2,15+(i-1)*2,16+(i-1)*2) + i,2,2,4,.true.) end do do i = 1, n_passive_scalars - Call Initialize_Diffusivity(kappa_chi_p(i,:),dlnkappa_chi_p(i,:),& + Call Initialize_Scalar_Diffusivity(kappa_chi_p(i,:),dlnkappa_chi_p(i,:),& kappa_chi_p_top(i),kappa_chi_p_type(i),kappa_chi_p_power(i),& - 12+(n_active_scalars+i-1)*2,15+(n_active_scalars+i-1)*2,16+(n_active_scalars+i-1)*2) + i,2,2,4,.false.) end do If (viscous_heating) Then @@ -1861,22 +2416,6 @@ Subroutine Initialize_Transport_Coefficients() temp_constants(6) = ra_constants(6) Endif - do i = 0, n_active_scalars-1 - If (kappa_chi_a_type(i+1) .eq. 3) Then - temp_functions(:,15+i*2) = ra_functions(:,15+i*2) - temp_functions(:,16+i*2) = ra_functions(:,16+i*2) - temp_constants(12+i*2) = ra_constants(12+i*2) - Endif - end do - - do i = 0, n_passive_scalars-1 - If (kappa_chi_p_type(i+1) .eq. 3) Then - temp_functions(:,15+(n_active_scalars+i)*2) = ra_functions(:,15+(n_active_scalars+i)*2) - temp_functions(:,16+(n_active_scalars+i)*2) = ra_functions(:,16+(n_active_scalars+i)*2) - temp_constants(12+(n_active_scalars+i)*2) = ra_constants(12+(n_active_scalars+i)*2) - Endif - end do - If (nu_type .eq. 3) Then temp_functions(:,3) = ra_functions(:,3) temp_functions(:,11) = ra_functions(:,11) @@ -1886,9 +2425,35 @@ Subroutine Initialize_Transport_Coefficients() ra_constants(1:n_ra_constants) = temp_constants(:) ra_functions(:,:) = temp_functions(:,:) DeAllocate(temp_functions, temp_constants) + endif + if (restore_chi_a) then - Endif + do i = 1, n_active_scalars + If (kappa_chi_a_type(i) .eq. 3) Then + temp_chi_a_functions(:,2,i) = chi_a_functions(:,2,i) + temp_chi_a_functions(:,4,i) = chi_a_functions(:,4,i) + temp_chi_a_constants(i,2) = chi_a_constants(i,2) + Endif + end do + + chi_a_functions(:,:,:) = temp_chi_a_functions(:,:,:) + chi_a_constants(1:n_active_scalars, 1:n_chi_constants) = temp_chi_a_constants(:,:) + endif + + if (restore_chi_p) then + + do i = 1, n_passive_scalars + If (kappa_chi_p_type(i) .eq. 3) Then + temp_chi_p_functions(:,2,i) = chi_p_functions(:,2,i) + temp_chi_p_functions(:,4,i) = chi_p_functions(:,4,i) + temp_chi_p_constants(i,2) = chi_p_constants(i,2) + Endif + end do + + chi_p_functions(:,:,:) = temp_chi_p_functions(:,:,:) + chi_p_constants(1:n_passive_scalars, 1:n_chi_constants) = temp_chi_p_constants(:,:) + endif Call Compute_Diffusion_Coefs() @@ -1976,6 +2541,96 @@ Subroutine Initialize_Diffusivity(x,dlnx,xtop,xtype,xpower,ci,fi,dlnfi) End Subroutine Initialize_Diffusivity + Subroutine Initialize_Scalar_Diffusivity(x,dlnx,xtop,xtype,xpower,si,ci,fi,dlnfi,active) + Implicit None + Real*8, Intent(InOut) :: x(:), dlnx(:) + Real*8, Intent(InOut) :: xtop + Integer, Intent(In) :: si, ci, fi, dlnfi, xtype + Logical, Intent(In) :: active + Real*8, Intent(In) :: xpower + Character(len=2) :: ind, sind + Real*8 :: norm + ! variables to point at different global parameters depending on whether we are reading in active or passive scalar info + Real*8, pointer :: chi_l_constants(:,:), chi_l_functions(:,:,:) + Integer, pointer :: chi_l_constant_set(:,:), chi_l_function_set(:,:) + Character*7 :: chi_l_type_str + + if (active) then + chi_l_constant_set => chi_a_constant_set + chi_l_function_set => chi_a_function_set + chi_l_constants => chi_a_constants + chi_l_functions => chi_a_functions + chi_l_type_str = 'active' + else + chi_l_constant_set => chi_p_constant_set + chi_l_function_set => chi_p_function_set + chi_l_constants => chi_p_constants + chi_l_functions => chi_p_functions + chi_l_type_str = 'passive' + end if + + If (reference_type .eq. 4) Then + If (xtop .le. 0) Then + If (chi_l_constant_set(ci,si) .eq. 0) Then + If (my_rank .eq. 0) Then + Write(ind, '(I2)') ci + Write(sind, '(I2)') si + Call stdout%print('ERROR: constant c_'//Trim(Adjustl(ind))//' for ') + Call stdout%print(Adjustl(chi_l_type_str)//' scalar field '//Adjustl(sind)) + Call stdout%print('must be set in the custom reference file') + Endif + Else + xtop = chi_l_constants(si,ci) + Endif + Endif + Endif + + Select Case(xtype) + Case(1) + x(:) = xtop + dlnx(:) = 0.0d0 + Case(2) + Call vary_with_density(x,dlnx,xtop, xpower) + Case(3) + If ((chi_l_function_set(fi,si) .eq. 1) .and. (chi_l_constant_set(ci,si) .eq. 1)) Then + + x(:) = chi_l_constants(si,ci)*chi_l_functions(:,fi,si) + dlnx(:) = chi_l_functions(:,dlnfi,si) + xtop = x(1) + ! Nothing to be done here for functions and constants -- completely set + ElseIf ((chi_l_function_set(fi,si) .eq. 1) .and. (chi_l_constant_set(ci,si) .eq. 0)) Then + x(:) = xtop*chi_l_functions(:,fi,si) + dlnx(:) = chi_l_functions(:,dlnfi,si) + Else + If (my_rank .eq. 0) Then + Write(ind, '(I2)') fi + Write(sind, '(I2)') si + Call stdout%print('ERROR: function f_'//Trim(Adjustl(ind))//' for') + Call stdout%print(Adjustl(chi_l_type_str)//' scalar field '//Adjustl(sind)) + Call stdout%print('must be set in the custom reference file') + EndIf + EndIf + + End Select + + ! We potentially need to renormalize for the Polytropic_ReferenceND_General() non-dimensionalization + If (reference_type .eq. 5) Then + If (ND_Inner_Radius) Then + ! We want "xtop" to be the value at the inner radius instead + norm = x(N_R) + x = xtop * (x/norm) + xtop = x(1) + Elseif (ND_Volume_Average) Then + ! We want "xtop" to be the value averaged over the whole shell instead + Call Integrate_in_radius(x,norm) + norm = four_pi*norm/shell_volume + x = xtop * (x/norm) + xtop = x(1) + Endif + Endif + + End Subroutine Initialize_Scalar_Diffusivity + Subroutine Vary_With_Density(coeff, dln, coeff_top, coeff_power) Implicit None Real*8, Intent(InOut) :: coeff(:), dln(:) @@ -2218,17 +2873,17 @@ Subroutine Set_Diffusivity_Equation_Coefficients Do i = 1, n_active_scalars If (kappa_chi_a_type(i) .ne. 3) Then - ra_constants(12+(i-1)*2) = kappa_chi_a_norm(i) - ra_functions(:,15+(i-1)*2) = kappa_chi_a(i,:)/kappa_chi_a_norm(i) - ra_functions(:,16+(i-1)*2) = dlnkappa_chi_a(i,:) + chi_a_constants(i,2) = kappa_chi_a_norm(i) + chi_a_functions(:,2,i) = kappa_chi_a(i,:)/kappa_chi_a_norm(i) + chi_a_functions(:,4,i) = dlnkappa_chi_a(i,:) Endif Enddo Do i = 1, n_passive_scalars If (kappa_chi_p_type(i) .ne. 3) Then - ra_constants(12+(n_active_scalars+i-1)*2) = kappa_chi_p_norm(i) - ra_functions(:,15+(n_active_scalars+i-1)*2) = kappa_chi_p(i,:)/kappa_chi_p_norm(i) - ra_functions(:,16+(n_active_scalars+i-1)*2) = dlnkappa_chi_p(i,:) + chi_p_constants(i,2) = kappa_chi_p_norm(i) + chi_p_functions(:,2,i) = kappa_chi_p(i,:)/kappa_chi_p_norm(i) + chi_p_functions(:,4,i) = dlnkappa_chi_p(i,:) Endif Enddo diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90 index adc4fb9b..5744bc1d 100755 --- a/src/Physics/Sphere_Linear_Terms.F90 +++ b/src/Physics/Sphere_Linear_Terms.F90 @@ -144,10 +144,16 @@ Subroutine Initialize_Linear_System do i = 1, n_active_scalars Call Initialize_Equation_Coefficients(weq,chiavar(i), 0,lp) Call Initialize_Equation_Coefficients(chiaeq(i),chiavar(i), 2,lp) + if (chi_a_advect_reference_state(i)) then + Call Initialize_Equation_Coefficients(chiaeq(i),wvar, 0,lp) + end if end do do i = 1, n_passive_scalars Call Initialize_Equation_Coefficients(chipeq(i),chipvar(i), 2,lp) + if (chi_p_advect_reference_state(i)) then + Call Initialize_Equation_Coefficients(chipeq(i),wvar, 0,lp) + end if end do If (advect_reference_state) Then @@ -239,7 +245,7 @@ Subroutine Load_Linear_Coefficients() Call add_implicit_term(peq, tvar, 0, amp,lp, static = .true.) ! Gravity --- Need LHS_Only Flag do i = 1, n_active_scalars - amp = -ref%chi_buoyancy_coeff(i,:) + amp = -ref%chi_buoyancy_coeff(:,i) Call add_implicit_term(peq, chiavar(i), 0, amp,lp, static = .true.) ! Gravity --- Need LHS_Only Flag end do @@ -298,7 +304,7 @@ Subroutine Load_Linear_Coefficients() ! Chi do i = 1, n_active_scalars - amp = -ref%chi_buoyancy_coeff(i,:)/H_Laplacian + amp = -ref%chi_buoyancy_coeff(:,i)/H_Laplacian If (pseudo_incompressible) Then amp = amp*ref%exp_entropy Endif @@ -413,7 +419,6 @@ Subroutine Load_Linear_Coefficients() Call add_implicit_term(teq,wvar,0,amp,lp) Endif - !CHIVAR equation ---- (no ell =0) do i = 1, n_active_scalars amp = 1.0d0 @@ -430,6 +435,11 @@ Subroutine Load_Linear_Coefficients() amp = H_Laplacian*kappa_chi_a(i,:)*diff_factor Call add_implicit_term(chiaeq(i),chiavar(i), 0, amp,lp) + + if (chi_a_advect_reference_state(i)) then + amp = H_Laplacian*ref%dchirefadr(:,i) + Call add_implicit_term(chiaeq(i),wvar, 0, amp,lp) + endif end do do i = 1, n_passive_scalars @@ -447,6 +457,11 @@ Subroutine Load_Linear_Coefficients() amp = H_Laplacian*kappa_chi_p(i,:)*diff_factor Call add_implicit_term(chipeq(i),chipvar(i), 0, amp,lp) + + if (chi_p_advect_reference_state(i)) then + amp = H_Laplacian*ref%dchirefpdr(:,i) + Call add_implicit_term(chipeq(i),wvar, 0, amp,lp) + endif end do !===================================================== diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index e1a85588..5921c213 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -183,11 +183,11 @@ Subroutine physical_space() do i = 1, n_active_scalars Call chi_Advection(chiavar(i), dchiadr(i), dchiadt(i), dchiadp(i)) - Call chi_Source_function(chiavar(i)) + Call chi_Source_function(chiavar(i), ref%chi_a_source(:,i)) end do do i = 1, n_passive_scalars Call chi_Advection(chipvar(i), dchipdr(i), dchipdt(i), dchipdp(i)) - Call chi_Source_function(chipvar(i)) + Call chi_Source_function(chipvar(i), ref%chi_p_source(:,i)) end do If (viscous_heating) Call Compute_Viscous_Heating() @@ -317,16 +317,17 @@ Subroutine Volumetric_Heating() End Subroutine Volumetric_Heating - Subroutine chi_Source_Function(chivar) + Subroutine chi_Source_Function(chivar, source) Implicit None - Integer :: chivar + Integer, intent(in) :: chivar + Real*8, intent(in) :: source(:) Integer :: t,r,k !$OMP PARALLEL DO PRIVATE(t,r,k) Do t = my_theta%min, my_theta%max Do r = my_r%min, my_r%max Do k =1, n_phi - wsp%p3b(k,r,t,chivar) = wsp%p3b(k,r,t,chivar)+0.0d0 ! This is where you would put a source function + wsp%p3b(k,r,t,chivar) = wsp%p3b(k,r,t,chivar)+source(r) Enddo Enddo Enddo diff --git a/src/Physics/Sphere_Spectral_Space.F90 b/src/Physics/Sphere_Spectral_Space.F90 index 53c997ab..413b2423 100755 --- a/src/Physics/Sphere_Spectral_Space.F90 +++ b/src/Physics/Sphere_Spectral_Space.F90 @@ -351,11 +351,17 @@ Subroutine Post_Solve() Call Add_Derivative(chiaeq(i),chiavar(i),2,wsp%p1b,wsp%p1a,d2chiadr2(i)) Call Add_Derivative(weq,chiavar(i),0, wsp%p1b,wsp%p1a,chiavar(i)) ! gravity + If (chi_a_advect_reference_state(i)) Then + Call Add_Derivative(chiaeq(i),wvar,0,wsp%p1b,wsp%p1a,wvar) + Endif end do do i = 1, n_passive_scalars Call Add_Derivative(chipeq(i),chipvar(i),0,wsp%p1b,wsp%p1a,chipvar(i)) Call Add_Derivative(chipeq(i),chipvar(i),1,wsp%p1b,wsp%p1a,dchipdr(i)) Call Add_Derivative(chipeq(i),chipvar(i),2,wsp%p1b,wsp%p1a,d2chipdr2(i)) + If (chi_p_advect_reference_state(i)) Then + Call Add_Derivative(chipeq(i),wvar,0,wsp%p1b,wsp%p1a,wvar) + Endif end do !/////////////////////////////// diff --git a/tests/chi_custom_reference/T.full/gen_ref_state.py b/tests/chi_custom_reference/T.full/gen_ref_state.py new file mode 100644 index 00000000..315cc5cf --- /dev/null +++ b/tests/chi_custom_reference/T.full/gen_ref_state.py @@ -0,0 +1,36 @@ +import reference_tools as rt +import numpy +# Grid Parameters +nr = 192 # Number of radial points - this is intentionally not the same as main_input as I want to test the interpolation of the reference state + +#aspect ratio +ar = 0.35e0 +# shell depth +d=1.e0 + +#outer radial boundary +ro=d/(1-ar) + +#inner radial boundary +ri=ar*ro + +#radial grid +radius=numpy.linspace(ri,ro,nr) + +ones = numpy.ones(nr,dtype='float64') +Q0 = 10.0 + +# Initialize an equation_coefficients structure consistent with the case 0 benchmark test +my_ref = rt.equation_coefficients(radius) + +my_ref.set_function(ones,'heating') # heating function +my_ref.set_constant(Q0,'luminosity') # multiplies the Heating function +my_ref.write('ref_state.dat') # Write the data file + +chi_ref = rt.scalar_equation_coefficients(1, radius) + +chi_ref.set_function(ones,'source_chi') +chi_ref.set_constant(Q0, 'source_chi_scale') +chi_ref.write('chi_ref_state.dat') + + diff --git a/tests/chi_custom_reference/T.full/main_input b/tests/chi_custom_reference/T.full/main_input new file mode 100644 index 00000000..14926658 --- /dev/null +++ b/tests/chi_custom_reference/T.full/main_input @@ -0,0 +1,77 @@ +&problemsize_namelist + n_r = 48 + n_theta = 64 + nprow = 2 + npcol = 2 + aspect_ratio = 0.35d0 + shell_depth = 1.0d0 +/ +&numerical_controls_namelist +/ +&physical_controls_namelist + benchmark_mode = 0 + rotation = .True. + magnetism = .false. + viscous_heating = .false. + ohmic_heating = .false. + n_active_scalars = 1 + n_passive_scalars = 0 +/ +&temporal_controls_namelist + max_time_step = 7.5e-5 + max_iterations = 17 + checkpoint_interval = 10000 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist +globalavg_values = 401,501,503,505,507,509,511,513,515,517,519,521,523,10001,10002,10003,10004,10005,10006,10007,10008,10009,10010,10011,10012 +globalavg_frequency = 17 +globalavg_nrec = 1 +!full3d_values = 1,2,3,501,507,10001 +!full3d_frequency = 5 +/ +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +fix_tvar_top = .false. +fix_tvar_bottom = .false. +fix_dtdr_top = .true. +fix_dtdr_bottom = .true. +dtdr_top_file = '../bench_dtdr_rmax_init' +dtdr_bottom_file = '../bench_dtdr_rmin_init' +fix_chivar_a_top(1) = .false. +fix_chivar_a_bottom(1) = .false. +fix_dchidr_a_top(1) = .true. +fix_dchidr_a_bottom(1) = .true. +dchidr_a_top_file(1) = '../bench_dtdr_rmax_init' +dchidr_a_bottom_file(1) = '../bench_dtdr_rmin_init' +/ +&Initial_Conditions_Namelist +init_type = 8 +t_init_file = '../bench_t_init' +chi_a_init_file(1) = '../bench_t_init' +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 1.0d5 +chi_a_rayleigh_number(1) = 0.0 +Prandtl_Number = 1.0d0 +chi_a_prandtl_number(1) = 1.0d0 +reference_type = 1 +heating_type = 2 ! overwritten by reference state +gravity_power = 1.0d0 ! g ~ radius +custom_reference_file = 'ref_state.dat' +with_custom_reference = .true. +with_custom_constants = 10 +with_custom_functions = 6 +chi_a_custom_reference_file = 'chi_ref_state.dat' +with_custom_chi_a_constants(1,1) = 3 +with_custom_chi_a_functions(1,1) = 3 +/ +&Transport_Namelist +/ diff --git a/tests/chi_custom_reference/T.split/gen_ref_state.py b/tests/chi_custom_reference/T.split/gen_ref_state.py new file mode 100644 index 00000000..800cd6f9 --- /dev/null +++ b/tests/chi_custom_reference/T.split/gen_ref_state.py @@ -0,0 +1,45 @@ +import reference_tools as rt +import numpy +# Grid Parameters +nr = 192 # Number of radial points - this is intentionally not the same as main_input as I want to test the interpolation of the reference state + +#aspect ratio +ar = 0.35e0 +# shell depth +d=1.e0 + +#outer radial boundary +ro=d/(1-ar) + +#inner radial boundary +ri=ar*ro + +#radial grid +radius=numpy.linspace(ri,ro,nr) + +ones = numpy.ones(nr,dtype='float64') +Q0 = 10.0 +dTrefdr = -ro*ri*(radius**-2) + +# Initialize an equation_coefficients structure consistent with the case 0 benchmark test +my_ref = rt.equation_coefficients(radius) + +my_ref.set_function(ones,'heating') # heating function +my_ref.set_constant(Q0,'luminosity') # multiplies the Heating function + +my_ref.set_function(dTrefdr,'ds_dr') +my_ref.set_constant(1.0, 'dsdr_scale') + +my_ref.write('ref_state.dat') # Write the data file + +chi_ref = rt.scalar_equation_coefficients(1, radius) + +chi_ref.set_function(ones,'source_chi') +chi_ref.set_constant(Q0, 'source_chi_scale') + +chi_ref.set_function(dTrefdr, 'dchirefdr') +chi_ref.set_constant(1.0, 'dchirefdr_scale') + +chi_ref.write('chi_ref_state.dat') + + diff --git a/tests/chi_custom_reference/T.split/main_input b/tests/chi_custom_reference/T.split/main_input new file mode 100644 index 00000000..94128c18 --- /dev/null +++ b/tests/chi_custom_reference/T.split/main_input @@ -0,0 +1,79 @@ +&problemsize_namelist + n_r = 48 + n_theta = 64 + nprow = 2 + npcol = 2 + aspect_ratio = 0.35d0 + shell_depth = 1.0d0 +/ +&numerical_controls_namelist +/ +&physical_controls_namelist + benchmark_mode = 0 + rotation = .True. + magnetism = .false. + viscous_heating = .false. + ohmic_heating = .false. + n_active_scalars = 1 + n_passive_scalars = 0 +/ +&temporal_controls_namelist + max_time_step = 7.5e-5 + max_iterations = 17 + checkpoint_interval = 10000 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist +globalavg_values = 401,501,503,505,507,509,511,513,515,517,519,521,523,10001,10002,10003,10004,10005,10006,10007,10008,10009,10010,10011,10012 +globalavg_frequency = 17 +globalavg_nrec = 1 +!full3d_values = 1,2,3,501,507,10001 +!full3d_frequency = 10 +/ +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +fix_tvar_top = .false. +fix_tvar_bottom = .false. +fix_dtdr_top = .true. +fix_dtdr_bottom = .true. +dtdr_top_file = '../bench_dtpdr_rmax_init' +dtdr_bottom_file = '../bench_dtpdr_rmin_init' +fix_chivar_a_top(1) = .false. +fix_chivar_a_bottom(1) = .false. +fix_dchidr_a_top(1) = .true. +fix_dchidr_a_bottom(1) = .true. +dchidr_a_top_file(1) = '../bench_dtpdr_rmax_init' +dchidr_a_bottom_file(1) = '../bench_dtpdr_rmin_init' +/ +&Initial_Conditions_Namelist +init_type = 8 +t_init_file = '../bench_tp_init' +chi_a_init_file(1) = '../bench_tp_init' +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 1.0d5 +chi_a_rayleigh_number(1) = 0.0 +Prandtl_Number = 1.0d0 +chi_a_prandtl_number(1) = 1.0d0 +reference_type = 1 +heating_type = 2 ! overwritten by reference state +gravity_power = 1.0d0 ! g ~ radius +custom_reference_file = 'ref_state.dat' +with_custom_reference = .true. +with_custom_constants = 10, 11 +with_custom_functions = 6, 14 +chi_a_custom_reference_file = 'chi_ref_state.dat' +with_custom_chi_a_constants(1,1) = 1 +with_custom_chi_a_constants(1,2) = 3 +with_custom_chi_a_functions(1,1) = 1 +with_custom_chi_a_functions(1,2) = 3 +/ +&Transport_Namelist +/ diff --git a/tests/chi_custom_reference/chi.full/gen_ref_state.py b/tests/chi_custom_reference/chi.full/gen_ref_state.py new file mode 100644 index 00000000..315cc5cf --- /dev/null +++ b/tests/chi_custom_reference/chi.full/gen_ref_state.py @@ -0,0 +1,36 @@ +import reference_tools as rt +import numpy +# Grid Parameters +nr = 192 # Number of radial points - this is intentionally not the same as main_input as I want to test the interpolation of the reference state + +#aspect ratio +ar = 0.35e0 +# shell depth +d=1.e0 + +#outer radial boundary +ro=d/(1-ar) + +#inner radial boundary +ri=ar*ro + +#radial grid +radius=numpy.linspace(ri,ro,nr) + +ones = numpy.ones(nr,dtype='float64') +Q0 = 10.0 + +# Initialize an equation_coefficients structure consistent with the case 0 benchmark test +my_ref = rt.equation_coefficients(radius) + +my_ref.set_function(ones,'heating') # heating function +my_ref.set_constant(Q0,'luminosity') # multiplies the Heating function +my_ref.write('ref_state.dat') # Write the data file + +chi_ref = rt.scalar_equation_coefficients(1, radius) + +chi_ref.set_function(ones,'source_chi') +chi_ref.set_constant(Q0, 'source_chi_scale') +chi_ref.write('chi_ref_state.dat') + + diff --git a/tests/chi_custom_reference/chi.full/main_input b/tests/chi_custom_reference/chi.full/main_input new file mode 100644 index 00000000..414d5087 --- /dev/null +++ b/tests/chi_custom_reference/chi.full/main_input @@ -0,0 +1,77 @@ +&problemsize_namelist + n_r = 48 + n_theta = 64 + nprow = 2 + npcol = 2 + aspect_ratio = 0.35d0 + shell_depth = 1.0d0 +/ +&numerical_controls_namelist +/ +&physical_controls_namelist + benchmark_mode = 0 + rotation = .True. + magnetism = .false. + viscous_heating = .false. + ohmic_heating = .false. + n_active_scalars = 1 + n_passive_scalars = 0 +/ +&temporal_controls_namelist + max_time_step = 7.5e-5 + max_iterations = 17 + checkpoint_interval = 10000 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist +globalavg_values = 401,501,503,505,507,509,511,513,515,517,519,521,523,10001,10002,10003,10004,10005,10006,10007,10008,10009,10010,10011,10012 +globalavg_frequency = 17 +globalavg_nrec = 1 +!full3d_values = 1,2,3,501,507,10001 +!full3d_frequency = 5 +/ +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +fix_tvar_top = .false. +fix_tvar_bottom = .false. +fix_dtdr_top = .true. +fix_dtdr_bottom = .true. +dtdr_top_file = '../bench_dtdr_rmax_init' +dtdr_bottom_file = '../bench_dtdr_rmin_init' +fix_chivar_a_top(1) = .false. +fix_chivar_a_bottom(1) = .false. +fix_dchidr_a_top(1) = .true. +fix_dchidr_a_bottom(1) = .true. +dchidr_a_top_file(1) = '../bench_dtdr_rmax_init' +dchidr_a_bottom_file(1) = '../bench_dtdr_rmin_init' +/ +&Initial_Conditions_Namelist +init_type = 8 +t_init_file = '../bench_t_init' +chi_a_init_file(1) = '../bench_t_init' +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 0.0 +chi_a_rayleigh_number(1) = -1.0d5 +Prandtl_Number = 1.0d0 +chi_a_prandtl_number(1) = 1.0d0 +reference_type = 1 +heating_type = 2 ! overwritten by reference state +gravity_power = 1.0d0 ! g ~ radius +custom_reference_file = 'ref_state.dat' +with_custom_reference = .true. +with_custom_constants = 10 +with_custom_functions = 6 +chi_a_custom_reference_file = 'chi_ref_state.dat' +with_custom_chi_a_constants(1,1) = 3 +with_custom_chi_a_functions(1,1) = 3 +/ +&Transport_Namelist +/ diff --git a/tests/chi_custom_reference/chi.split/gen_ref_state.py b/tests/chi_custom_reference/chi.split/gen_ref_state.py new file mode 100644 index 00000000..800cd6f9 --- /dev/null +++ b/tests/chi_custom_reference/chi.split/gen_ref_state.py @@ -0,0 +1,45 @@ +import reference_tools as rt +import numpy +# Grid Parameters +nr = 192 # Number of radial points - this is intentionally not the same as main_input as I want to test the interpolation of the reference state + +#aspect ratio +ar = 0.35e0 +# shell depth +d=1.e0 + +#outer radial boundary +ro=d/(1-ar) + +#inner radial boundary +ri=ar*ro + +#radial grid +radius=numpy.linspace(ri,ro,nr) + +ones = numpy.ones(nr,dtype='float64') +Q0 = 10.0 +dTrefdr = -ro*ri*(radius**-2) + +# Initialize an equation_coefficients structure consistent with the case 0 benchmark test +my_ref = rt.equation_coefficients(radius) + +my_ref.set_function(ones,'heating') # heating function +my_ref.set_constant(Q0,'luminosity') # multiplies the Heating function + +my_ref.set_function(dTrefdr,'ds_dr') +my_ref.set_constant(1.0, 'dsdr_scale') + +my_ref.write('ref_state.dat') # Write the data file + +chi_ref = rt.scalar_equation_coefficients(1, radius) + +chi_ref.set_function(ones,'source_chi') +chi_ref.set_constant(Q0, 'source_chi_scale') + +chi_ref.set_function(dTrefdr, 'dchirefdr') +chi_ref.set_constant(1.0, 'dchirefdr_scale') + +chi_ref.write('chi_ref_state.dat') + + diff --git a/tests/chi_custom_reference/chi.split/main_input b/tests/chi_custom_reference/chi.split/main_input new file mode 100644 index 00000000..0fa3d122 --- /dev/null +++ b/tests/chi_custom_reference/chi.split/main_input @@ -0,0 +1,79 @@ +&problemsize_namelist + n_r = 48 + n_theta = 64 + nprow = 2 + npcol = 2 + aspect_ratio = 0.35d0 + shell_depth = 1.0d0 +/ +&numerical_controls_namelist +/ +&physical_controls_namelist + benchmark_mode = 0 + rotation = .True. + magnetism = .false. + viscous_heating = .false. + ohmic_heating = .false. + n_active_scalars = 1 + n_passive_scalars = 0 +/ +&temporal_controls_namelist + max_time_step = 7.5e-5 + max_iterations = 17 + checkpoint_interval = 10000 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist +globalavg_values = 401,501,503,505,507,509,511,513,515,517,519,521,523,10001,10002,10003,10004,10005,10006,10007,10008,10009,10010,10011,10012 +globalavg_frequency = 17 +globalavg_nrec = 1 +!full3d_values = 1,2,3,501,507,10001 +!full3d_frequency = 10 +/ +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +fix_tvar_top = .false. +fix_tvar_bottom = .false. +fix_dtdr_top = .true. +fix_dtdr_bottom = .true. +dtdr_top_file = '../bench_dtpdr_rmax_init' +dtdr_bottom_file = '../bench_dtpdr_rmin_init' +fix_chivar_a_top(1) = .false. +fix_chivar_a_bottom(1) = .false. +fix_dchidr_a_top(1) = .true. +fix_dchidr_a_bottom(1) = .true. +dchidr_a_top_file(1) = '../bench_dtpdr_rmax_init' +dchidr_a_bottom_file(1) = '../bench_dtpdr_rmin_init' +/ +&Initial_Conditions_Namelist +init_type = 8 +t_init_file = '../bench_tp_init' +chi_a_init_file(1) = '../bench_tp_init' +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 0.0 +chi_a_rayleigh_number(1) = -1.0d5 +Prandtl_Number = 1.0d0 +chi_a_prandtl_number(1) = 1.0d0 +reference_type = 1 +heating_type = 2 ! overwritten by reference state +gravity_power = 1.0d0 ! g ~ radius +custom_reference_file = 'ref_state.dat' +with_custom_reference = .true. +with_custom_constants = 10, 11 +with_custom_functions = 6, 14 +chi_a_custom_reference_file = 'chi_ref_state.dat' +with_custom_chi_a_constants(1,1) = 1 +with_custom_chi_a_constants(1,2) = 3 +with_custom_chi_a_functions(1,1) = 1 +with_custom_chi_a_functions(1,2) = 3 +/ +&Transport_Namelist +/ diff --git a/tests/chi_custom_reference/run_test.sh b/tests/chi_custom_reference/run_test.sh new file mode 100755 index 00000000..bb9961a1 --- /dev/null +++ b/tests/chi_custom_reference/run_test.sh @@ -0,0 +1,62 @@ +#!/usr/bin/env bash + +set -x + +cd tests/chi_custom_reference + +# generate initial conditions +# ../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o const_chi_init \ +# -e '1.0' + +# ../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o diff_t_init \ +# -e 'rmax*rmin/radius - rmin' + +../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_t_init \ + -e 'import numpy as np; x = 2*radius - rmin - rmax; rmax*rmin/radius - rmin + 210*0.1*(1 - 3*x*x + 3*(x**4) - x**6)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + +# ../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_dtdr_init \ +# -e 'import numpy as np; x = 2*radius - rmin - rmax; xp = 2; -rmax*rmin*(radius**-2) + 210*0.1*(- 6*x*xp + 12*(x**3)*xp - 6*(x**5)*xp)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + +../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_dtdr_rmin_init \ + -e 'import numpy as np; x = rmin - rmax; xp = 2; -rmax/rmin + 210*0.1*(- 6*x*xp + 12*(x**3)*xp - 6*(x**5)*xp)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + +../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_dtdr_rmax_init \ + -e 'import numpy as np; x = rmax - rmin; xp = 2; -rmin/rmax + 210*0.1*(- 6*x*xp + 12*(x**3)*xp - 6*(x**5)*xp)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + +../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_tp_init \ + -e 'import numpy as np; x = 2*radius - rmin - rmax; 210*0.1*(1 - 3*x*x + 3*(x**4) - x**6)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + +# ../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_dtpdr_init \ +# -e 'import numpy as np; x = 2*radius - rmin - rmax; xp = 2; 210*0.1*(- 6*x*xp + 12*(x**3)*xp - 6*(x**5)*xp)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + +../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_dtpdr_rmin_init \ + -e 'import numpy as np; x = rmin - rmax; xp = 2; 210*0.1*(- 6*x*xp + 12*(x**3)*xp - 6*(x**5)*xp)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + +../../pre_processing/rayleigh_spectral_input.py -ar 0.35 -sd 1.0 -nt 64 -nr 48 -o bench_dtpdr_rmax_init \ + -e 'import numpy as np; x = rmax - rmin; xp = 2; 210*0.1*(- 6*x*xp + 12*(x**3)*xp - 6*(x**5)*xp)*(np.sin(theta)**4)*np.cos(4*phi)/np.sqrt(17920*np.pi)' + + +cd T.full +PYTHONPATH=../../../post_processing:../../../pre_processing:$PYTHONPATH python3 gen_ref_state.py +mpirun -np 4 $RAYLEIGH_TEST_MPI_PARAMS ../../../bin/rayleigh.dbg +cd .. + +cd T.split +PYTHONPATH=../../../post_processing:../../../pre_processing:$PYTHONPATH python3 gen_ref_state.py +mpirun -np 4 $RAYLEIGH_TEST_MPI_PARAMS ../../../bin/rayleigh.dbg +cd .. + +cd chi.full +PYTHONPATH=../../../post_processing:../../../pre_processing:$PYTHONPATH python3 gen_ref_state.py +mpirun -np 4 $RAYLEIGH_TEST_MPI_PARAMS ../../../bin/rayleigh.dbg +cd .. + +cd chi.split +PYTHONPATH=../../../post_processing:../../../pre_processing:$PYTHONPATH python3 gen_ref_state.py +mpirun -np 4 $RAYLEIGH_TEST_MPI_PARAMS ../../../bin/rayleigh.dbg +cd .. + + +# after both versions have run, we test the output for errors +PYTHONPATH=../../post_processing:../../pre_processing:$PYTHONPATH python3 test_output.py + diff --git a/tests/chi_custom_reference/test_output.py b/tests/chi_custom_reference/test_output.py new file mode 100644 index 00000000..b8cee115 --- /dev/null +++ b/tests/chi_custom_reference/test_output.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 + +from rayleigh_diagnostics import G_Avgs, build_file_list +import numpy as np +import os +import sys + +#aspect ratio +ar = 0.35e0 +# shell depth +d=1.e0 + +#outer radial boundary +rmax=d/(1-ar) + +#inner radial boundary +rmin=ar*rmax + +vals = { + "kesgav" : (401, 0.0), + "Tsgav" : (501, (rmax**3*rmin/6 - rmax*rmin**3/2 + rmin**4/3)*3/(rmax**3-rmin**3)), + "Tpsgav" : (503, 0.0), + "Tmsgav" : (505, (rmax**3*rmin/6 - rmax*rmin**3/2 + rmin**4/3)*3/(rmax**3-rmin**3)), + "Tdrsgav" : (507, -rmax*rmin*(rmax-rmin)*3/(rmax**3-rmin**3)), + "Tpdrsgav" : (509, 0.0), + "Tmdrsgav" : (511, -rmax*rmin*(rmax-rmin)*3/(rmax**3-rmin**3)), + "Tdtsgav" : (513, 0.0), + "Tpdtsgav" : (515, 0.0), + "Tmdtsgav" : (517, 0.0), + "Tdpsgav" : (519, 0.0), + "Tpdpsgav" : (521, 0.0), + "Tmdpsgav" : (523, 0.0), + "chia1sgav" : (10001, (rmax**3*rmin/6 - rmax*rmin**3/2 + rmin**4/3)*3/(rmax**3-rmin**3)), + "chia1psgav" : (10002, 0.0), + "chia1msgav" : (10003, (rmax**3*rmin/6 - rmax*rmin**3/2 + rmin**4/3)*3/(rmax**3-rmin**3)), + "chia1drsgav" : (10004, -rmax*rmin*(rmax-rmin)*3/(rmax**3-rmin**3)), + "chia1pdrsgav" : (10005, 0.0), + "chia1mdrsgav" : (10006, -rmax*rmin*(rmax-rmin)*3/(rmax**3-rmin**3)), + "chia1dtsgav" : (10007, 0.0), + "chia1pdtsgav" : (10008, 0.0), + "chia1mdtsgav" : (10009, 0.0), + "chia1dpsgav" : (10010, 0.0), + "chia1pdpsgav" : (10011, 0.0), + "chia1mdpsgav" : (10012, 0.0), + } + +def check_results(dirs, tol=1.e-10): + error = False + results = { + "times" : [], + } + for k in vals.keys(): results[k] = [] + for d in dirs: + files = build_file_list(0,1000000,path=os.path.join(d,'G_Avgs')) + a = G_Avgs(files[-1], path='') + results["times"].append(a.time[-1]) + for k,(v,off) in vals.items(): + try: + val = a.vals[-1, a.lut[v]] + if d.endswith('full'): val -= off + results[k].append(val) + except IndexError: + import ipdb; ipdb.set_trace() + + for k,v in results.items(): + print(k+":\t", v) + if np.any(np.abs(v-v[0]) > tol): + print("ERROR: different "+k+" produced between runs (within a tolerance of "+repr(tol)+")!") + error = True + if k.startswith("chi"): + Tk = "T"+ k[5:] + if np.any(np.abs(np.asarray(v)-np.asarray(results[Tk])) > tol): + print("ERROR: different "+k+" and "+Tk+" produced (with a tolerance of "+repr(tol)+")!") + error = True + + return error + +error = check_results(["T.full", "chi.full"]) or check_results(["T.split", "chi.split"]) + +if error: sys.exit(1) +sys.exit(0) +