Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ref%dsdr should now be scalable by constant c_11 #489

Merged
merged 4 commits into from
Dec 5, 2023
Merged

ref%dsdr should now be scalable by constant c_11 #489

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4c76a13
ref%dsdr should now be scalable by constant c_11
Nov 25, 2023
ce25463
forgot about with_custom stuff
Nov 25, 2023
d3cec9f
I had accidentally removed reading of fset
Nov 26, 2023
9fe898a
added c_11 to documentation
illorenzo7 Dec 2, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 8 additions & 5 deletions doc/source/User_Guide/physics_math_overview.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,7 @@ respectively. The evolution of :math:`\Theta` is described by
.. math::
:label: theta_evol

\mathrm{f}_1(r)\,\mathrm{f}_4(r)\left[\frac{\partial \Theta}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta + \mathrm{f}_{14}(r)v_r\right] =\
\mathrm{f}_1(r)\,\mathrm{f}_4(r)\left[\frac{\partial \Theta}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta + c_{11}\,\mathrm{f}_{14}(r)v_r\right] =\
c_6\,\boldsymbol{\nabla}\cdot\left[\mathrm{f}_1(r)\,\mathrm{f}_4(r)\,\mathrm{f}_5(r)\,\boldsymbol{\nabla}\Theta \right] \\
+\ c_{10}\,\mathrm{f}_6(r)
+ c_8\,\Phi(r,\theta,\phi)
Expand Down Expand Up @@ -235,7 +235,8 @@ assigning the following to the functions :math:`\mathrm{f}_i(r)` and the constan
\mathrm{f}_7(r) &\rightarrow \tilde{\eta}(r)\; &c_7 &\rightarrow \frac{1}{Pm} \\
&\vdots &c_8&\rightarrow 0\\
&\vdots &c_9&\rightarrow 0 \\
\mathrm{f}_{14}(r)&\rightarrow 0\; &c_{10}&\rightarrow 0.\end{aligned}
&\vdots &c_{10}&\rightarrow 0 \\
\mathrm{f}_{14}(r)&\rightarrow 0\; &c_{11}&\rightarrow 0.\end{aligned}

Here the tildes denote nondimensional radial profiles, e.g., :math:`\tilde{\nu}(r) = \nu(r)/\nu_o`.

Expand Down Expand Up @@ -296,7 +297,8 @@ When run in dimensional, anelastic mode (cgs units; **reference_type=2**
\mathrm{f}_7(r) &\rightarrow \eta(r)\; &c_7 &\rightarrow 1 \\
&\vdots &c_8&\rightarrow 1\\
&\vdots &c_9&\rightarrow \frac{1}{4\pi} \\
\mathrm{f}_{14}(r)&\rightarrow \frac{d\hat{S}}{dr }&c_{10}&\rightarrow L_*.\end{aligned}
&\vdots &c_{10}&\rightarrow L_* \\
\mathrm{f}_{14}(r)&\rightarrow \frac{d\hat{S}}{dr }&c_{11}&\rightarrow 1.\end{aligned}

Here :math:`\hat{\rho}(r)`, :math:`\hat{T}(r)`, and :math:`d\hat{S}/dr` are the spherically symmetric, time-independent reference-state
density, temperature, and entropy gradient, respectively. The thermal variables satisfy the
Expand Down Expand Up @@ -392,7 +394,8 @@ choices result in the functions :math:`\mathrm{f}_i` and the constants
\mathrm{f}_7(r) &\rightarrow \tilde{\eta}(r) \; &c_7 &\rightarrow \frac{\mathrm{E}}{\mathrm{Pm}} \\
&\vdots &c_8&\rightarrow \frac{\mathrm{E}\,\mathrm{Di}}{\mathrm{Ra}^*}\\
&\vdots &c_9&\rightarrow \frac{\mathrm{E}^2\,\mathrm{Di}}{\mathrm{Pm}^2\mathrm{Ra}^*}\\
\mathrm{f}_{14}(r)&\rightarrow \frac{d\tilde{S}}{dr }&c_{10}&\rightarrow L_*.\end{aligned}
&\vdots &c_{10}&\rightarrow L_* \\
\mathrm{f}_{14}(r)&\rightarrow 0&c_{11}&\rightarrow 0.\end{aligned}

As in the Boussinesq case, the nondimensional diffusivities are defined according to, e.g., :math:`\tilde{\nu}(r) \equiv \nu(r)/\nu_o`. The nondimensional heating :math:`\tilde{Q}(r)` is defined such that its volume integral equals the nondimensional **luminosity** or **heating_integral** set in the *main_input* file. As in the dimensional anelastic case, the volume integral of :math:`\mathrm{f}_6(r)` equals unity, and :math:`\mathrm{c}_{10} = L_*`. The unit for luminosity in this nondimensionalization (to get a dimensional luminosity from the nondimensional :math:`L_*`) is :math:`\rho_oL^3T_o\Delta s\Omega_0`.

Expand Down Expand Up @@ -428,7 +431,7 @@ We thus arrive at the following nondimensionalized equations:
+ \mathrm{E}\boldsymbol{\nabla}\cdot\boldsymbol{\mathcal{D}}\\
%
%
\tilde{\rho}(r)\,\tilde{T}(r)\left[\frac{\partial \Theta}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta + v_r\frac{d\hat{S}}{dr}\right] =\;
\tilde{\rho}(r)\,\tilde{T}(r)\left[\frac{\partial \Theta}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla}\Theta\right] =\;
&\frac{\mathrm{E}}{\mathrm{Pr}}\boldsymbol{\nabla}\cdot\left[\tilde{\kappa}(r)\tilde{\rho}(r)\,\tilde{T}(r)\,\boldsymbol{\nabla}\Theta \right] % diffusion
+ \tilde{Q}(r) % Internal heating
\\
Expand Down
51 changes: 36 additions & 15 deletions src/Physics/PDE_Coefficients.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ Module PDE_Coefficients
Real*8 :: Angular_Velocity = -1.0d0 ! Frame rotation rate (sets Coriolis force)

! Custom reference-state variables (reference_type = 4)
Integer, Parameter :: max_ra_constants = 10 + 2*n_scalar_max
Integer, Parameter :: max_ra_constants = 11 + 2*n_scalar_max
Integer, Parameter :: max_ra_functions = 14 + 2*n_scalar_max
Integer :: n_ra_constants
Integer :: n_ra_functions
Expand Down Expand Up @@ -254,7 +254,7 @@ End Subroutine Initialize_Reference
Subroutine Allocate_Reference_State
Implicit None

n_ra_constants = 10 + 2*(n_active_scalars + n_passive_scalars)
n_ra_constants = 11 + 2*(n_active_scalars + n_passive_scalars)
n_ra_functions = 14 + 2*(n_active_scalars + n_passive_scalars)

Allocate(ref%density(1:N_R))
Expand Down Expand Up @@ -1037,7 +1037,7 @@ Subroutine Augment_Reference()
Call stdout%print('Only heating, buoyancy, or background dTdr/dSdr may be modified.')
Call stdout%print('Heating requires both c_10 and f_6 to be set.')
Call stdout%print('Buoyancy requires both c_2 and f_2 to be set.')
Call stdout%print('dTdr/dSdr requires f_14 to be set.')
Call stdout%print('dTdr/dSdr requires f_14 to be set (sets c_11 to 1 if unspecified).')
Call stdout%print('Reading from: '//Trim(custom_reference_file))
Endif

Expand Down Expand Up @@ -1080,10 +1080,11 @@ Subroutine Augment_Reference()
If (use_custom_function(14)) Then
If (my_rank .eq. 0) Then
Call stdout%print('Background thermal gradient is set to:')
Call stdout%print('f_14')
Call stdout%print('c_11*f_14')
Call stdout%print(' ')
Endif
ref%dsdr(:) = ra_functions(:,14)
ref%dsdr(:) = ra_constants(11)*ra_functions(:,14)
temp_constants(11) = ra_constants(11)
temp_functions(:,14) = ra_functions(:,14)
Endif

Expand Down Expand Up @@ -1196,7 +1197,7 @@ Subroutine Get_Custom_Reference()
ref%Lorentz_Coeff = ra_constants(4)
ref%ohmic_amp(:) = ra_constants(9)/(ref%density(:)*ref%temperature(:))

ref%dsdr(:) = ra_functions(:,14)
ref%dsdr(:) = ra_constants(11)*ra_functions(:,14)

End Subroutine Get_Custom_Reference

Expand Down Expand Up @@ -1277,9 +1278,20 @@ Subroutine Read_Custom_Reference_File(filename)

! Read in constants and their 'set' flags
Read(15) eqversion
Read(15) cset(1:n_ra_constants)
Read(15) fset(1:n_ra_functions)
Read(15) input_constants(1:n_ra_constants)
If (eqversion .eq. 1) Then
!Read(15) cset(1:n_ra_constants-1) ! c_11 didn't exist yet
Read(15) cset(1:10) ! equation_coefficients couldn't write the custom active/passive scalar constants yet, and c_11 didn't exist yet
Read(15) fset(1:n_ra_functions)
Read(15) input_constants(1:10)
cset(11) = 1 ! treat this as if c_11 = 1 was specified in the custom reference file
input_constants(11) = 1.0d0
If (my_rank .eq. 0) Call stdout%print('got here')
Else
Read(15) cset(1:n_ra_constants)
Read(15) fset(1:n_ra_functions)
Read(15) input_constants(1:n_ra_constants)
Call stdout%print('actualy got here')
Endif

! Cset(i) is 1 if a constant(i) was set; it is 0 otherwise.
! The logic below deals with a constant set in both the reference
Expand Down Expand Up @@ -1680,12 +1692,12 @@ Subroutine Initialize_Transport_Coefficients()
do i = 1, n_active_scalars
Call Initialize_Diffusivity(kappa_chi_a(i,:),dlnkappa_chi_a(i,:),&
kappa_chi_a_top(i),kappa_chi_a_type(i),kappa_chi_a_power(i),&
11+(i-1)*2,15+(i-1)*2,16+(i-1)*2)
12+(i-1)*2,15+(i-1)*2,16+(i-1)*2)
end do
do i = 1, n_passive_scalars
Call Initialize_Diffusivity(kappa_chi_p(i,:),dlnkappa_chi_p(i,:),&
kappa_chi_p_top(i),kappa_chi_p_type(i),kappa_chi_p_power(i),&
11+(n_active_scalars+i-1)*2,15+(n_active_scalars+i-1)*2,16+(n_active_scalars+i-1)*2)
12+(n_active_scalars+i-1)*2,15+(n_active_scalars+i-1)*2,16+(n_active_scalars+i-1)*2)
end do

If (viscous_heating) Then
Expand Down Expand Up @@ -1723,15 +1735,15 @@ Subroutine Initialize_Transport_Coefficients()
If (kappa_chi_a_type(i+1) .eq. 3) Then
temp_functions(:,15+i*2) = ra_functions(:,15+i*2)
temp_functions(:,16+i*2) = ra_functions(:,16+i*2)
temp_constants(11+i*2) = ra_constants(11+i*2)
temp_constants(12+i*2) = ra_constants(12+i*2)
Endif
end do

do i = 0, n_passive_scalars-1
If (kappa_chi_p_type(i+1) .eq. 3) Then
temp_functions(:,15+(n_active_scalars+i)*2) = ra_functions(:,15+(n_active_scalars+i)*2)
temp_functions(:,16+(n_active_scalars+i)*2) = ra_functions(:,16+(n_active_scalars+i)*2)
temp_constants(11+(n_active_scalars+i)*2) = ra_constants(11+(n_active_scalars+i)*2)
temp_constants(12+(n_active_scalars+i)*2) = ra_constants(12+(n_active_scalars+i)*2)
Endif
end do

Expand Down Expand Up @@ -1964,6 +1976,15 @@ Subroutine Set_Reference_Equation_Coefficients
ra_constants(4) = ref%Lorentz_Coeff
ra_constants(8) = ref%viscous_amp(1)*ref%temperature(1)/2.0d0
ra_constants(9) = ref%ohmic_amp(1)*ref%density(1)*ref%temperature(1)
Select Case(reference_type)
Case(1,2)
ra_constants(11) = 0.0d0
Case(3)
ra_constants(11) = 1.0d0
Case(5)
ra_constants(11) = Prandtl_Number*Buoyancy_Number_Visc/Rayleigh_Number
End Select


ra_functions(:,1) = ref%density
ra_functions(:,4) = ref%temperature
Expand Down Expand Up @@ -2034,13 +2055,13 @@ Subroutine Set_Diffusivity_Equation_Coefficients
Endif ! if no magnetism, all of the above are already zero

Do i = 1, n_active_scalars
ra_constants(11+(i-1)*2) = kappa_chi_a_norm(i)
ra_constants(12+(i-1)*2) = kappa_chi_a_norm(i)
ra_functions(:,15+(i-1)*2) = kappa_chi_a(i,:)/kappa_chi_a_norm(i)
ra_functions(:,16+(i-1)*2) = dlnkappa_chi_a(i,:)
Enddo

Do i = 1, n_passive_scalars
ra_constants(11+(n_active_scalars+i-1)*2) = kappa_chi_p_norm(i)
ra_constants(12+(n_active_scalars+i-1)*2) = kappa_chi_p_norm(i)
ra_functions(:,15+(n_active_scalars+i-1)*2) = kappa_chi_p(i,:)/kappa_chi_p_norm(i)
ra_functions(:,16+(n_active_scalars+i-1)*2) = dlnkappa_chi_p(i,:)
Enddo
Expand Down
Loading