Pleiades Installation Instruction Update

doc/source/User_Guide/getting_started.rst

@@ -415,9 +415,8 @@ Installation on NASA's Pleiades cluster is similarly straightforward.  After clo
 .. code-block:: bash
 
    module purge
-   module load comp-intel
-   module load mpi-hpe
-   ./configure --FC=mpif90 --CC=icc  # select 'ALL'
+   module load comp-intel mpi-hpe
+   ./configure --FC=mpif90 --CC=icc -mkl  # select 'ALL'
    make -j
    make install
 

doc/source/User_Guide/newtonian_cooling.txt

@@ -0,0 +1,42 @@
+.. _newtonian_cooling:
+
+Newtonian Cooling
+-----------------------
+
+We have added an initial implementation of Newtonian cooling to the code that adds a term to the thermal energy equation of the form
+
+.. math::
+   :label: newton
+
+   \frac{\partial \Theta}{\partial t}  = \frac{\Delta\Theta_{\mathrm{eq}}-\Theta}{\tau_\mathrm{eq}}.
+
+Here, :math:`\Delta\Theta_\mathrm{eq}` is a target temperature/entropy variation about the background state and :math:`\tau_\mathrm{eq}` is the cooling timescale.  Newtonian cooling can be turned on by setting ``newtonian_cooling=.true.`` in the ``physical_controls_namelist``, and the cooling time is similarly controlled by specifying the value of ``newtonian_cooling_time``.  
+
+At present, :math:`\Delta\Theta_\mathrm{eq}` is allowed to take one of two forms.  These are controlled by setting the ``newtonian_cooling_type`` variable in ``physical_controls_namelist``.  A value of 1 yields
+
+.. math::
+   :label: ncteq1
+
+   \Delta\Theta_\mathrm{eq} = A,
+
+and a value of 2 yields
+
+.. math::
+   :label: ncteq2
+
+   \Delta\Theta_\mathrm{eq} = A\mathrm{cos}(\theta)\mathrm{sin}(\phi)\mathrm{log}(\hat{P}/\hat{P}_\mathrm{bottom})/\mathrm{log}(\hat{P}_\mathrm{top}/\hat{P}_\mathrm{bottom}),
+
+where :math:`\hat{P}` is the reference state pressure.
+
+The amplitude :math:`A` is controlled by setting the variable  ``newtonian_cooling_tvar_amp``.   As an example, to use Newtonian cooling, one might add the following four lines to the ``physical_controls_namelist``.
+
+::
+
+   &physical_controls_namelist
+    newtonian_cooling = .true.
+    newtonian_cooling_type = 1
+    newtonian_cooling_time = 1.0d0
+    newtonian_cooling_tvar_amp = 1.0d0
+   /
+
+

doc/source/User_Guide/under_development.rst

@@ -16,8 +16,6 @@ Arbitrary Scalar Fields
 
 .. include:: arbitrary_scalar_fields.txt 
 
-
-
-
+.. include:: newtonian_cooling.txt
 
 

src/Physics/Controls.F90

@@ -69,16 +69,25 @@ Module Controls
     Integer :: n_active_scalars = 0         ! number of active scalar fields
     Integer :: n_passive_scalars = 0        ! number of passive scalar fields
 
+
     ! --- This flag determines if the code is run in benchmark mode
     !     0 (default) is no benchmarking.  1-5 are various accuracy benchmarks (see documentation)
     Integer :: benchmark_mode = 0
     Integer :: benchmark_integration_interval = -1 ! manual override of integration_interval
     Integer :: benchmark_report_interval = -1      ! and report interval in Benchmarking.F90 (for debugging)
 
+    ! --- Newtonian Cooling Variables
+    Logical :: newtonian_cooling = .false.  ! Turn newtonian_cooling on/off
+    Integer :: newtonian_cooling_type = 1
+    Real*8  :: newtonian_cooling_time = 1.0d22
+    Real*8  :: newtonian_cooling_tvar_amp = 0.0d0
+
     Namelist /Physical_Controls_Namelist/ magnetism, nonlinear, rotation, lorentz_forces, &
                 & viscous_heating, ohmic_heating, advect_reference_state, benchmark_mode, &
                 & benchmark_integration_interval, benchmark_report_interval, &
-                & momentum_advection, inertia, n_active_scalars, n_passive_scalars
+                & momentum_advection, inertia, n_active_scalars, n_passive_scalars, &
+                & newtonian_cooling, newtonian_cooling_type, newtonian_cooling_time, &
+                & newtonian_cooling_tvar_amp
 
     !///////////////////////////////////////////////////////////////////////////
     !   Temporal Controls

src/Physics/Sphere_Driver.F90

@@ -21,7 +21,7 @@
 Module Sphere_Driver
     Use ClockInfo
     Use Sphere_Hybrid_Space,   Only : rlm_spacea, rlm_spaceb, hybrid_init
-    Use Sphere_Physical_Space, Only : physical_space, ohmic_heating_coeff
+    Use Sphere_Physical_Space, Only : physical_space, ohmic_heating_coeff, physical_space_init
     Use Sphere_Spectral_Space, Only : post_solve, advancetime, ctemp, post_solve_FD
     Use Diagnostics_Interface, Only : Reboot_Diagnostics
     Use Spherical_IO, Only : time_to_output
@@ -130,6 +130,7 @@ Subroutine Main_Loop_Sphere()
 
 
         Call Hybrid_Init()
+        Call Physical_Space_Init()
         Call StopWatch(loop_time)%StartClock()
         max_time_seconds = 60*max_time_minutes
 

src/Physics/Sphere_Physical_Space.F90

@@ -37,7 +37,55 @@ Module Sphere_Physical_Space
     Use Benchmarking, Only : benchmark_checkup
     Implicit None
 
+    Real*8, Allocatable :: tvar_eq(:,:,:)
+
 Contains
+
+    Subroutine Physical_Space_Init()
+        Implicit None
+        Integer :: k, r, t
+        Real*8, Allocatable :: press(:)
+
+
+
+        ! Any persistant arrays needs for physical space routines can be
+        ! initialized here.
+        If (newtonian_cooling) Then
+            Allocate(tvar_eq(1:n_phi, my_r%min:my_r%max, my_theta%min:my_theta%max))
+            tvar_eq(:,:,:) = 0.0d0
+
+            If (newtonian_cooling_type .eq. 1) Then
+                ! No angular variation
+                If (my_rank .eq. 0) Write(6,*) 'Newtonian cooling is active.  Type = 1'
+                Do t = my_theta%min, my_theta%max
+                    Do r = my_r%min, my_r%max
+                        Do k =1, n_phi
+                            tvar_eq(k,r,t) = newtonian_cooling_tvar_amp
+                        Enddo
+                    Enddo
+                Enddo
+            Endif
+
+
+            If (newtonian_cooling_type .eq. 2) Then
+                ! Angular variation (ell=1,m=1, motivated by hot Jupiters)
+
+                Allocate(press(1:N_R))
+                press = ref%density*ref%temperature
+                If (my_rank .eq. 0) Write(6,*) 'Newtonian cooling is active.  Type = 2'
+                Do t = my_theta%min, my_theta%max
+                    Do r = my_r%min, my_r%max
+                        Do k =1, n_phi
+                            tvar_eq(k,r,t) = newtonian_cooling_tvar_amp*sintheta(t)*sinphi(k)
+                            tvar_eq(k,r,t) = tvar_eq(k,r,t)*log(press(r)/press(N_R))/log(press(1)/press(N_R))
+                        Enddo
+                    Enddo
+                Enddo
+            Endif
+        Endif
+
+    End Subroutine Physical_Space_Init
+
     Subroutine physical_space()
         Implicit None
         Integer :: i
@@ -123,6 +171,7 @@ Subroutine physical_space()
         Call Temperature_Advection()
         Call Volumetric_Heating()
 
+
         do i = 1, n_active_scalars
           Call chi_Advection(chiavar(i), dchiadr(i), dchiadt(i), dchiadp(i))
           Call chi_Source_function(chiavar(i))
@@ -241,8 +290,24 @@ Subroutine Volumetric_Heating()
             Enddo
             !$OMP END PARALLEL DO
         Endif
+
+
+        If (newtonian_cooling) Then
+            ! Added a volumetric heating to the energy equation
+            !$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,tvar) = wsp%p3b(k,r,t,tvar) + &
+                                      (tvar_eq(k,r,t) -wsp%p3b(k,r,t,tvar))/newtonian_cooling_time
+                    Enddo
+                Enddo
+            Enddo
+            !$OMP END PARALLEL DO
+        Endif
     End Subroutine Volumetric_Heating
 
+
     Subroutine chi_Source_Function(chivar)
         Implicit None
         Integer :: chivar