From 3d7a2b453a20338370c364d0889c5fa0028d3b66 Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Thu, 20 Jun 2024 12:00:46 -0600 Subject: [PATCH 01/10] Initial commit of newtonian cooling modifications. 1) Added relevant variables to Controls.F90 2) Added skeleton loop to Volumetric_Heating routine --- src/Physics/Controls.F90 | 9 ++++++++- src/Physics/Sphere_Physical_Space.F90 | 16 ++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 1b314aae..13d90860 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -69,16 +69,23 @@ 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_timescale = 1.0d22 + 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 !/////////////////////////////////////////////////////////////////////////// ! Temporal Controls diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index 968bf644..20ee857d 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -123,6 +123,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 +242,23 @@ 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)+0.0d0 + Enddo + Enddo + Enddo + !$OMP END PARALLEL DO + Endif End Subroutine Volumetric_Heating + Subroutine chi_Source_Function(chivar) Implicit None Integer :: chivar From da4ac0d101a93fc853da2000ef70fe79f73b5168 Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Thu, 20 Jun 2024 12:09:49 -0600 Subject: [PATCH 02/10] Added physical_space_init subroutine (empty) and call (sphere_driver.F90) --- src/Physics/Sphere_Driver.F90 | 3 ++- src/Physics/Sphere_Physical_Space.F90 | 10 ++++++++++ 2 files changed, 12 insertions(+), 1 deletion(-) 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 20ee857d..49255e74 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -38,6 +38,16 @@ Module Sphere_Physical_Space Implicit None Contains + + Subroutine Physical_Space_Init() + Implicit None + + ! Any persistant arrays needs for physical space routines can be + ! initialized here. + + + End Subroutine Physical_Space_Init + Subroutine physical_space() Implicit None Integer :: i From e3b76329a880c52053780c74bcd54121183e68cd Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Thu, 20 Jun 2024 12:24:44 -0600 Subject: [PATCH 03/10] Completed newtonian_cooling loop in volumetric_heating routine --- src/Physics/Controls.F90 | 2 +- src/Physics/Sphere_Physical_Space.F90 | 13 +++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 13d90860..88316eea 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -79,7 +79,7 @@ Module Controls ! --- Newtonian Cooling Variables Logical :: newtonian_cooling = .false. ! Turn newtonian_cooling on/off Integer :: newtonian_cooling_type = 1 - Real*8 :: newtonian_cooling_timescale = 1.0d22 + Real*8 :: newtonian_cooling_time = 1.0d22 Namelist /Physical_Controls_Namelist/ magnetism, nonlinear, rotation, lorentz_forces, & & viscous_heating, ohmic_heating, advect_reference_state, benchmark_mode, & diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index 49255e74..5091fc5e 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -37,6 +37,8 @@ Module Sphere_Physical_Space Use Benchmarking, Only : benchmark_checkup Implicit None + Real*8, Allocatable :: tvar_eq(:,:,:) + Contains Subroutine Physical_Space_Init() @@ -44,7 +46,13 @@ Subroutine Physical_Space_Init() ! 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 + tvar_eq(:,:,:) = 0.0d0 + Endif + Endif End Subroutine Physical_Space_Init @@ -260,7 +268,8 @@ Subroutine Volumetric_Heating() 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)+0.0d0 + 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 From ea50650c21cf646d7a61a115f1e4561c0e0b9b32 Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Thu, 20 Jun 2024 12:41:12 -0600 Subject: [PATCH 04/10] Initial version of newtonian_cooling is now implemented (but not documented) --- src/Physics/Controls.F90 | 4 +++- src/Physics/Sphere_Physical_Space.F90 | 9 ++++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 88316eea..b8f4cb9c 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -80,12 +80,14 @@ Module Controls Logical :: newtonian_cooling = .false. ! Turn newtonian_cooling on/off Integer :: newtonian_cooling_type = 1 Real*8 :: newtonian_cooling_time = 1.0d22 + Real*8 :: newtonian_delta_tvar_Eq = 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, & - & newtonian_cooling + & newtonian_cooling, newtonian_cooling_type, newtonian_cooling_time, & + & newtonian_delta_tvar_eq !/////////////////////////////////////////////////////////////////////////// ! Temporal Controls diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index 5091fc5e..f0ad3079 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -43,6 +43,7 @@ Module Sphere_Physical_Space Subroutine Physical_Space_Init() Implicit None + Integer :: k, r, t ! Any persistant arrays needs for physical space routines can be ! initialized here. @@ -50,7 +51,13 @@ Subroutine Physical_Space_Init() 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 - tvar_eq(:,:,:) = 0.0d0 + 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_delta_tvar_eq*costheta(t)*sinphi(k) + Enddo + Enddo + Enddo Endif Endif From 48b6e72779c9a3fc1fe9f6aa4dbc66a7883b541c Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Thu, 20 Jun 2024 12:42:52 -0600 Subject: [PATCH 05/10] Added newtonian_cooling_type =2 to distinguish between angular and no angular variation --- src/Physics/Sphere_Physical_Space.F90 | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index f0ad3079..dc5e3ae6 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -50,7 +50,21 @@ Subroutine Physical_Space_Init() 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 + 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_delta_tvar_eq + Enddo + Enddo + Enddo + Endif + + + If (newtonian_cooling_type .eq. 2) Then + ! Angular variation (ell=1,m=1, motivated by hot Jupiters) Do t = my_theta%min, my_theta%max Do r = my_r%min, my_r%max Do k =1, n_phi From eb5015fcb9ff42ea389ae28fc1e3481c520b4118 Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Thu, 20 Jun 2024 13:31:03 -0600 Subject: [PATCH 06/10] Initial documentation added to 'under development' section. Changed amplitude variable name for cooling function. --- doc/source/User_Guide/newtonian_cooling.txt | 40 +++++++++++++++++++++ doc/source/User_Guide/under_development.rst | 4 +-- src/Physics/Controls.F90 | 4 +-- src/Physics/Sphere_Physical_Space.F90 | 6 ++-- 4 files changed, 47 insertions(+), 7 deletions(-) create mode 100644 doc/source/User_Guide/newtonian_cooling.txt diff --git a/doc/source/User_Guide/newtonian_cooling.txt b/doc/source/User_Guide/newtonian_cooling.txt new file mode 100644 index 00000000..88dbc143 --- /dev/null +++ b/doc/source/User_Guide/newtonian_cooling.txt @@ -0,0 +1,40 @@ +.. _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). + +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 + / + + diff --git a/doc/source/User_Guide/under_development.rst b/doc/source/User_Guide/under_development.rst index 8c4e31d4..b34f8d33 100644 --- a/doc/source/User_Guide/under_development.rst +++ b/doc/source/User_Guide/under_development.rst @@ -16,8 +16,6 @@ Arbitrary Scalar Fields .. include:: arbitrary_scalar_fields.txt - - - +.. include:: newtonian_cooling.txt diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index b8f4cb9c..1de01030 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -80,14 +80,14 @@ Module Controls Logical :: newtonian_cooling = .false. ! Turn newtonian_cooling on/off Integer :: newtonian_cooling_type = 1 Real*8 :: newtonian_cooling_time = 1.0d22 - Real*8 :: newtonian_delta_tvar_Eq = 0.0d0 + 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, & & newtonian_cooling, newtonian_cooling_type, newtonian_cooling_time, & - & newtonian_delta_tvar_eq + & newtonian_cooling_tvar_amp !/////////////////////////////////////////////////////////////////////////// ! Temporal Controls diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index dc5e3ae6..b6b9b9e9 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -53,10 +53,11 @@ Subroutine Physical_Space_Init() 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_delta_tvar_eq + tvar_eq(k,r,t) = newtonian_cooling_tvar_amp Enddo Enddo Enddo @@ -65,10 +66,11 @@ Subroutine Physical_Space_Init() 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 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_delta_tvar_eq*costheta(t)*sinphi(k) + tvar_eq(k,r,t) = newtonian_cooling_tvar_amp*costheta(t)*sinphi(k) Enddo Enddo Enddo From 1277d0bcd5a9d740dd098fb171291dd0f827fab1 Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Fri, 21 Jun 2024 10:41:05 -0600 Subject: [PATCH 07/10] Added radial variation to cooling_type =2 --- doc/source/User_Guide/newtonian_cooling.txt | 4 +++- src/Physics/Sphere_Physical_Space.F90 | 7 +++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/doc/source/User_Guide/newtonian_cooling.txt b/doc/source/User_Guide/newtonian_cooling.txt index 88dbc143..f25d1567 100644 --- a/doc/source/User_Guide/newtonian_cooling.txt +++ b/doc/source/User_Guide/newtonian_cooling.txt @@ -24,7 +24,9 @@ and a value of 2 yields .. math:: :label: ncteq2 - \Delta\Theta_\mathrm{eq} = A\mathrm{cos}(\theta)\mathrm{sin}(\phi). + \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``. diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index b6b9b9e9..bc895049 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -44,6 +44,9 @@ Module Sphere_Physical_Space 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. @@ -66,11 +69,15 @@ Subroutine Physical_Space_Init() 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*costheta(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 From fba424b42147c9fdc203b59cbf39b4b97eac8332 Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Fri, 21 Jun 2024 11:47:05 -0600 Subject: [PATCH 08/10] Changed to ell=1, m=1 --- src/Physics/Sphere_Physical_Space.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index bc895049..0e85c57f 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -76,7 +76,7 @@ Subroutine Physical_Space_Init() 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*costheta(t)*sinphi(k) + 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 From 497336cab3851cc58909ad75d433ca932476eb0a Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Fri, 21 Jun 2024 17:34:14 -0600 Subject: [PATCH 09/10] Implemented ability to read in custom cooling profile and updated documentation. --- doc/source/User_Guide/newtonian_cooling.txt | 9 +++---- src/Physics/Controls.F90 | 4 ++- src/Physics/Initial_Conditions.F90 | 7 ++++++ src/Physics/Sphere_Physical_Space.F90 | 27 ++++++++++++++------- 4 files changed, 32 insertions(+), 15 deletions(-) diff --git a/doc/source/User_Guide/newtonian_cooling.txt b/doc/source/User_Guide/newtonian_cooling.txt index f25d1567..04c8ce19 100644 --- a/doc/source/User_Guide/newtonian_cooling.txt +++ b/doc/source/User_Guide/newtonian_cooling.txt @@ -17,18 +17,16 @@ At present, :math:`\Delta\Theta_\mathrm{eq}` is allowed to take one of two forms .. math:: :label: ncteq1 - \Delta\Theta_\mathrm{eq} = A, + \Delta\Theta_\mathrm{eq} = A\,f_\mathrm{c}(r), 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}), + \Delta\Theta_\mathrm{eq} = A\,f_\mathrm{c}(r)\mathrm{sin}(\theta)\mathrm{sin}(\phi), -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``. +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``. :: @@ -37,6 +35,7 @@ The amplitude :math:`A` is controlled by setting the variable ``newtonian_cooli 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/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 1de01030..9ba23c53 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -81,13 +81,15 @@ Module Controls 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_tvar_amp, newtonian_cooling_profile_file !/////////////////////////////////////////////////////////////////////////// ! Temporal Controls 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_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index 0e85c57f..274bc523 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -44,9 +44,18 @@ Module Sphere_Physical_Space Subroutine Physical_Space_Init() Implicit None Integer :: k, r, t - Real*8, Allocatable :: press(:) - - + 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. @@ -56,11 +65,11 @@ Subroutine Physical_Space_Init() If (newtonian_cooling_type .eq. 1) Then ! No angular variation - If (my_rank .eq. 0) Write(6,*) 'Newtonian cooling is active. Type = 1' + 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 + tvar_eq(k,r,t) = newtonian_cooling_tvar_amp*newtonian_cooling_profile(r) Enddo Enddo Enddo @@ -70,20 +79,20 @@ Subroutine Physical_Space_Init() 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' + 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)*log(press(r)/press(N_R))/log(press(1)/press(N_R)) + 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() From 8cbb42fbc45d5aecbd9af9271a74adaefd3d397f Mon Sep 17 00:00:00 2001 From: Nick Featherstone Date: Mon, 15 Jul 2024 18:32:57 -0600 Subject: [PATCH 10/10] Update Controls.F90 Added comma and line continuation character that were missing following initial attempt at resolving merge conflict. --- src/Physics/Controls.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index dbcba2eb..9e9d767e 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -90,7 +90,7 @@ Module Controls & 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 + & newtonian_cooling_tvar_amp, newtonian_cooling_profile_file, & & pseudo_incompressible !///////////////////////////////////////////////////////////////////////////