diff --git a/doc/source/User_Guide/newtonian_cooling.txt b/doc/source/User_Guide/newtonian_cooling.txt new file mode 100644 index 00000000..04c8ce19 --- /dev/null +++ b/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' + / + + diff --git a/doc/source/User_Guide/under_development.rst b/doc/source/User_Guide/under_development.rst index 4ad4b9c3..d453af06 100644 --- a/doc/source/User_Guide/under_development.rst +++ b/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 + diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index d7f1d43c..9e9d767e 100755 --- a/src/Physics/Controls.F90 +++ b/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 !/////////////////////////////////////////////////////////////////////////// diff --git a/src/Physics/Initial_Conditions.F90 b/src/Physics/Initial_Conditions.F90 index 92e0638b..6be28159 100755 --- a/src/Physics/Initial_Conditions.F90 +++ b/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 diff --git a/src/Physics/Sphere_Driver.F90 b/src/Physics/Sphere_Driver.F90 index 49196b54..26076e73 100755 --- a/src/Physics/Sphere_Driver.F90 +++ b/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 diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index 992e8a7d..7d585696 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/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