From 7559b3f8b89cd33f1353e10d14ecd83a1112f8b4 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Tue, 13 Jun 2023 12:20:49 -0600 Subject: [PATCH 01/10] An initial implementation of coupled bcs. Untested but compiles. Think this can actually be reimplemented in a considerably simpler way but I would like to add a test to this version first. This may break current tests as it changes the default boundary condition logic (it is no longer possible to just set the bc type based on what it is not because there are more than two types of bc now). --- src/Physics/BoundaryConditions.F90 | 89 +++++++++++- src/Physics/Sphere_Linear_Terms.F90 | 205 +++++++++++++++++++++++++++- 2 files changed, 285 insertions(+), 9 deletions(-) diff --git a/src/Physics/BoundaryConditions.F90 b/src/Physics/BoundaryConditions.F90 index d242ac9c..0c80fac0 100755 --- a/src/Physics/BoundaryConditions.F90 +++ b/src/Physics/BoundaryConditions.F90 @@ -38,11 +38,21 @@ Module BoundaryConditions Logical :: Fix_dTdr_Top = .False. Logical :: Fix_dTdr_Bottom = .False. + Logical :: Couple_Tvar_Top = .False. + Logical :: Couple_Tvar_Bottom = .False. + Logical :: Couple_dTdr_Top = .False. + Logical :: Couple_dTdr_Bottom = .False. + Logical :: Fix_chivar_a_Top(1:n_scalar_max) = .True. Logical :: Fix_chivar_a_Bottom(1:n_scalar_max) = .True. Logical :: Fix_dchidr_a_Top(1:n_scalar_max) = .False. Logical :: Fix_dchidr_a_Bottom(1:n_scalar_max) = .False. + Logical :: Couple_chivar_a_Top(1:n_scalar_max) = .False. + Logical :: Couple_chivar_a_Bottom(1:n_scalar_max) = .False. + Logical :: Couple_dchidr_a_Top(1:n_scalar_max) = .False. + Logical :: Couple_dchidr_a_Bottom(1:n_scalar_max) = .False. + Logical :: Fix_chivar_p_Top(1:n_scalar_max) = .True. Logical :: Fix_chivar_p_Bottom(1:n_scalar_max) = .True. Logical :: Fix_dchidr_p_Top(1:n_scalar_max) = .False. @@ -63,11 +73,41 @@ Module BoundaryConditions Real*8 :: dTdr_Top = 0.0d0 Real*8 :: dTdr_Bottom = 0.0d0 + Real*8 :: T_chi_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dchidr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dTdr_coeff_top = 0.0d0 + Real*8 :: T_chi_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dchidr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dTdr_coeff_bottom = 0.0d0 + Real*8 :: dTdr_chi_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_dchidr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_T_coeff_top = 0.0d0 + Real*8 :: dTdr_chi_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_dchidr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_T_coeff_bottom = 0.0d0 + Real*8 :: chi_a_Bottom(1:n_scalar_max) = 1.0d0 Real*8 :: chi_a_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_a_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_a_Bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_chi_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_dchidr_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_T_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_chi_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_dchidr_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_T_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_chi_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_dchidr_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_T_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_chi_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_dchidr_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_T_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_p_Bottom(1:n_scalar_max) = 1.0d0 Real*8 :: chi_p_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_p_Top(1:n_scalar_max) = 0.0d0 @@ -120,7 +160,17 @@ Module BoundaryConditions fix_dchidr_a_bottom, fix_dchidr_a_top, dchidr_a_top, dchidr_a_bottom, & chi_p_top_file, chi_p_bottom_file, dchidr_p_top_file, dchidr_p_bottom_file, & fix_chivar_p_top, fix_chivar_p_bottom, chi_p_bottom, chi_p_top, & - fix_dchidr_p_bottom, fix_dchidr_p_top, dchidr_p_top, dchidr_p_bottom + fix_dchidr_p_bottom, fix_dchidr_p_top, dchidr_p_top, dchidr_p_bottom, & + couple_Tvar_top, couple_Tvar_bottom, couple_dtdr_top, couple_dtdr_bottom, & + T_chi_coeff_top, T_dchidr_coeff_top, T_dTdr_coeff_top, & + T_chi_coeff_bottom, T_dchidr_coeff_bottom, T_dTdr_coeff_bottom, & + dTdr_chi_coeff_top, dTdr_dchidr_coeff_top, dTdr_T_coeff_top, & + dTdr_chi_coeff_bottom, dTdr_dchidr_coeff_bottom, dTdr_T_coeff_bottom, & + couple_chivar_a_top, couple_chivar_a_bottom, couple_dchidr_a_top, couple_dchidr_a_bottom, & + chi_chi_coeff_top, chi_dchidr_coeff_top, chi_T_coeff_top, chi_dTdr_coeff_top, & + chi_chi_coeff_bottom, chi_dchidr_coeff_bottom, chi_T_coeff_bottom, chi_dTdr_coeff_bottom, & + dchidr_chi_coeff_top, dchidr_dchidr_coeff_top, dchidr_T_coeff_top, dchidr_dTdr_coeff_top, & + dchidr_chi_coeff_bottom, dchidr_dchidr_coeff_bottom, dchidr_T_coeff_bottom, dchidr_dTdr_coeff_bottom Contains @@ -129,18 +179,43 @@ Subroutine Initialize_Boundary_Conditions() Real*8 :: tilt_angle_radians,a,b Real*8 :: Ftop, rhotk_top, Fbottom Integer :: i + Integer :: n_active_bcs + + n_active_bcs = count( (/ fix_tvar_top, fix_dtdr_top, couple_tvar_top, couple_dtdr_top /) ) + if (n_active_bcs /= 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for tvar on top boundary.") + endif - fix_tvar_top = .not. fix_dtdr_top - fix_tvar_bottom = .not. fix_dtdr_bottom + n_active_bcs = count( (/ fix_tvar_bottom, fix_dtdr_bottom, couple_tvar_bottom, couple_dtdr_bottom /) ) + if (n_active_bcs /= 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for tvar on bottom boundary.") + endif do i = 1, n_active_scalars - fix_chivar_a_top(i) = .not. fix_dchidr_a_top(i) - fix_chivar_a_bottom(i) = .not. fix_dchidr_a_bottom(i) + n_active_bcs = count( (/ fix_chivar_a_top(i), fix_dchidr_a_top(i), couple_chivar_a_top(i), couple_dchidr_a_top(i) /) ) + if (n_active_bcs /= 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_a on top boundary.") + endif + + n_active_bcs = count( (/ fix_chivar_a_bottom(i), fix_dchidr_a_bottom(i), & + couple_chivar_a_bottom(i), couple_dchidr_a_bottom(i) /) ) + if (n_active_bcs /= 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_a on bottom boundary.") + endif end do do i = 1, n_passive_scalars - fix_chivar_p_top(i) = .not. fix_dchidr_p_top(i) - fix_chivar_p_bottom(i) = .not. fix_dchidr_p_bottom(i) + n_active_bcs = count( (/ fix_chivar_p_top(i), fix_dchidr_p_top(i) /) ) + if (n_active_bcs /= 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_p on top boundary.") + call stdout%print(" ") + endif + + n_active_bcs = count( (/ fix_chivar_p_bottom(i), fix_dchidr_p_bottom(i) /) ) + if (n_active_bcs /= 1) then + call stdout%print(" -- Error: Incompatible boundary conditions for chivar_p on bottom boundary.") + call stdout%print(" ") + endif end do If (no_slip_top .and. strict_L_Conservation) Then diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90 index f8b37e7c..63146d3b 100755 --- a/src/Physics/Sphere_Linear_Terms.F90 +++ b/src/Physics/Sphere_Linear_Terms.F90 @@ -521,7 +521,7 @@ Subroutine Set_Boundary_Conditions(mode_ind) Implicit None Real*8 :: samp,one Integer, Intent(In) :: mode_ind - Integer :: l, r,lp, i + Integer :: l, r,lp, i, j one = 1.0d0 lp = mode_ind @@ -558,6 +558,30 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_top) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_top) Then + do i = 1, n_active_scalars + samp = -T_chi_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_top) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif r = N_R If (fix_tvar_bottom) Then @@ -566,8 +590,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_bottom) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_bottom) Then + do i = 1, n_active_scalars + samp = -T_chi_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_bottom) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif - ! PASSIVE + ! active scalars do i = 1, n_active_scalars r = 1 If (fix_chivar_a_top(i)) Then @@ -576,6 +624,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_top(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_top(i)) Then + do j = 1, n_active_scalars + samp = -chi_chi_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_dchidr_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_top(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_chi_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_dchidr_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif r = N_R If (fix_chivar_a_bottom(i)) Then @@ -584,8 +658,35 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_bottom(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -chi_chi_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_dchidr_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_chi_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_dchidr_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif end do + ! passive scalars do i = 1, n_passive_scalars r = 1 If (fix_chivar_p_top(i)) Then @@ -663,6 +764,30 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_top) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_top) Then + do i = 1, n_active_scalars + samp = -T_chi_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_top) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_coeff_top(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_top + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif r = N_R If (fix_tvar_bottom) Then @@ -671,6 +796,30 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dtdr_bottom) Then Call Load_BC(lp,r,teq,tvar,one,1) Endif + If (couple_tvar_bottom) Then + do i = 1, n_active_scalars + samp = -T_chi_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -T_dchidr_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = -T_dTdr_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif + If (couple_dtdr_bottom) Then + do i = 1, n_active_scalars + samp = -dTdr_chi_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,0) + samp = -dTdr_dchidr_coeff_bottom(i) + Call Load_BC(lp,r,teq,chiavar(i),samp,1) + end do + samp = -dTdr_T_coeff_bottom + Call Load_BC(lp,r,teq,tvar,samp,0) + samp = 1.0 + Call Load_BC(lp,r,teq,tvar,samp,1) + Endif ! chivar BC do i = 1, n_active_scalars @@ -681,6 +830,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_top(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_top(i)) Then + do j = 1, n_active_scalars + samp = -chi_chi_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_dchidr_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_top(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_chi_coeff_top(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_dchidr_coeff_top(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_T_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_dTdr_coeff_top(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif r = N_R If (fix_chivar_a_bottom(i)) Then @@ -689,6 +864,32 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (fix_dchidr_a_bottom(i)) Then Call Load_BC(lp,r,chiaeq(i),chiavar(i),one,1) Endif + If (couple_chivar_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -chi_chi_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -chi_dchidr_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -chi_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -chi_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif + If (couple_dchidr_a_bottom(i)) Then + do j = 1, n_active_scalars + samp = -dchidr_chi_coeff_bottom(i,j) + if (i==j) samp = 1.0 + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) + samp = -dchidr_dchidr_coeff_bottom(i,j) + Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) + end do + samp = -dchidr_T_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) + samp = -dchidr_dTdr_coeff_bottom(i) + Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) + Endif end do do i = 1, n_passive_scalars From 398f4b1c2f01e0c60c8b256e639131056cb08c6c Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 15 Jun 2023 15:54:35 -0600 Subject: [PATCH 02/10] Extend the length of the string used to name 3D output. This allows generic scalar fields to be output which have a length of at least 5 (not just 4). --- src/IO/Spherical_IO.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/IO/Spherical_IO.F90 b/src/IO/Spherical_IO.F90 index 51bb7d88..bb52c65b 100755 --- a/src/IO/Spherical_IO.F90 +++ b/src/IO/Spherical_IO.F90 @@ -663,12 +663,12 @@ Subroutine Write_Full_3D(qty) Implicit None Real*8, Intent(In) :: qty(:,my_rmin:,my_theta_min:) Integer :: i, funit - Character*4 :: qstring + Character*5 :: qstring Character*120 :: iterstring, data_file, grid_file ! Write the data file (Parallel I/O) Write(iterstring,i_ofmt) current_iteration - Write(qstring,'(i4.4)') current_qval + Write(qstring,'(i5.5)') current_qval data_file = trim(local_file_path)//'Spherical_3D/'//trim(iterstring)//'_'//qstring Call full_3d_buffer%cache_data(qty) Call full_3d_buffer%write_data(filename=data_file) From c0620493ebc63d1e9f9c5bda83774259927bbcfe Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 15 Jun 2023 15:55:55 -0600 Subject: [PATCH 03/10] Correction to which coefficient is set to 1 for dchidr coupled bcs. --- src/Physics/Sphere_Linear_Terms.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90 index 63146d3b..00f12d85 100755 --- a/src/Physics/Sphere_Linear_Terms.F90 +++ b/src/Physics/Sphere_Linear_Terms.F90 @@ -640,9 +640,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (couple_dchidr_a_top(i)) Then do j = 1, n_active_scalars samp = -dchidr_chi_coeff_top(i,j) - if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) samp = -dchidr_dchidr_coeff_top(i,j) + if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do samp = -dchidr_T_coeff_top(i) @@ -674,9 +674,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (couple_dchidr_a_bottom(i)) Then do j = 1, n_active_scalars samp = -dchidr_chi_coeff_bottom(i,j) - if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) samp = -dchidr_dchidr_coeff_bottom(i,j) + if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do samp = -dchidr_T_coeff_bottom(i) @@ -846,9 +846,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (couple_dchidr_a_top(i)) Then do j = 1, n_active_scalars samp = -dchidr_chi_coeff_top(i,j) - if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) samp = -dchidr_dchidr_coeff_top(i,j) + if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do samp = -dchidr_T_coeff_top(i) @@ -880,9 +880,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) If (couple_dchidr_a_bottom(i)) Then do j = 1, n_active_scalars samp = -dchidr_chi_coeff_bottom(i,j) - if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) samp = -dchidr_dchidr_coeff_bottom(i,j) + if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do samp = -dchidr_T_coeff_bottom(i) From 89397853448c53ceccf585ce9077552d91848997 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 15 Jun 2023 15:56:26 -0600 Subject: [PATCH 04/10] Turn off all bcs by default and set values to 0.0 where appropriate. This may break tests that rely on the old defaults. --- src/Physics/BoundaryConditions.F90 | 34 +++++++++++++++--------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/Physics/BoundaryConditions.F90 b/src/Physics/BoundaryConditions.F90 index 0c80fac0..d2f7df2a 100755 --- a/src/Physics/BoundaryConditions.F90 +++ b/src/Physics/BoundaryConditions.F90 @@ -33,8 +33,8 @@ Module BoundaryConditions Implicit None - Logical :: Fix_Tvar_Top = .True. - Logical :: Fix_Tvar_Bottom = .True. + Logical :: Fix_Tvar_Top = .False. + Logical :: Fix_Tvar_Bottom = .False. Logical :: Fix_dTdr_Top = .False. Logical :: Fix_dTdr_Bottom = .False. @@ -43,8 +43,8 @@ Module BoundaryConditions Logical :: Couple_dTdr_Top = .False. Logical :: Couple_dTdr_Bottom = .False. - Logical :: Fix_chivar_a_Top(1:n_scalar_max) = .True. - Logical :: Fix_chivar_a_Bottom(1:n_scalar_max) = .True. + Logical :: Fix_chivar_a_Top(1:n_scalar_max) = .False. + Logical :: Fix_chivar_a_Bottom(1:n_scalar_max) = .False. Logical :: Fix_dchidr_a_Top(1:n_scalar_max) = .False. Logical :: Fix_dchidr_a_Bottom(1:n_scalar_max) = .False. @@ -53,8 +53,8 @@ Module BoundaryConditions Logical :: Couple_dchidr_a_Top(1:n_scalar_max) = .False. Logical :: Couple_dchidr_a_Bottom(1:n_scalar_max) = .False. - Logical :: Fix_chivar_p_Top(1:n_scalar_max) = .True. - Logical :: Fix_chivar_p_Bottom(1:n_scalar_max) = .True. + Logical :: Fix_chivar_p_Top(1:n_scalar_max) = .False. + Logical :: Fix_chivar_p_Bottom(1:n_scalar_max) = .False. Logical :: Fix_dchidr_p_Top(1:n_scalar_max) = .False. Logical :: Fix_dchidr_p_Bottom(1:n_scalar_max) = .False. @@ -340,7 +340,7 @@ Subroutine Generate_Boundary_Mask() If (.not. full_restart) Then - If (fix_tvar_top) Then + If (fix_tvar_top .or. couple_tvar_top) Then if (trim(T_top_file) .eq. '__nothing__') then bc_val= T_Top*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, uind) @@ -349,7 +349,7 @@ Subroutine Generate_Boundary_Mask() end if Endif - If (fix_tvar_bottom) Then + If (fix_tvar_bottom .or. couple_tvar_bottom) Then if (trim(T_bottom_file) .eq. '__nothing__') then bc_val= T_bottom*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, lind) @@ -358,7 +358,7 @@ Subroutine Generate_Boundary_Mask() end if Endif - If (fix_dtdr_top) Then + If (fix_dtdr_top .or. couple_dtdr_top) Then if (trim(dTdr_top_file) .eq. '__nothing__') then bc_val= dtdr_top*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, uind) @@ -367,7 +367,7 @@ Subroutine Generate_Boundary_Mask() end if Endif - If (fix_dtdr_bottom) Then + If (fix_dtdr_bottom .or. couple_dtdr_bottom) Then if (trim(dTdr_bottom_file) .eq. '__nothing__') then bc_val= dtdr_bottom*sqrt(four_pi) Call Set_BC(bc_val,0,0, teq,real_ind, lind) @@ -397,7 +397,7 @@ Subroutine Generate_Boundary_Mask() Endif do i = 1, n_active_scalars - If (fix_chivar_a_top(i)) Then + If (fix_chivar_a_top(i) .or. couple_chivar_a_top(i)) Then if (trim(chi_a_top_file(i)) .eq. '__nothing__') then bc_val= chi_a_Top(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, uind) @@ -419,7 +419,7 @@ Subroutine Generate_Boundary_Mask() end do do i = 1, n_active_scalars - If (fix_chivar_a_bottom(i)) Then + If (fix_chivar_a_bottom(i) .or. couple_chivar_a_bottom(i)) Then if (trim(chi_a_bottom_file(i)) .eq. '__nothing__') then bc_val= chi_a_bottom(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, lind) @@ -441,7 +441,7 @@ Subroutine Generate_Boundary_Mask() end do do i = 1, n_active_scalars - If (Fix_dchidr_a_Top(i)) Then + If (Fix_dchidr_a_Top(i) .or. couple_dchidr_a_top(i)) Then if (trim(dchidr_a_top_file(i)) .eq. '__nothing__') then bc_val= dchidr_a_top(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, uind) @@ -463,7 +463,7 @@ Subroutine Generate_Boundary_Mask() end do do i = 1, n_active_scalars - If (Fix_dchidr_a_Bottom(i)) Then + If (Fix_dchidr_a_Bottom(i) .or. couple_dchidr_a_bottom(i)) Then if (trim(dchidr_a_bottom_file(i)) .eq. '__nothing__') then bc_val= dchidr_a_bottom(i)*sqrt(four_pi) Call Set_BC(bc_val,0,0, chiaeq(i),real_ind, lind) @@ -571,8 +571,8 @@ End Subroutine Transport_Dependencies Subroutine Restore_BoundaryCondition_Defaults() Implicit None - Fix_Tvar_Top = .True. - Fix_Tvar_Bottom = .True. + Fix_Tvar_Top = .False. + Fix_Tvar_Bottom = .False. Fix_dTdr_Top = .False. Fix_dTdr_Bottom = .False. Fix_divrt_top = .False. @@ -583,7 +583,7 @@ Subroutine Restore_BoundaryCondition_Defaults() Fix_poloidalfield_bottom = .False. Impose_Dipole_Field = .False. - T_Bottom = 1.0d0 + T_Bottom = 0.0d0 T_Top = 0.0d0 dTdr_Top = 0.0d0 dTdr_Bottom = 0.0d0 From 3b8319725daa6be388e9c98e00ada3a44a4351ae Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Thu, 15 Jun 2023 15:58:13 -0600 Subject: [PATCH 05/10] Add a test for coupled_bcs. This just runs for 2 timesteps and checks that the shell slices produced match the expected values on 2 coupled fields (1 temperature, 1 chi, both coupled to another chi). --- .github/workflows/main.yml | 4 ++ tests/coupled_bcs/const/generate_input.py | 14 ++++ tests/coupled_bcs/const/main_input | 88 +++++++++++++++++++++++ tests/coupled_bcs/run_test.sh | 13 ++++ tests/coupled_bcs/test_output.py | 77 ++++++++++++++++++++ 5 files changed, 196 insertions(+) create mode 100644 tests/coupled_bcs/const/generate_input.py create mode 100644 tests/coupled_bcs/const/main_input create mode 100755 tests/coupled_bcs/run_test.sh create mode 100644 tests/coupled_bcs/test_output.py diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 5d280a88..5e5b28af 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -74,6 +74,10 @@ jobs: cd "$GITHUB_WORKSPACE" sh ./tests/custom_reference/run_test.sh + # coupled bc test + cd "$GITHUB_WORKSPACE" + sh ./tests/coupled_bcs/run_test.sh + git diff > changes.diff git diff --exit-code --name-only diff --git a/tests/coupled_bcs/const/generate_input.py b/tests/coupled_bcs/const/generate_input.py new file mode 100644 index 00000000..c85ad352 --- /dev/null +++ b/tests/coupled_bcs/const/generate_input.py @@ -0,0 +1,14 @@ +### Use rayleigh_spectral_input.py to generate generic input. + +from rayleigh_spectral_input import * + +rmin = 0.5; rmax = 1.0 + +si = SpectralInput(n_theta=1,n_r=48) +si.transform_from_rtp_function(lambda radius: 2.0 - (rmax/radius)*(radius - rmin)/(rmax - rmin), rmin=rmin, rmax=rmax) +si.write('radial_prof_init') + +si = SpectralInput() +si.add_mode(0.0+0.j, 0, 0, 0) +si.write('constant_init') + diff --git a/tests/coupled_bcs/const/main_input b/tests/coupled_bcs/const/main_input new file mode 100644 index 00000000..380f9fc7 --- /dev/null +++ b/tests/coupled_bcs/const/main_input @@ -0,0 +1,88 @@ +! This is a smaller version of the input_examples/benchmark_diagnostics_input +! model that only runs for two time-steps so that we can check the boundary conditions. + +&problemsize_namelist + n_r = 48 + n_theta = 64 + nprow = 2 + npcol = 2 + rmin = 0.5 + rmax = 1.0 +/ +&numerical_controls_namelist + chebyshev = .true. +/ +&physical_controls_namelist + rotation = .True. + magnetism = .false. + viscous_heating = .false. + ohmic_heating = .false. + n_active_scalars = 2 +/ +&temporal_controls_namelist + max_time_step = 1.0d-4 + max_iterations = 2 + alpha_implicit = 0.50001d0 + checkpoint_interval = 500000 + quicksave_interval = 1000000 + num_quicksaves = -1 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist +full3d_values = 501 507 10001 10004 10201 10204 +full3d_frequency = 2 + +shellavg_values = 501 507 10001 10004 10201 10204 +shellavg_frequency = 2 +shellavg_nrec = 1 + +shellslice_levels = 1,48 +shellslice_values = 501 507 10001 10004 10201 10204 ! T dTdr chi1 dchi1dr chi2 dchi2dr +shellslice_frequency = 2 +shellslice_nrec = 1 +/ + +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +dTdr_Top = 10.0d0 +T_Bottom = 10.0d0 +couple_dtdr_top = .true. +dTdr_chi_coeff_top(1) = 1.0 +couple_tvar_bottom = .true. +T_chi_coeff_bottom(1) = 10.0 +T_dchidr_coeff_bottom(1) = -5.0 +chi_a_Top(1) = 1.0d0 +chi_a_Bottom(1) = 2.0d0 +fix_chivar_a_top(1) = .true. +fix_chivar_a_bottom(1) = .true. +chi_a_Top(2) = 1.0d0 +dchidr_a_Bottom(2) = 20.0d0 +couple_chivar_a_top(2) = .true. +chi_chi_coeff_top(2,1) = 10.0 +chi_dchidr_coeff_top(2,1) = 5.0 +couple_dchidr_a_bottom(2) = .true. +dchidr_chi_coeff_bottom(2,1) = 1.0 +/ +&Initial_Conditions_Namelist +init_type = 8 ! File init +t_init_file = 'constant_init' +chi_a_init_file(1) = 'radial_prof_init' +chi_a_init_file(2) = 'constant_init' +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 0.0 +Prandtl_Number = 1.0d0 +Magnetic_Prandtl_Number = 5.0d0 +reference_type = 1 +heating_type = 0 ! No heating +gravity_power = 1.0d0 ! g ~ radius +/ +&Transport_Namelist +/ diff --git a/tests/coupled_bcs/run_test.sh b/tests/coupled_bcs/run_test.sh new file mode 100755 index 00000000..40e09c34 --- /dev/null +++ b/tests/coupled_bcs/run_test.sh @@ -0,0 +1,13 @@ +#!/usr/bin/env bash + +cd tests/coupled_bcs + +# we run a const case, where the bcs on the field that everything is coupled to are just constant +cd const +PYTHONPATH=../../../pre_processing:$PYTHONPATH python3 generate_input.py +mpirun -np 4 $RAYLEIGH_TEST_MPI_PARAMS ../../../bin/rayleigh.dbg +cd .. + +# after running, we test the output for errors +PYTHONPATH=../../post_processing:../../pre_processing:$PYTHONPATH python3 test_output.py + diff --git a/tests/coupled_bcs/test_output.py b/tests/coupled_bcs/test_output.py new file mode 100644 index 00000000..719d489a --- /dev/null +++ b/tests/coupled_bcs/test_output.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 + +from rayleigh_diagnostics import Shell_Slices, main_input +import numpy as np +import sys, os + +error = False + +def test_bcs(case): + """ + For a given case check if the values contained in the shell slices match the expected values based on the coefficients in + main_input + """ + error = False + + mi = main_input(os.path.join(case, 'main_input')) + bcs = mi.vals['boundary_conditions'] + slices = Shell_Slices('00000002', path=os.path.join(case, 'Shell_Slices/')) + + # values produced by the code + chi1_rmin = slices.vals[:,:,1,slices.lut[10001],0] + chi1_rmax = slices.vals[:,:,0,slices.lut[10001],0] + T_rmin = slices.vals[:,:,1,slices.lut[501],0] + T_rmax = slices.vals[:,:,0,slices.lut[501],0] + chi2_rmin = slices.vals[:,:,1,slices.lut[10201],0] + chi2_rmax = slices.vals[:,:,0,slices.lut[10201],0] + + # derivatives produced by the code + dchi1dr_rmin = slices.vals[:,:,1,slices.lut[10004],0] + dchi1dr_rmax = slices.vals[:,:,0,slices.lut[10004],0] + dTdr_rmin = slices.vals[:,:,1,slices.lut[507],0] + dTdr_rmax = slices.vals[:,:,0,slices.lut[507],0] + dchi2dr_rmin = slices.vals[:,:,1,slices.lut[10204],0] + dchi2dr_rmax = slices.vals[:,:,0,slices.lut[10204],0] + + # coefficients input + chi1bc_rmax = float(bcs['chi_a_top(1)'].replace('d', 'e')) + chi1bc_rmin = float(bcs['chi_a_bottom(1)'].replace('d', 'e')) + + dTdrbc_rmax = float(bcs['dtdr_top'].replace('d', 'e')) + dtdr_chi1_coeff_rmax = float(bcs['dtdr_chi_coeff_top(1)'].replace('d', 'e')) + + Tbc_rmin = float(bcs['t_bottom'].replace('d', 'e')) + t_chi1_coeff_rmin = float(bcs['t_chi_coeff_bottom(1)'].replace('d', 'e')) + t_dchi1dr_coeff_rmin = float(bcs['t_dchidr_coeff_bottom(1)'].replace('d', 'e')) + + chi2bc_rmax = float(bcs['chi_a_top(2)'].replace('d', 'e')) + chi2_chi1_coeff_rmax = float(bcs['chi_chi_coeff_top(2,1)'].replace('d', 'e')) + chi2_dchi1dr_coeff_rmax = float(bcs['chi_dchidr_coeff_top(2,1)'].replace('d', 'e')) + + dchi2drbc_rmin = float(bcs['dchidr_a_bottom(2)'].replace('d', 'e')) + dchi2dr_chi1_coeff_rmin = float(bcs['dchidr_chi_coeff_bottom(2,1)'].replace('d', 'e')) + + # check if values produced match expected values based on input coefficients + def failed(values, expected, tol=1.e-10): + return np.any(np.abs(values - expected) > tol) + + error = error or failed(chi1_rmax, chi1bc_rmax) + error = error or failed(chi1_rmin, chi1bc_rmin) + + error = error or failed(dTdr_rmax, dtdr_chi1_coeff_rmax*chi1_rmax + dTdrbc_rmax) + error = error or failed(T_rmin, t_chi1_coeff_rmin*chi1_rmin + t_dchi1dr_coeff_rmin*dchi1dr_rmin + Tbc_rmin) + + error = error or failed(chi2_rmax, chi2_chi1_coeff_rmax*chi1_rmax + chi2_dchi1dr_coeff_rmax*dchi1dr_rmax + chi2bc_rmax) + error = error or failed(dchi2dr_rmin, dchi2dr_chi1_coeff_rmin*chi1_rmin + dchi2drbc_rmin) + + return error + +lerror = test_bcs('const') +if lerror: + print("ERROR: const case produced unexpected output relative to input coefficients!") +error = error or lerror + +if error: + sys.exit(1) +sys.exit(0) + From e712b00ecdd2718e7a73238d4d10b4a0cdbfb2ae Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 16 Jun 2023 11:31:48 -0600 Subject: [PATCH 06/10] Randomize the numbers in main_input a bit to make the test more robust. --- tests/coupled_bcs/const/main_input | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/coupled_bcs/const/main_input b/tests/coupled_bcs/const/main_input index 380f9fc7..d6605d55 100644 --- a/tests/coupled_bcs/const/main_input +++ b/tests/coupled_bcs/const/main_input @@ -48,24 +48,24 @@ shellslice_nrec = 1 &Boundary_Conditions_Namelist no_slip_boundaries = .true. strict_L_Conservation = .false. -dTdr_Top = 10.0d0 -T_Bottom = 10.0d0 +dTdr_Top = 32.31 +T_Bottom = 17.865 couple_dtdr_top = .true. -dTdr_chi_coeff_top(1) = 1.0 +dTdr_chi_coeff_top(1) = 70.2 couple_tvar_bottom = .true. -T_chi_coeff_bottom(1) = 10.0 -T_dchidr_coeff_bottom(1) = -5.0 +T_chi_coeff_bottom(1) = 34.1 +T_dchidr_coeff_bottom(1) = -51.87 chi_a_Top(1) = 1.0d0 chi_a_Bottom(1) = 2.0d0 fix_chivar_a_top(1) = .true. fix_chivar_a_bottom(1) = .true. -chi_a_Top(2) = 1.0d0 -dchidr_a_Bottom(2) = 20.0d0 +chi_a_Top(2) = 1.0342d0 +dchidr_a_Bottom(2) = -27.6d0 couple_chivar_a_top(2) = .true. -chi_chi_coeff_top(2,1) = 10.0 -chi_dchidr_coeff_top(2,1) = 5.0 +chi_chi_coeff_top(2,1) = 107.01 +chi_dchidr_coeff_top(2,1) = 7.064 couple_dchidr_a_bottom(2) = .true. -dchidr_chi_coeff_bottom(2,1) = 1.0 +dchidr_chi_coeff_bottom(2,1) = 9.872 / &Initial_Conditions_Namelist init_type = 8 ! File init From 338b4e647c60a3cc1dd2e859f306dca59c8a26d0 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 16 Jun 2023 11:32:06 -0600 Subject: [PATCH 07/10] If no bcs are set default to Dirichlet. But unclear what to do about the default values. --- src/Physics/BoundaryConditions.F90 | 56 +++++++++++++++++++++++++----- 1 file changed, 48 insertions(+), 8 deletions(-) diff --git a/src/Physics/BoundaryConditions.F90 b/src/Physics/BoundaryConditions.F90 index d2f7df2a..d33a6db7 100755 --- a/src/Physics/BoundaryConditions.F90 +++ b/src/Physics/BoundaryConditions.F90 @@ -68,7 +68,7 @@ Module BoundaryConditions Logical :: Dipole_Field_Bottom = .False. Logical :: adjust_dTdr_Top = .True. - Real*8 :: T_Bottom = 1.0d0 + Real*8 :: T_Bottom = 0.0d0 Real*8 :: T_Top = 0.0d0 Real*8 :: dTdr_Top = 0.0d0 Real*8 :: dTdr_Bottom = 0.0d0 @@ -86,7 +86,7 @@ Module BoundaryConditions Real*8 :: dTdr_dchidr_coeff_bottom(1:n_scalar_max) = 0.0d0 Real*8 :: dTdr_T_coeff_bottom = 0.0d0 - Real*8 :: chi_a_Bottom(1:n_scalar_max) = 1.0d0 + Real*8 :: chi_a_Bottom(1:n_scalar_max) = 0.0d0 Real*8 :: chi_a_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_a_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_a_Bottom(1:n_scalar_max) = 0.0d0 @@ -182,39 +182,79 @@ Subroutine Initialize_Boundary_Conditions() Integer :: n_active_bcs n_active_bcs = count( (/ fix_tvar_top, fix_dtdr_top, couple_tvar_top, couple_dtdr_top /) ) - if (n_active_bcs /= 1) then + if (n_active_bcs > 1) then call stdout%print(" -- Error: Incompatible boundary conditions for tvar on top boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for tvar on top boundary.") + call stdout%print(" Defaulting to fix_tvar_top.") + call stdout%print(" ") + fix_tvar_top = .true. endif n_active_bcs = count( (/ fix_tvar_bottom, fix_dtdr_bottom, couple_tvar_bottom, couple_dtdr_bottom /) ) - if (n_active_bcs /= 1) then + if (n_active_bcs > 1) then call stdout%print(" -- Error: Incompatible boundary conditions for tvar on bottom boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for tvar on bottom boundary.") + call stdout%print(" Defaulting to fix_tvar_bottom.") + call stdout%print(" ") + fix_tvar_bottom = .true. endif do i = 1, n_active_scalars n_active_bcs = count( (/ fix_chivar_a_top(i), fix_dchidr_a_top(i), couple_chivar_a_top(i), couple_dchidr_a_top(i) /) ) - if (n_active_bcs /= 1) then + if (n_active_bcs > 1) then call stdout%print(" -- Error: Incompatible boundary conditions for chivar_a on top boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on top boundary.") + call stdout%print(" Defaulting to fix_chivar_a_top.") + call stdout%print(" ") + fix_chivar_a_top(i) = .true. endif n_active_bcs = count( (/ fix_chivar_a_bottom(i), fix_dchidr_a_bottom(i), & couple_chivar_a_bottom(i), couple_dchidr_a_bottom(i) /) ) - if (n_active_bcs /= 1) then + if (n_active_bcs > 1) then call stdout%print(" -- Error: Incompatible boundary conditions for chivar_a on bottom boundary.") + call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on bottom boundary.") + call stdout%print(" Defaulting to fix_chivar_a_bottom.") + call stdout%print(" ") + fix_chivar_a_bottom(i) = .true. endif end do do i = 1, n_passive_scalars n_active_bcs = count( (/ fix_chivar_p_top(i), fix_dchidr_p_top(i) /) ) - if (n_active_bcs /= 1) then + if (n_active_bcs > 1) then call stdout%print(" -- Error: Incompatible boundary conditions for chivar_p on top boundary.") call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on top boundary.") + call stdout%print(" Defaulting to fix_chivar_p_top.") + call stdout%print(" ") + fix_chivar_p_top(i) = .true. endif n_active_bcs = count( (/ fix_chivar_p_bottom(i), fix_dchidr_p_bottom(i) /) ) - if (n_active_bcs /= 1) then + if (n_active_bcs > 1) then call stdout%print(" -- Error: Incompatible boundary conditions for chivar_p on bottom boundary.") call stdout%print(" ") + call pfi%exit(1) + else if (n_active_bcs == 0) then + call stdout%print(" -- Warning: no boundary conditions set for chivar_a on bottom boundary.") + call stdout%print(" Defaulting to fix_chivar_p_bottom.") + call stdout%print(" ") + fix_chivar_p_bottom(i) = .true. endif end do From f246561ffdf7c19cd1872a035f186aed191bfb4d Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 21 Jun 2024 12:16:04 -0600 Subject: [PATCH 08/10] Adding documentation for coupled boundary conditions. Renaming the chi variables to explicitly state that we are only coupling chi_a. --- .../boundary_conditions_namelist.txt | 6 +- doc/source/User_Guide/coupled_bcs.txt | 181 ++++++++++++++++++ doc/source/User_Guide/under_development.rst | 8 + src/Physics/BoundaryConditions.F90 | 70 +++---- src/Physics/Sphere_Linear_Terms.F90 | 96 +++++----- 5 files changed, 276 insertions(+), 85 deletions(-) create mode 100644 doc/source/User_Guide/coupled_bcs.txt diff --git a/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt b/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt index 5dd45a39..cee8c4e7 100644 --- a/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt +++ b/doc/source/Namelist_Definitions/boundary_conditions_namelist.txt @@ -1,7 +1,7 @@ **fix_tvar_top** - Logical flag indicating whether thermal variable (T,S) should be fixed on the upper boundary. Default = .true. + Logical flag indicating whether thermal variable (T,S) should be fixed on the upper boundary. Default = .false. **fix_tvar_bottom** - Logical flag indicating whether thermal variable (T,S) should be fixed on the lower boundary. Default = .true. + Logical flag indicating whether thermal variable (T,S) should be fixed on the lower boundary. Default = .false. **fix_dtdr_top** Logical flag indicating whether the radial derivative of thermal variable (T,S) should be fixed on the upper boundary. Default = .false. **fix_dtdr_bottom** @@ -9,7 +9,7 @@ **T_top** Value of thermal variable (T,S) at the upper boundary. Default = 0. **T_bottom** - Value of thermal variable (T,S) at the lower boundary. Default = 1. + Value of thermal variable (T,S) at the lower boundary. Default = 0. **dTdr_top** Value of radial derivative of thermal variable (T,S) at the upper boundary. Default = 0. **dTdr_bottom** diff --git a/doc/source/User_Guide/coupled_bcs.txt b/doc/source/User_Guide/coupled_bcs.txt new file mode 100644 index 00000000..ef9cb6ac --- /dev/null +++ b/doc/source/User_Guide/coupled_bcs.txt @@ -0,0 +1,181 @@ +Rayleigh can couple the values and gradients of scalar fields (including the temperature/entropy, :math:`\Theta`, and active, +:math:`\chi_{a_i}`, scalar fields) at the +boundaries through this initial implementation of coupled boundary conditions. + +These boundary conditions take the form: + +.. math:: + :label: coupled_bcs_T + + \Theta = b_{\Theta,\Theta} + b_{\Theta,\frac{\partial \Theta}{\partial r}}\frac{\partial \Theta}{\partial r} \ + + \sum_{j} b_{\Theta,\chi_{a_j}} \chi_{a_j} \ + + \sum_{j} b_{\Theta,\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling temperature/entropy to its gradient and/or to other scalar fields, + +.. math:: + :label: coupled_bcs_dTdr + + \frac{\partial \Theta}{\partial r} = b_{\frac{\partial \Theta}{\partial r},\frac{\partial \Theta}{\partial r}} + b_{\frac{\partial \Theta}{\partial r},\Theta}\Theta \ + + \sum_{j} b_{\frac{\partial \Theta}{\partial r},\chi_{a_j}} \chi_{a_j} \ + + \sum_{j} b_{\frac{\partial \Theta}{\partial r},\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling the derivative of temperature/entropy to its value and/or to other scalar fields, + +.. math:: + :label: coupled_bcs_chia + + \chi_{a_i} = b_{\chi_{a_i},\chi_{a_i}} \ + + b_{\chi_{a_i},\Theta}\Theta \ + + b_{\chi_{a_i},\frac{\partial \Theta}{\partial r}}\frac{\partial \Theta}{\partial r} \ + + \sum_{j (j \ne i)} b_{\chi_{a_i},\chi_{a_j}} \chi_{a_j} \ + + \sum_{j} b_{\chi_{a_i},\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling active scalar i to its gradient and/or to other scalar fields including temperature/entropy, and finally + +.. math:: + :label: coupled_bcs_dchiadr + + \frac{\partial \chi_{a_i}}{\partial r} = b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}} \ + + b_{\frac{\partial\chi_{a_i}}{\partial r},\Theta}\Theta \ + + b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial \Theta}{\partial r}}\frac{\partial \Theta}{\partial r} \ + + \sum_{j} b_{\frac{\partial\chi_{a_i}}{\partial r},\chi_{a_j}} \chi_{a_j} \ + + \sum_{j (j \ne i)} b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_j}}{\partial r}} \frac{\partial\chi_{a_j}}{\partial r} + +for coupling the gradient of active scalar i to its value and/or to other scalar fields including temperature/entropy. + +The values of the :math:`b_{i,i}` terms can be constant or spatially varying. The values of the :math:`b_{i,j}` coefficients are constant. All values can be set using options in the boundary conditions namespace. + +Setup +^^^^^ + +Boundary conditions +******************* + +Model parameters for the scalar fields follow the same convention as temperature but using the prefix `chi_a` or `chi_p` for active and passive +scalars respectively. + +**couple_tvar_top** + Logical flag indicating whether thermal variable (T,S) should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_tvar_bottom** + Logical flag indicating whether thermal variable (T,S) should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**couple_dtdr_top** + Logical flag indicating whether radial gradient of thermal variable (T,S) should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_dtdr_bottom** + Logical flag indicating whether radial gradient of thermal variable (T,S) should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**couple_chivar_a_top** + Logical flag indicating whether active scalar i should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_chivar_a_bottom** + Logical flag indicating whether active scalar i should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**couple_dchidr_a_top** + Logical flag indicating whether radial gradient of active scalar i should be coupled to other scalar fields or their gradients on the upper boundary. Default = .false. +**couple_dchidr_a_bottom** + Logical flag indicating whether radial gradient of active scalar i should be coupled to other scalar fields or their gradients on the lower boundary. Default = .false. + +**T_top** + Set the thermal variable, :math:`b_{\Theta,\Theta}`, at the top of the domain (overloaded) +**T_top_file** + Set a spatially varying thermal variable, :math:`b_{\Theta,\Theta}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**T_dTdr_coeff_top** + Set the coupling coefficient between the thermal variable and the radial derivative of the thermal variable, :math:`b_{\Theta,\frac{\partial \Theta}{\partial r}}`, at the top of the domain +**T_chi_a_coeff_top(i)** + Set the coupling coefficient between the thermal variable and active scalar field i, :math:`b_{\Theta,\chi_{a_i}}`, at the top of the domain +**T_dchidr_a_coeff_top(i)** + Set the coupling coefficient between the thermal variable and the radial derivative of active scalar field i, :math:`b_{\Theta,\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain + +**T_bottom** + Set the thermal variable, :math:`b_{\Theta,\Theta}`, at the base of the domain (overloaded) +**T_bottom_file** + Set a spatially varying thermal variable, :math:`b_{\Theta,\Theta}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**T_dTdr_coeff_bottom** + Set the coupling coefficient between the thermal variable and the radial derivative of the thermal variable, :math:`b_{\Theta,\frac{\partial \Theta}{\partial r}}`, at the base of the domain +**T_chi_a_coeff_bottom(i)** + Set the coupling coefficient between the thermal variable and active scalar field i, :math:`b_{\Theta,\chi_{a_i}}`, at the base of the domain +**T_dchidr_a_coeff_bottom(i)** + Set the coupling coefficient between the thermal variable and the radial derivative of active scalar field i, :math:`b_{\Theta,\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain + +**dTdr_top** + Set the radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the top of the domain (overloaded) +**dTdr_top_file** + Set a spatially varying radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**dTdr_T_coeff_top** + Set the coupling coefficient between the radial derivative of the thermal variable and the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\Theta}`, at the top of the domain +**dTdr_chi_a_coeff_top(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain +**dTdr_dchidr_a_coeff_top(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain + +**dTdr_bottom** + Set the radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the base of the domain (overloaded) +**dTdr_bottom_file** + Set a spatially varying radial derivative of the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**dTdr_T_coeff_bottom** + Set the coupling coefficient between the radial derivative of the thermal variable and the thermal variable, :math:`b_{\frac{\partial\Theta}{\partial r},\Theta}`, at the base of the domain +**dTdr_chi_a_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain +**dTdr_dchidr_a_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of the thermal variable and the radial derivative of active scalar field i, :math:`b_{\frac{\partial\Theta}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain + + +**chi_a_top(i)** + Set active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the top of the domain (overloaded) +**chi_a_top_file(i)** + Set a spatially varying active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**chi_a_T_coeff_top(i)** + Set the coupling coefficient between active scalar field i and the thermal variable, :math:`b_{\chi_{a_i},\Theta}`, at the top of the domain +**chi_a_dTdr_coeff_top(i)** + Set the coupling coefficient between active scalar field i and the radial derivative of the thermal variable, :math:`b_{\chi_{a_i},\frac{\partial\Theta}{\partial r}}`, at the top of the domain +**chi_a_chi_a_coeff_top(i,j)** + Set the coupling coefficient between active scalar field i and active scalar field j, :math:`b_{\chi_{a_i},\chi_{a_j}}`, at the top of the domain. :math:`b_{\chi_{a_i},\chi_{a_i}}` is ignored +**chi_a_dchidr_a_coeff_top(i,j)** + Set the coupling coefficient between active scalar field i and the radial derivative of active scalar field j, :math:`b_{\chi_{a_i},\frac{\partial \chi_{a_j}}{\partial r}}`, at the top of the domain + +**chi_a_bottom(i)** + Set active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the base of the domain (overloaded) +**chi_a_bottom_file(i)** + Set a spatially varying active scalar field i, :math:`b_{\chi_{a_i},\chi_{a_i}}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**chi_a_T_coeff_bottom(i)** + Set the coupling coefficient between active scalar field i and the thermal variable, :math:`b_{\chi_{a_i},\Theta}`, at the base of the domain +**chi_a_dTdr_coeff_bottom(i)** + Set the coupling coefficient between active scalar field i and the radial derivative of the thermal variable, :math:`b_{\chi_{a_i},\frac{\partial\Theta}{\partial r}}`, at the base of the domain +**chi_a_chi_a_coeff_bottom(i,j)** + Set the coupling coefficient between active scalar field i and active scalar field j, :math:`b_{\chi_{a_i},\chi_{a_j}}`, at the base of the domain. :math:`b_{\chi_{a_i},\chi_{a_i}}` is ignored +**chi_a_dchidr_a_coeff_bottom(i,j)** + Set the coupling coefficient between active scalar field i and the radial derivative of active scalar field j, :math:`b_{\chi_{a_i},\frac{\partial \chi_{a_j}}{\partial r}}`, at the base of the domain + +**dchidr_a_top(i)** + Set the radial derivative of active scalar field i, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the top of the domain (overloaded) +**dchidr_a_top_file(i)** + Set a spatially varying radial derivative of active scalar field i, :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}`, at the top of the domain by specifying a generic input filename (overloaded, untested) +**dchidr_a_T_coeff_top(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\Theta}`, at the top of the domain +**dchidr_a_dTdr_coeff_top(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the top of the domain +**dchidr_a_chi_a_coeff_top(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\chi_{a_j}}`, at the top of the domain +**dchidr_a_dchidr_a_coeff_top(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_j}}{\partial r}}`, at the top of the domain. :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}` is ignored + +**dchidr_a_bottom(i)** + Set the radial derivative of active scalar field i, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_i}}{\partial r}}`, at the base of the domain (overloaded) +**dchidr_a_bottom_file(i)** + Set a spatially varying radial derivative of active scalar field i, :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}`, at the base of the domain by specifying a generic input filename (overloaded, untested) +**dchidr_a_T_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\Theta}`, at the base of the domain +**dchidr_a_dTdr_coeff_bottom(i)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of the thermal variable, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial\Theta}{\partial r}}`, at the base of the domain +**dchidr_a_chi_a_coeff_bottom(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\chi_{a_j}}`, at the base of the domain +**dchidr_a_dchidr_a_coeff_bottom(i,j)** + Set the coupling coefficient between the radial derivative of active scalar field i and the radial derivative of active scalar field j, :math:`b_{\frac{\partial \chi_{a_i}}{\partial r},\frac{\partial \chi_{a_j}}{\partial r}}`, at the base of the domain. :math:`b_{\frac{\partial\chi_{a_i}}{\partial r},\frac{\partial\chi_{a_i}}{\partial r}}` is ignored + +Further Information +^^^^^^^^^^^^^^^^^^^ + +See `tests/coupled_bcs` for example input files. + + + diff --git a/doc/source/User_Guide/under_development.rst b/doc/source/User_Guide/under_development.rst index 8c4e31d4..2ae76b40 100644 --- a/doc/source/User_Guide/under_development.rst +++ b/doc/source/User_Guide/under_development.rst @@ -17,6 +17,14 @@ Arbitrary Scalar Fields .. include:: arbitrary_scalar_fields.txt +.. _coupled_bcs: + +Coupled Boundary Conditions +--------------------------- + +.. include:: coupled_bcs.txt + + diff --git a/src/Physics/BoundaryConditions.F90 b/src/Physics/BoundaryConditions.F90 index d33a6db7..7a19ed83 100755 --- a/src/Physics/BoundaryConditions.F90 +++ b/src/Physics/BoundaryConditions.F90 @@ -73,17 +73,17 @@ Module BoundaryConditions Real*8 :: dTdr_Top = 0.0d0 Real*8 :: dTdr_Bottom = 0.0d0 - Real*8 :: T_chi_coeff_top(1:n_scalar_max) = 0.0d0 - Real*8 :: T_dchidr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: T_chi_a_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dchidr_a_coeff_top(1:n_scalar_max) = 0.0d0 Real*8 :: T_dTdr_coeff_top = 0.0d0 - Real*8 :: T_chi_coeff_bottom(1:n_scalar_max) = 0.0d0 - Real*8 :: T_dchidr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: T_chi_a_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: T_dchidr_a_coeff_bottom(1:n_scalar_max) = 0.0d0 Real*8 :: T_dTdr_coeff_bottom = 0.0d0 - Real*8 :: dTdr_chi_coeff_top(1:n_scalar_max) = 0.0d0 - Real*8 :: dTdr_dchidr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_chi_a_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_dchidr_a_coeff_top(1:n_scalar_max) = 0.0d0 Real*8 :: dTdr_T_coeff_top = 0.0d0 - Real*8 :: dTdr_chi_coeff_bottom(1:n_scalar_max) = 0.0d0 - Real*8 :: dTdr_dchidr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_chi_a_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dTdr_dchidr_a_coeff_bottom(1:n_scalar_max) = 0.0d0 Real*8 :: dTdr_T_coeff_bottom = 0.0d0 Real*8 :: chi_a_Bottom(1:n_scalar_max) = 0.0d0 @@ -91,22 +91,22 @@ Module BoundaryConditions Real*8 :: dchidr_a_Top(1:n_scalar_max) = 0.0d0 Real*8 :: dchidr_a_Bottom(1:n_scalar_max) = 0.0d0 - Real*8 :: chi_chi_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: chi_dchidr_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: chi_T_coeff_top(1:n_scalar_max) = 0.0d0 - Real*8 :: chi_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 - Real*8 :: chi_chi_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: chi_dchidr_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: chi_T_coeff_bottom(1:n_scalar_max) = 0.0d0 - Real*8 :: chi_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_chi_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_dchidr_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_T_coeff_top(1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_chi_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_dchidr_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_T_coeff_bottom(1:n_scalar_max) = 0.0d0 - Real*8 :: dchidr_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_chi_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dchidr_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_T_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_chi_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dchidr_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_T_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: chi_a_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_chi_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dchidr_a_coeff_top(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_T_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dTdr_coeff_top(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_chi_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dchidr_a_coeff_bottom(1:n_scalar_max,1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_T_coeff_bottom(1:n_scalar_max) = 0.0d0 + Real*8 :: dchidr_a_dTdr_coeff_bottom(1:n_scalar_max) = 0.0d0 Real*8 :: chi_p_Bottom(1:n_scalar_max) = 1.0d0 Real*8 :: chi_p_Top(1:n_scalar_max) = 0.0d0 @@ -162,15 +162,15 @@ Module BoundaryConditions fix_chivar_p_top, fix_chivar_p_bottom, chi_p_bottom, chi_p_top, & fix_dchidr_p_bottom, fix_dchidr_p_top, dchidr_p_top, dchidr_p_bottom, & couple_Tvar_top, couple_Tvar_bottom, couple_dtdr_top, couple_dtdr_bottom, & - T_chi_coeff_top, T_dchidr_coeff_top, T_dTdr_coeff_top, & - T_chi_coeff_bottom, T_dchidr_coeff_bottom, T_dTdr_coeff_bottom, & - dTdr_chi_coeff_top, dTdr_dchidr_coeff_top, dTdr_T_coeff_top, & - dTdr_chi_coeff_bottom, dTdr_dchidr_coeff_bottom, dTdr_T_coeff_bottom, & + T_chi_a_coeff_top, T_dchidr_a_coeff_top, T_dTdr_coeff_top, & + T_chi_a_coeff_bottom, T_dchidr_a_coeff_bottom, T_dTdr_coeff_bottom, & + dTdr_chi_a_coeff_top, dTdr_dchidr_a_coeff_top, dTdr_T_coeff_top, & + dTdr_chi_a_coeff_bottom, dTdr_dchidr_a_coeff_bottom, dTdr_T_coeff_bottom, & couple_chivar_a_top, couple_chivar_a_bottom, couple_dchidr_a_top, couple_dchidr_a_bottom, & - chi_chi_coeff_top, chi_dchidr_coeff_top, chi_T_coeff_top, chi_dTdr_coeff_top, & - chi_chi_coeff_bottom, chi_dchidr_coeff_bottom, chi_T_coeff_bottom, chi_dTdr_coeff_bottom, & - dchidr_chi_coeff_top, dchidr_dchidr_coeff_top, dchidr_T_coeff_top, dchidr_dTdr_coeff_top, & - dchidr_chi_coeff_bottom, dchidr_dchidr_coeff_bottom, dchidr_T_coeff_bottom, dchidr_dTdr_coeff_bottom + chi_a_chi_a_coeff_top, chi_a_dchidr_a_coeff_top, chi_a_T_coeff_top, chi_a_dTdr_coeff_top, & + chi_a_chi_a_coeff_bottom, chi_a_dchidr_a_coeff_bottom, chi_a_T_coeff_bottom, chi_a_dTdr_coeff_bottom, & + dchidr_a_chi_a_coeff_top, dchidr_a_dchidr_a_coeff_top, dchidr_a_T_coeff_top, dchidr_a_dTdr_coeff_top, & + dchidr_a_chi_a_coeff_bottom, dchidr_a_dchidr_a_coeff_bottom, dchidr_a_T_coeff_bottom, dchidr_a_dTdr_coeff_bottom Contains @@ -188,9 +188,10 @@ Subroutine Initialize_Boundary_Conditions() call pfi%exit(1) else if (n_active_bcs == 0) then call stdout%print(" -- Warning: no boundary conditions set for tvar on top boundary.") - call stdout%print(" Defaulting to fix_tvar_top.") + call stdout%print(" Defaulting to fix_tvar_top with T_top = 0.0.") call stdout%print(" ") fix_tvar_top = .true. + T_top = 0.0 endif n_active_bcs = count( (/ fix_tvar_bottom, fix_dtdr_bottom, couple_tvar_bottom, couple_dtdr_bottom /) ) @@ -200,9 +201,10 @@ Subroutine Initialize_Boundary_Conditions() call pfi%exit(1) else if (n_active_bcs == 0) then call stdout%print(" -- Warning: no boundary conditions set for tvar on bottom boundary.") - call stdout%print(" Defaulting to fix_tvar_bottom.") + call stdout%print(" Defaulting to fix_tvar_bottom with T_bottom = 1.0.") call stdout%print(" ") fix_tvar_bottom = .true. + T_bottom = 1.0 endif do i = 1, n_active_scalars diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90 index 00f12d85..fb35ad02 100755 --- a/src/Physics/Sphere_Linear_Terms.F90 +++ b/src/Physics/Sphere_Linear_Terms.F90 @@ -560,9 +560,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_tvar_top) Then do i = 1, n_active_scalars - samp = -T_chi_coeff_top(i) + samp = -T_chi_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -T_dchidr_coeff_top(i) + samp = -T_dchidr_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = 1.0 @@ -572,9 +572,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_dtdr_top) Then do i = 1, n_active_scalars - samp = -dTdr_chi_coeff_top(i) + samp = -dTdr_chi_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -dTdr_dchidr_coeff_top(i) + samp = -dTdr_dchidr_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = -dTdr_T_coeff_top @@ -592,9 +592,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_tvar_bottom) Then do i = 1, n_active_scalars - samp = -T_chi_coeff_bottom(i) + samp = -T_chi_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -T_dchidr_coeff_bottom(i) + samp = -T_dchidr_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = 1.0 @@ -604,9 +604,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_dtdr_bottom) Then do i = 1, n_active_scalars - samp = -dTdr_chi_coeff_bottom(i) + samp = -dTdr_chi_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -dTdr_dchidr_coeff_bottom(i) + samp = -dTdr_dchidr_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = -dTdr_T_coeff_bottom @@ -626,28 +626,28 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_chivar_a_top(i)) Then do j = 1, n_active_scalars - samp = -chi_chi_coeff_top(i,j) + samp = -chi_a_chi_a_coeff_top(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -chi_dchidr_coeff_top(i,j) + samp = -chi_a_dchidr_a_coeff_top(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -chi_T_coeff_top(i) + samp = -chi_a_T_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -chi_dTdr_coeff_top(i) + samp = -chi_a_dTdr_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif If (couple_dchidr_a_top(i)) Then do j = 1, n_active_scalars - samp = -dchidr_chi_coeff_top(i,j) + samp = -dchidr_a_chi_a_coeff_top(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -dchidr_dchidr_coeff_top(i,j) + samp = -dchidr_a_dchidr_a_coeff_top(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -dchidr_T_coeff_top(i) + samp = -dchidr_a_T_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -dchidr_dTdr_coeff_top(i) + samp = -dchidr_a_dTdr_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif @@ -660,28 +660,28 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_chivar_a_bottom(i)) Then do j = 1, n_active_scalars - samp = -chi_chi_coeff_bottom(i,j) + samp = -chi_a_chi_a_coeff_bottom(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -chi_dchidr_coeff_bottom(i,j) + samp = -chi_a_dchidr_a_coeff_bottom(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -chi_T_coeff_bottom(i) + samp = -chi_a_T_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -chi_dTdr_coeff_bottom(i) + samp = -chi_a_dTdr_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif If (couple_dchidr_a_bottom(i)) Then do j = 1, n_active_scalars - samp = -dchidr_chi_coeff_bottom(i,j) + samp = -dchidr_a_chi_a_coeff_bottom(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -dchidr_dchidr_coeff_bottom(i,j) + samp = -dchidr_a_dchidr_a_coeff_bottom(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -dchidr_T_coeff_bottom(i) + samp = -dchidr_a_T_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -dchidr_dTdr_coeff_bottom(i) + samp = -dchidr_a_dTdr_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif end do @@ -766,9 +766,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_tvar_top) Then do i = 1, n_active_scalars - samp = -T_chi_coeff_top(i) + samp = -T_chi_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -T_dchidr_coeff_top(i) + samp = -T_dchidr_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = 1.0 @@ -778,9 +778,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_dtdr_top) Then do i = 1, n_active_scalars - samp = -dTdr_chi_coeff_top(i) + samp = -dTdr_chi_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -dTdr_dchidr_coeff_top(i) + samp = -dTdr_dchidr_a_coeff_top(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = -dTdr_T_coeff_top @@ -798,9 +798,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_tvar_bottom) Then do i = 1, n_active_scalars - samp = -T_chi_coeff_bottom(i) + samp = -T_chi_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -T_dchidr_coeff_bottom(i) + samp = -T_dchidr_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = 1.0 @@ -810,9 +810,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_dtdr_bottom) Then do i = 1, n_active_scalars - samp = -dTdr_chi_coeff_bottom(i) + samp = -dTdr_chi_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,0) - samp = -dTdr_dchidr_coeff_bottom(i) + samp = -dTdr_dchidr_a_coeff_bottom(i) Call Load_BC(lp,r,teq,chiavar(i),samp,1) end do samp = -dTdr_T_coeff_bottom @@ -832,28 +832,28 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_chivar_a_top(i)) Then do j = 1, n_active_scalars - samp = -chi_chi_coeff_top(i,j) + samp = -chi_a_chi_a_coeff_top(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -chi_dchidr_coeff_top(i,j) + samp = -chi_a_dchidr_a_coeff_top(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -chi_T_coeff_top(i) + samp = -chi_a_T_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -chi_dTdr_coeff_top(i) + samp = -chi_a_dTdr_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif If (couple_dchidr_a_top(i)) Then do j = 1, n_active_scalars - samp = -dchidr_chi_coeff_top(i,j) + samp = -dchidr_a_chi_a_coeff_top(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -dchidr_dchidr_coeff_top(i,j) + samp = -dchidr_a_dchidr_a_coeff_top(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -dchidr_T_coeff_top(i) + samp = -dchidr_a_T_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -dchidr_dTdr_coeff_top(i) + samp = -dchidr_a_dTdr_coeff_top(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif @@ -866,28 +866,28 @@ Subroutine Set_Boundary_Conditions(mode_ind) Endif If (couple_chivar_a_bottom(i)) Then do j = 1, n_active_scalars - samp = -chi_chi_coeff_bottom(i,j) + samp = -chi_a_chi_a_coeff_bottom(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -chi_dchidr_coeff_bottom(i,j) + samp = -chi_a_dchidr_a_coeff_bottom(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -chi_T_coeff_bottom(i) + samp = -chi_a_T_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -chi_dTdr_coeff_bottom(i) + samp = -chi_a_dTdr_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif If (couple_dchidr_a_bottom(i)) Then do j = 1, n_active_scalars - samp = -dchidr_chi_coeff_bottom(i,j) + samp = -dchidr_a_chi_a_coeff_bottom(i,j) Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,0) - samp = -dchidr_dchidr_coeff_bottom(i,j) + samp = -dchidr_a_dchidr_a_coeff_bottom(i,j) if (i==j) samp = 1.0 Call Load_BC(lp,r,chiaeq(i),chiavar(j),samp,1) end do - samp = -dchidr_T_coeff_bottom(i) + samp = -dchidr_a_T_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,0) - samp = -dchidr_dTdr_coeff_bottom(i) + samp = -dchidr_a_dTdr_coeff_bottom(i) Call Load_BC(lp,r,chiaeq(i),tvar,samp,1) Endif end do From 8ecfae563c9cf3d8c85e7f03444416ff898a174f Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 21 Jun 2024 12:21:46 -0600 Subject: [PATCH 09/10] Update the changelog. --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 700b079a..b61ae319 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm ## [Unreleased] ### Added +- Allow the boundary conditions of scalar fields (both the thermal temperature/entropy field and active scalar fields) to be linearly coupled together. This is documented in the "Under Development" section of the User Guide and has an example test. \[Cian Wilson; 6-21-2024; [#460](https://github.com/geodynamics/Rayleigh/pull/460)\] + - The documentation for solving for active and passive scalar fields has been expanded in the "Under Development" section of the User Guide. \[Cian Wilson; 6-20-2024; [#541](https://github.com/geodynamics/Rayleigh/pull/541)\] - A docker image for Stampede3 and Frontera (TACC) has been added to the repository. \[Rene Gassmoeller; 6-20-2024; [#402](https://github.com/geodynamics/Rayleigh/pull/402)\] From 4d52a8dee7ca991a1066462d63aa4ed1149bf113 Mon Sep 17 00:00:00 2001 From: Cian Wilson Date: Fri, 21 Jun 2024 14:43:17 -0600 Subject: [PATCH 10/10] Update the test to the new options names. --- tests/coupled_bcs/const/main_input | 12 ++++++------ tests/coupled_bcs/test_output.py | 24 ++++++++++++------------ 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/tests/coupled_bcs/const/main_input b/tests/coupled_bcs/const/main_input index d6605d55..988215e0 100644 --- a/tests/coupled_bcs/const/main_input +++ b/tests/coupled_bcs/const/main_input @@ -51,10 +51,10 @@ strict_L_Conservation = .false. dTdr_Top = 32.31 T_Bottom = 17.865 couple_dtdr_top = .true. -dTdr_chi_coeff_top(1) = 70.2 +dTdr_chi_a_coeff_top(1) = 70.2 couple_tvar_bottom = .true. -T_chi_coeff_bottom(1) = 34.1 -T_dchidr_coeff_bottom(1) = -51.87 +T_chi_a_coeff_bottom(1) = 34.1 +T_dchidr_a_coeff_bottom(1) = -51.87 chi_a_Top(1) = 1.0d0 chi_a_Bottom(1) = 2.0d0 fix_chivar_a_top(1) = .true. @@ -62,10 +62,10 @@ fix_chivar_a_bottom(1) = .true. chi_a_Top(2) = 1.0342d0 dchidr_a_Bottom(2) = -27.6d0 couple_chivar_a_top(2) = .true. -chi_chi_coeff_top(2,1) = 107.01 -chi_dchidr_coeff_top(2,1) = 7.064 +chi_a_chi_a_coeff_top(2,1) = 107.01 +chi_a_dchidr_a_coeff_top(2,1) = 7.064 couple_dchidr_a_bottom(2) = .true. -dchidr_chi_coeff_bottom(2,1) = 9.872 +dchidr_a_chi_a_coeff_bottom(2,1) = 9.872 / &Initial_Conditions_Namelist init_type = 8 ! File init diff --git a/tests/coupled_bcs/test_output.py b/tests/coupled_bcs/test_output.py index 719d489a..ea995193 100644 --- a/tests/coupled_bcs/test_output.py +++ b/tests/coupled_bcs/test_output.py @@ -34,22 +34,22 @@ def test_bcs(case): dchi2dr_rmax = slices.vals[:,:,0,slices.lut[10204],0] # coefficients input - chi1bc_rmax = float(bcs['chi_a_top(1)'].replace('d', 'e')) - chi1bc_rmin = float(bcs['chi_a_bottom(1)'].replace('d', 'e')) + chi1bc_rmax = float(bcs['chi_a_top(1)'.lower()].replace('d', 'e')) + chi1bc_rmin = float(bcs['chi_a_bottom(1)'.lower()].replace('d', 'e')) - dTdrbc_rmax = float(bcs['dtdr_top'].replace('d', 'e')) - dtdr_chi1_coeff_rmax = float(bcs['dtdr_chi_coeff_top(1)'].replace('d', 'e')) + dTdrbc_rmax = float(bcs['dTdr_top'.lower()].replace('d', 'e')) + dtdr_chi1_coeff_rmax = float(bcs['dTdr_chi_a_coeff_top(1)'.lower()].replace('d', 'e')) - Tbc_rmin = float(bcs['t_bottom'].replace('d', 'e')) - t_chi1_coeff_rmin = float(bcs['t_chi_coeff_bottom(1)'].replace('d', 'e')) - t_dchi1dr_coeff_rmin = float(bcs['t_dchidr_coeff_bottom(1)'].replace('d', 'e')) + Tbc_rmin = float(bcs['T_bottom'.lower()].replace('d', 'e')) + t_chi1_coeff_rmin = float(bcs['T_chi_a_coeff_bottom(1)'.lower()].replace('d', 'e')) + t_dchi1dr_coeff_rmin = float(bcs['T_dchidr_a_coeff_bottom(1)'.lower()].replace('d', 'e')) - chi2bc_rmax = float(bcs['chi_a_top(2)'].replace('d', 'e')) - chi2_chi1_coeff_rmax = float(bcs['chi_chi_coeff_top(2,1)'].replace('d', 'e')) - chi2_dchi1dr_coeff_rmax = float(bcs['chi_dchidr_coeff_top(2,1)'].replace('d', 'e')) + chi2bc_rmax = float(bcs['chi_a_top(2)'.lower()].replace('d', 'e')) + chi2_chi1_coeff_rmax = float(bcs['chi_a_chi_a_coeff_top(2,1)'.lower()].replace('d', 'e')) + chi2_dchi1dr_coeff_rmax = float(bcs['chi_a_dchidr_a_coeff_top(2,1)'.lower()].replace('d', 'e')) - dchi2drbc_rmin = float(bcs['dchidr_a_bottom(2)'].replace('d', 'e')) - dchi2dr_chi1_coeff_rmin = float(bcs['dchidr_chi_coeff_bottom(2,1)'].replace('d', 'e')) + dchi2drbc_rmin = float(bcs['dchidr_a_bottom(2)'.lower()].replace('d', 'e')) + dchi2dr_chi1_coeff_rmin = float(bcs['dchidr_a_chi_a_coeff_bottom(2,1)'.lower()].replace('d', 'e')) # check if values produced match expected values based on input coefficients def failed(values, expected, tol=1.e-10):