Adding arbitrary (up to 50 of each) active and scalar fields

.github/workflows/main.yml

@@ -49,6 +49,10 @@ jobs:
         cd "$GITHUB_WORKSPACE"
         sh ./tests/generic_input/run_test.sh
 
+        # chi scalar test
+        cd "$GITHUB_WORKSPACE"
+        sh ./tests/chi_scalar/run_test.sh
+
         git diff > changes.diff
         git diff --exit-code --name-only
 

doc/source/User_Guide/index.rst

@@ -12,4 +12,5 @@ User Guide
    analyze_output
    ../../../CONTRIBUTING.md
    troubleshooting
-   references
+   references
+   under_development

doc/source/User_Guide/under_development.rst

@@ -0,0 +1,41 @@
+.. raw:: latex
+
+   \clearpage
+
+.. _under_development:
+
+Under Development
+=================
+
+.. _scalar_fields:
+
+Arbitrary Scalar Fields
+-----------------------
+
+Rayleigh can solve for additional active, :math:`\chi_{a_i}`, (coupled to the momentum equation through buoyancy) or
+passive, :math:`\chi_{p_i}`, scalar fields (where :math:`i` can range up to 50 for each type of scalar).  Both types of field follow a simple advection-diffusion equation:
+
+.. math::
+   :label: scalar_evol
+
+   \frac{\partial \chi_{a,p_i}}{\partial t}  + \boldsymbol{v}\cdot\boldsymbol{\nabla}\chi_{a,p_i}  = 0
+
+The number of each type of field can be set using, e.g.:
+
+::
+
+   &physical_controls_namelist
+    n_active_scalars = 2
+    n_passive_scalars = 2
+   /
+
+Other model parameters follow the same convention as temperature but using the prefix `chi_a` or `chi_p` for active and passive
+scalars respectively.
+
+See `tests/chi_scalar` for example input files.
+
+
+
+
+
+

post_processing/rayleigh_diagnostics.py

@@ -24,7 +24,7 @@
 import glob
 from collections import OrderedDict
 
-maxq = 4000
+maxq = 30000
 
 def get_lut(quantities):
     """return the lookup table based on the quantity codes"""

src/Diagnostics/Diagnostics_Base.F90

@@ -89,6 +89,7 @@ Module Diagnostics_Base
 
     include "turbKE_codes.F"
     include "axial_field_codes.F"
+    include "scalars_field_codes.F"
 
 
     !///////////////////////////////////

src/Diagnostics/Diagnostics_Interface.F90

@@ -61,6 +61,8 @@ Module Diagnostics_Interface
 
     Use Diagnostics_Custom
 
+    Use Diagnostics_Scalars
+
     Implicit None
 
 
@@ -228,6 +230,8 @@ Subroutine PS_Output(buffer,iteration, current_time)
                 Call Compute_Misc_Diagnostics(buffer)
                 Call Custom_Hydro_Diagnostics(buffer)
 
+                Call Compute_Scalars(buffer)
+
                 !////// Magnetic Quantities
                 If (magnetism) Then
                     Call Compute_BField_Components(buffer)

src/Diagnostics/Diagnostics_Scalars.F90

@@ -0,0 +1,145 @@
+!
+!  Copyright (C) 2018 by the authors of the RAYLEIGH code.
+!
+!  This file is part of RAYLEIGH.
+!
+!  RAYLEIGH is free software; you can redistribute it and/or modify
+!  it under the terms of the GNU General Public License as published by
+!  the Free Software Foundation; either version 3, or (at your option)
+!  any later version.
+!
+!  RAYLEIGH is distributed in the hope that it will be useful,
+!  but WITHOUT ANY WARRANTY; without even the implied warranty of
+!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+!  GNU General Public License for more details.
+!
+!  You should have received a copy of the GNU General Public License
+!  along with RAYLEIGH; see the file LICENSE.  If not see
+!  <http://www.gnu.org/licenses/>.
+!
+
+#include "indices.F"
+Module Diagnostics_Scalars
+    Use Diagnostics_Base
+    Use Controls, Only: n_active_scalars, n_passive_scalars
+    Implicit None
+
+Contains
+
+    Subroutine Compute_Scalars(buffer)
+        Implicit None
+        Real*8, Intent(InOut) :: buffer(1:,my_r%min:,my_theta%min:,1:)
+        Integer :: r,k, t, ii, ind, scoff, n_scalars
+        Integer :: chivar, dchidr, dchidt, dchidp, d2chidr2
+
+        n_scalars = n_active_scalars + n_passive_scalars
+        Do ii = 1, n_scalars
+            if (ii .le. n_active_scalars) then
+                 ind   = ii
+                 scoff = a_scalar_offset + (ind-1)*scalar_skip
+                 chivar   = chiavar(ind)
+                 dchidr   = dchiadr(ind)
+                 dchidt   = dchiadt(ind)
+                 dchidp   = dchiadp(ind)
+                 d2chidr2 = d2chiadr2(ind)
+            else
+                 ind = ii - n_active_scalars
+                 scoff = p_scalar_offset + (ind-1)*scalar_skip
+                 chivar   = chipvar(ind)
+                 dchidr   = dchipdr(ind)
+                 dchidt   = dchipdt(ind)
+                 dchidp   = dchipdp(ind)
+                 d2chidr2 = d2chipdr2(ind)
+            end if
+
+            !  Chi: field
+            If (compute_quantity(chi+scoff)) Then
+                DO_PSI
+                    qty(PSI) = buffer(PSI,chivar)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_p+scoff)) Then
+                DO_PSI
+                    qty(PSI) = fbuffer(PSI,chivar)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_m+scoff)) Then
+                DO_PSI
+                    qty(PSI) = m0_values(PSI2,chivar)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            ! Chi:  radial derivatives
+            If (compute_quantity(chi_dr+scoff)) Then
+                DO_PSI
+                    qty(PSI) = buffer(PSI,dchidr)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_p_dr+scoff)) Then
+                DO_PSI
+                    qty(PSI) = fbuffer(PSI,dchidr)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_m_dr+scoff)) Then
+                DO_PSI
+                    qty(PSI) = m0_values(PSI2,dchidr)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            ! Chi:  theta derivatives
+            If (compute_quantity(chi_dtheta+scoff)) Then
+                DO_PSI
+                    qty(PSI) = buffer(PSI,dchidt)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_p_dtheta+scoff)) Then
+                DO_PSI
+                    qty(PSI) = fbuffer(PSI,dchidt)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_m_dtheta+scoff)) Then
+                DO_PSI
+                    qty(PSI) = m0_values(PSI2,dchidt)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            ! Chi:  phi derivatives
+            If (compute_quantity(chi_dphi+scoff)) Then
+                DO_PSI
+                    qty(PSI) = buffer(PSI,dchidp)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_p_dphi+scoff)) Then
+                DO_PSI
+                    qty(PSI) = fbuffer(PSI,dchidp)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+
+            If (compute_quantity(chi_m_dphi+scoff)) Then
+                DO_PSI
+                    qty(PSI) = m0_values(PSI2,dchidp)
+                END_DO
+                Call Add_Quantity(qty)
+            Endif
+        Enddo
+    End Subroutine Compute_Scalars
+
+End Module Diagnostics_Scalars

src/Diagnostics/scalars_field_codes.F

@@ -0,0 +1,56 @@
+!
+!  Copyright (C) 2018 by the authors of the RAYLEIGH code.
+!
+!  This file is part of RAYLEIGH.
+!
+!  RAYLEIGH is free software; you can redistribute it and/or modify
+!  it under the terms of the GNU General Public License as published by
+!  the Free Software Foundation; either version 3, or (at your option)
+!  any later version.
+!
+!  RAYLEIGH is distributed in the hope that it will be useful,
+!  but WITHOUT ANY WARRANTY; without even the implied warranty of
+!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+!  GNU General Public License for more details.
+!
+!  You should have received a copy of the GNU General Public License
+!  along with RAYLEIGH; see the file LICENSE.  If not see
+!  <http://www.gnu.org/licenses/>.
+!
+
+    !//////////////////////////////////////////////////////////////////
+    !  Composition/Scalar (Active)
+    Integer, Parameter :: scalar_skip     = 200
+    Integer, Parameter :: a_scalar_offset = 10000 ! :OFFSET CODE:
+    Integer, Parameter :: p_scalar_offset = 20000 ! :OFFSET CODE:
+
+    !------------ Fields ---------------------!
+    ! Full
+    Integer, Parameter :: chi   = 1 !  :tex:  $\chi$
+    ! Fluctuating
+    Integer, Parameter :: chi_p = 2 !  :tex:  $\chi'$
+    ! Mean
+    Integer, Parameter :: chi_m = 3 !  :tex:  $\overline{\chi}$
+
+    !------------ Radial Derivatives --------------!
+    !  Full
+    Integer, Parameter :: chi_dr   = 4  ! :tex:  $\frac{\partial \chi}{\partial r}$
+    ! Fluctuating
+    Integer, Parameter :: chi_p_dr = 5  ! :tex:  $\frac{\partial \chi'}{\partial r}$
+    ! Mean
+    Integer, Parameter :: chi_m_dr = 6 ! :tex:  $\frac{\partial \overline{\chi}}{\partial r}$
+
+    !------------ Theta Derivatives ---------------!
+    Integer, Parameter :: chi_dtheta   = 7  ! :tex:  $\frac{\partial \chi}{\partial \theta}$
+
+    Integer, Parameter :: chi_p_dtheta = 8  ! :tex:  $\frac{\partial \chi'}{\partial \theta}$
+
+    Integer, Parameter :: chi_m_dtheta = 9  ! :tex:  $\frac{\partial \overline{\chi}}{\partial \theta}$
+
+    !------------ Phi Derivatives -----------------!
+    Integer, Parameter :: chi_dphi     = 10  ! :tex:  $\frac{\partial \chi}{\partial \phi}$
+
+    Integer, Parameter :: chi_p_dphi   = 11  ! :tex:  $\frac{\partial \chi'}{\partial \phi}$
+
+    Integer, Parameter :: chi_m_dphi   = 12  ! :tex:  $\frac{\partial \overline{\chi}}{\partial \phi}$
+

src/IO/Spherical_IO.F90

@@ -49,7 +49,7 @@ Module Spherical_IO
 
 
     !////////////////////////////////////////////
-    Integer, Parameter :: nqmax=4000, nshellmax=2048, nmeridmax=8192, nmodemax=2048
+    Integer, Parameter :: nqmax=30000, nshellmax=2048, nmeridmax=8192, nmodemax=2048
     Integer, Parameter :: nprobemax=4096
     Integer, Parameter :: endian_tag = 314      ! first 4 bits of each diagnostic file - used for assessing endianness on read-in
     Integer, Parameter :: reallybig = 90000000

src/Makefile.fdeps

@@ -10,20 +10,21 @@ Controls.o : Controls.F90 BufferedOutput.o
 Diagnostics_ADotGradB.o : Diagnostics_ADotGradB.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Angular_Momentum.o : Diagnostics_Angular_Momentum.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Axial_Field.o : Diagnostics_Axial_Field.F90 indices.F Diagnostics_Base.o 
-Diagnostics_Base.o : Diagnostics_Base.F90 axial_field_codes.F turbKE_codes.F me_equation_codes.F ke_equation_codes.F amom_equation_codes.F induction_equation_codes.F thermal_equation_codes.F momentum_equation_codes.F magnetic_energy_codes.F current_density_codes.F magnetic_field_codes.F thermal_energy_codes.F thermal_field_codes.F kinetic_energy_codes.F vorticity_field_codes.F mass_flux_codes.F velocity_field_codes.F indices.F PDE_Coefficients.o Math_Constants.o Fields.o Spherical_IO.o ProblemSize.o 
+Diagnostics_Base.o : Diagnostics_Base.F90 scalars_field_codes.F axial_field_codes.F turbKE_codes.F me_equation_codes.F ke_equation_codes.F amom_equation_codes.F induction_equation_codes.F thermal_equation_codes.F momentum_equation_codes.F magnetic_energy_codes.F current_density_codes.F magnetic_field_codes.F thermal_energy_codes.F thermal_field_codes.F kinetic_energy_codes.F vorticity_field_codes.F mass_flux_codes.F velocity_field_codes.F indices.F PDE_Coefficients.o Math_Constants.o Fields.o Spherical_IO.o ProblemSize.o 
 Diagnostics_Current_Density.o : Diagnostics_Current_Density.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Custom.o : Diagnostics_Custom.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Energies.o : Diagnostics_Energies.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Induction.o : Diagnostics_Induction.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
 Diagnostics_Inertial_Forces.o : Diagnostics_Inertial_Forces.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
-Diagnostics_Interface.o : Diagnostics_Interface.F90 indices.F Diagnostics_Custom.o Diagnostics_Miscellaneous.o Diagnostics_Poynting_Flux.o Diagnostics_Induction.o Diagnostics_Axial_Field.o Diagnostics_TurbKE_Budget.o Diagnostics_KE_Flux.o Diagnostics_Lorentz_Forces.o Diagnostics_Angular_Momentum.o Diagnostics_Inertial_Forces.o Diagnostics_Linear_Forces.o Diagnostics_Current_Density.o Diagnostics_Vorticity_Field.o Diagnostics_Thermal_Equation.o Diagnostics_Thermal_Energies.o Diagnostics_Thermodynamic_Gradients.o Diagnostics_Energies.o Diagnostics_Magnetic_Field.o Diagnostics_Velocity_Field.o Diagnostics_Mean_Correction.o Diagnostics_Second_Derivatives.o Diagnostics_Base.o Math_Constants.o PDE_Coefficients.o Legendre_Polynomials.o Fields.o Spherical_IO.o Controls.o ProblemSize.o 
+Diagnostics_Interface.o : Diagnostics_Interface.F90 indices.F Diagnostics_Scalars.o Diagnostics_Custom.o Diagnostics_Miscellaneous.o Diagnostics_Poynting_Flux.o Diagnostics_Induction.o Diagnostics_Axial_Field.o Diagnostics_TurbKE_Budget.o Diagnostics_KE_Flux.o Diagnostics_Lorentz_Forces.o Diagnostics_Angular_Momentum.o Diagnostics_Inertial_Forces.o Diagnostics_Linear_Forces.o Diagnostics_Current_Density.o Diagnostics_Vorticity_Field.o Diagnostics_Thermal_Equation.o Diagnostics_Thermal_Energies.o Diagnostics_Thermodynamic_Gradients.o Diagnostics_Energies.o Diagnostics_Magnetic_Field.o Diagnostics_Velocity_Field.o Diagnostics_Mean_Correction.o Diagnostics_Second_Derivatives.o Diagnostics_Base.o Math_Constants.o PDE_Coefficients.o Legendre_Polynomials.o Fields.o Spherical_IO.o Controls.o ProblemSize.o 
 Diagnostics_KE_Flux.o : Diagnostics_KE_Flux.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Linear_Forces.o : Diagnostics_Linear_Forces.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Lorentz_Forces.o : Diagnostics_Lorentz_Forces.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Magnetic_Field.o : Diagnostics_Magnetic_Field.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Mean_Correction.o : Diagnostics_Mean_Correction.F90 indices.F Diagnostics_ADotGradB.o Spherical_IO.o Diagnostics_Base.o 
 Diagnostics_Miscellaneous.o : Diagnostics_Miscellaneous.F90 indices.F Diagnostics_Base.o 
 Diagnostics_Poynting_Flux.o : Diagnostics_Poynting_Flux.F90 indices.F Diagnostics_Base.o 
+Diagnostics_Scalars.o : Diagnostics_Scalars.F90 indices.F Controls.o Diagnostics_Base.o 
 Diagnostics_Second_Derivatives.o : Diagnostics_Second_Derivatives.F90 indices.F Spectral_Derivatives.o Structures.o Diagnostics_Base.o 
 Diagnostics_Thermal_Energies.o : Diagnostics_Thermal_Energies.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
 Diagnostics_Thermal_Equation.o : Diagnostics_Thermal_Equation.F90 indices.F Diagnostics_ADotGradB.o Diagnostics_Base.o 
@@ -89,6 +90,7 @@ magnetic_field_codes.o : magnetic_field_codes.F
 mass_flux_codes.o : mass_flux_codes.F 
 me_equation_codes.o : me_equation_codes.F 
 momentum_equation_codes.o : momentum_equation_codes.F 
+scalars_field_codes.o : scalars_field_codes.F 
 thermal_energy_codes.o : thermal_energy_codes.F 
 thermal_equation_codes.o : thermal_equation_codes.F 
 thermal_field_codes.o : thermal_field_codes.F