Merge pull request #542 from feathern/cooling_june20

Newtonian Cooling

doc/source/User_Guide/newtonian_cooling.txt

@@ -0,0 +1,41 @@
+.. _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\,f_\mathrm{c}(r),
+
+and a value of 2 yields
+
+.. math::
+   :label: ncteq2
+
+   \Delta\Theta_\mathrm{eq} = A\,f_\mathrm{c}(r)\mathrm{sin}(\theta)\mathrm{sin}(\phi),
+
+where :math:`f_\mathrm{c}(r)` is the radial cooling profile.  It is 1 by default, but the user can specify a file from which to read a custom profile by setting the ``newtonian_cooling_profile_file`` variable in the ``physical_controls_namelist``.  The amplitude :math:`A` is controlled by setting the variable  ``newtonian_cooling_tvar_amp``.   As an example, to use Newtonian cooling, one can and and modify the following lines in 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
+    newtonian_cooling_profile_file = 'my_cooling_profile.dat'
+   /
+
+

doc/source/User_Guide/under_development.rst

@@ -17,6 +17,9 @@ Arbitrary Scalar Fields
 .. include:: arbitrary_scalar_fields.txt 
 
 
+.. include:: newtonian_cooling.txt
+
+
 .. _coupled_bcs:
 
 Coupled Boundary Conditions
@@ -55,3 +58,4 @@ The continuity equation is still a solenoidal constraint, but instead of the mas
 
 
 
+

src/Physics/Controls.F90

@@ -70,16 +70,27 @@ 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
+    Character*120 :: newtonian_cooling_profile_file = '__nothing__'
+    Real*8, Allocatable :: newtonian_cooling_profile(:)
+
     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, &
+                & newtonian_cooling, newtonian_cooling_type, newtonian_cooling_time, &
+                & newtonian_cooling_tvar_amp, newtonian_cooling_profile_file, &
                 & pseudo_incompressible
 
     !///////////////////////////////////////////////////////////////////////////

src/Physics/Initial_Conditions.F90

@@ -216,6 +216,13 @@ Subroutine Initialize_Fields()
                 Endif
             Endif
         Endif
+
+        If ((newtonian_cooling) .and. (newtonian_cooling_profile_file .ne. '__nothing__')) Then
+            Allocate(newtonian_cooling_profile(1:N_R))
+            newtonian_cooling_profile(:) = 0.0d0
+            Call Load_Radial_Profile(newtonian_cooling_profile_file,newtonian_cooling_profile)
+        Endif
+
         ! Fields are now initialized and loaded into the RHS.
         ! We are ready to enter the main loop
         If (my_rank .eq. 0) Then

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,64 @@ 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 :: cooling_profile(:)
+
+        Allocate(cooling_profile(1:N_R))
+        If (newtonian_cooling .and. (newtonian_cooling_profile_file .ne. '__nothing__')) Then
+            cooling_profile(:) = newtonian_cooling_profile(:)
+            If (my_rank .eq. 0) Then
+                Write(6,*) 'Newtonian cooling is active.'
+                Write(6,*) 'Cooling profile set from: ',newtonian_cooling_profile_file
+            Endif
+        Else
+            cooling_profile(:) = 1.0d0
+        Endif
+
+        ! 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 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*newtonian_cooling_profile(r)
+                        Enddo
+                    Enddo
+                Enddo
+            Endif
+
+
+            If (newtonian_cooling_type .eq. 2) Then
+                ! Angular variation (ell=1,m=1, motivated by hot Jupiters)
+
+                If (my_rank .eq. 0) Write(6,*) 'Newtonian cooling 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)*newtonian_cooling_profile(r)
+                        Enddo
+                    Enddo
+                Enddo
+            Endif
+        Endif
+
+        DeAllocate(cooling_profile)
+
+    End Subroutine Physical_Space_Init
+
     Subroutine physical_space()
         Implicit None
         Integer :: i
@@ -123,6 +180,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 +299,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