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.

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

Set f_14 depending on convention; don't reset c's and f's for diffusion types = 3 #505

Merged
merged 3 commits into from
May 30, 2024
Merged
Changes from all commits
Commits
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
94 changes: 68 additions & 26 deletions src/Physics/PDE_Coefficients.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1063,6 +1063,9 @@ End Subroutine Initialize_Reference_Heating
Subroutine Augment_Reference()
Implicit None
Real*8, Allocatable :: temp_functions(:,:), temp_constants(:)
Character(len=2) :: intstring
Character*12 :: dstring
Character*8 :: dofmt = '(ES12.5)'

If (my_rank .eq. 0) Then
Call stdout%print('Reference state will be augmented.')
Expand Down Expand Up @@ -1111,10 +1114,29 @@ Subroutine Augment_Reference()
temp_constants(2) = ra_constants(2)
Endif

If (use_custom_function(14)) Then
If (with_custom_reference .and. use_custom_function(14)) Then
! Set c_11 = 1 by default
If (.not. use_custom_constant(11)) Then
If (my_rank .eq. 0) Then
Call stdout%print('User didn''t set c_11; now setting c_11 to 1.')
Endif
ra_constants(11) = 1.0d0
Else
If (my_rank .eq. 0) Then
Write(dstring,dofmt) ra_constants(11)
Call stdout%print('User set c_11 to '//Trim(dstring))
Endif
Endif
If (my_rank .eq. 0) Then
Call stdout%print('Background thermal gradient is set to:')
Call stdout%print('c_11*f_14')
! Make a warning if the user is trying to change a dsdr that physically was assumed = 0
If (any(reference_type .eq. (/1, 2/))) Then
Write(intstring, '(I2)') reference_type
Call stdout%print('WARNING: reference_type = '//Adjustl(intstring)//' assumes dsdr = 0 by default.')
Call stdout%print('make sure the combination of')
Call stdout%print('reference_type = '//Adjustl(intstring)//' and nonzero dsdr was intended.')
Endif
Call stdout%print(' ')
Endif
ref%dsdr(:) = ra_constants(11)*ra_functions(:,14)
Expand Down Expand Up @@ -2013,22 +2035,31 @@ 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

! the combination c_11 and f_14 cannot be backed out from the "ref" structure
! (which only has ref%dsdr = c_11*f_14)
! and will depend on the convention of the different reference_types
If (.not. use_custom_function(14)) Then
! If user set dsdr via "with custom", then c_11 and f_14 have already been set
Select Case(reference_type)
Case(1,2)
ra_constants(11) = 0.0d0
ra_functions(:,14) = 0.0d0
Case(3)
ra_constants(11) = 1.0d0
ra_functions(:,14) = ref%dsdr
! For reference_type = 4, c_11 and f_14 have been set already
Case(5)
ra_constants(11) = Prandtl_Number*Buoyancy_Number_Visc/Rayleigh_Number
ra_functions(:,14) = ref%dsdr/ra_constants(11)
End Select
Endif

ra_functions(:,1) = ref%density
ra_functions(:,4) = ref%temperature
ra_functions(:,8) = ref%dlnrho
ra_functions(:,9) = ref%d2lnrho
ra_functions(:,10) = ref%dlnT
ra_functions(:,14) = ref%dsdr

End Subroutine Set_Reference_Equation_Coefficients

Expand Down Expand Up @@ -2077,30 +2108,41 @@ Subroutine Set_Diffusivity_Equation_Coefficients
Endif
Endif

ra_constants(5) = nu_norm
ra_functions(:,3) = nu(:)/nu_norm
ra_functions(:,11) = dlnu(:)
If (nu_type .ne. 3) Then
! if the type is 3, the constant and function were set directly
ra_constants(5) = nu_norm
ra_functions(:,3) = nu(:)/nu_norm
ra_functions(:,11) = dlnu(:)
Endif

ra_constants(6) = kappa_norm
ra_functions(:,5) = kappa(:)/kappa_norm
ra_functions(:,12) = dlnkappa(:)

If (kappa_type .ne. 3) Then
ra_constants(6) = kappa_norm
ra_functions(:,5) = kappa(:)/kappa_norm
ra_functions(:,12) = dlnkappa(:)
Endif

If (magnetism) Then
ra_constants(7) = eta_norm
ra_functions(:,7) = eta(:)/eta_norm
ra_functions(:,13) = dlneta(:)
If (eta_type .ne. 3) Then
ra_constants(7) = eta_norm
ra_functions(:,7) = eta(:)/eta_norm
ra_functions(:,13) = dlneta(:)
Endif
Endif ! if no magnetism, all of the above are already zero

Do i = 1, n_active_scalars
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,:)
If (kappa_chi_a_type(i) .ne. 3) Then
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,:)
Endif
Enddo

Do i = 1, n_passive_scalars
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,:)
If (kappa_chi_p_type(i) .ne. 3) Then
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,:)
Endif
Enddo

End Subroutine Set_Diffusivity_Equation_Coefficients
Expand Down
Loading