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

Fix bugs for custom diffusions #429

Merged
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
41 changes: 33 additions & 8 deletions src/Physics/PDE_Coefficients.F90
Original file line number Diff line number Diff line change
Expand Up @@ -828,9 +828,9 @@ Subroutine Get_Custom_Reference()
ref%dlnrho(:) = ra_functions(:,8)
ref%d2lnrho(:) = ra_functions(:,9)
ref%buoyancy_coeff(:) = ra_constants(2)*ra_functions(:,2)
do i = 0, n_active_scalars-1
ref%chi_buoyancy_coeff(i,:) = ra_constants(12+i*2)*ra_functions(:,2)
end do
Do i = 1, n_active_scalars
ref%chi_buoyancy_coeff(i,:) = ra_constants(12+(i-1)*2)*ra_functions(:,2)
Enddo

ref%temperature(:) = ra_functions(:,4)
ref%dlnT(:) = ra_functions(:,10)
Expand Down Expand Up @@ -1280,22 +1280,47 @@ Subroutine Initialize_Transport_Coefficients()
Implicit None
Integer :: i
Real*8, Allocatable :: temp_functions(:,:), temp_constants(:)
Logical :: restore
Logical :: restore, need_custom

restore = .false.

Call Allocate_Transport_Coefficients

If ((.not. custom_reference_read) .and. &
((nu_type .eq. 3) .or. (kappa_type .eq. 3) .or. (eta_type .eq. 3))) Then
! Figure out if we need to read anything from the custom file
! (many "types" to check now because of the new scalar diffusion coefficients)
need_custom = .false.
If ( (nu_type .eq. 3) .or. (kappa_type .eq. 3) .or. (eta_type .eq. 3) ) Then
need_custom = .true.
EndIf

Do i = 1, n_active_scalars
If (kappa_chi_a_type(i) .eq. 3) Then
need_custom = .true.
Endif
Enddo

Do i = 1, n_passive_scalars
If (kappa_chi_p_type(i) .eq. 3) Then
need_custom = .true.
Endif
Enddo

If ((.not. custom_reference_read) .and. need_custom) Then
Allocate(temp_functions(1:n_r, 1:n_ra_functions))
Allocate(temp_constants(1:n_ra_constants))
temp_functions(:,:) = ra_functions(:,:)
temp_constants(:) = ra_constants(:)
! Note that ra_constants is allocated up to max_ra_constants,
! which could be more than n_ra_constants
temp_constants(:) = ra_constants(1:n_ra_constants)
restore = .true.
! If we read the custom file, we may overwrite things besides the diffusion coefficients
! We "back up" the current reference state in temp_constants and temp_functions
! Below, we modify only the "temp" equation coefficients associated with custom diffusions
! Then we restore ra_constants and ra_functions from the "temp" arrays
Call Read_Custom_Reference_File(custom_reference_file)
EndIf


Call Initialize_Diffusivity(nu,dlnu,nu_top,nu_type,nu_power,5,3,11)
Call Initialize_Diffusivity(kappa,dlnkappa,kappa_top,kappa_type,kappa_power,6,5,12)
do i = 1, n_active_scalars
Expand Down Expand Up @@ -1362,7 +1387,7 @@ Subroutine Initialize_Transport_Coefficients()
temp_constants(5) = ra_constants(5)
Endif

ra_constants(:) = temp_constants(:)
ra_constants(1:n_ra_constants) = temp_constants(:)
ra_functions(:,:) = temp_functions(:,:)
DeAllocate(temp_functions, temp_constants)

Expand Down