From c5568c8fcf56ea96f585e835a86f8e5042f08784 Mon Sep 17 00:00:00 2001 From: BWHindman Date: Fri, 21 Jun 2024 19:02:32 -0600 Subject: [PATCH 1/3] Controls and part of PDE_Coefficients modified. --- src/Physics/Controls.F90 | 5 ++- src/Physics/PDE_Coefficients.F90 | 52 +++++++++++++++++++++++++++++++- 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 1b314aae..9df1d77e 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -63,6 +63,7 @@ Module Controls Logical :: lorentz_forces = .true. ! Turn Lorentz forces on or off (default is on - as long as magnetism is on) Logical :: viscous_heating = .true. ! Turns viscous heating on/off Logical :: ohmic_heating = .true. + Logical :: pseudo_incompressible = .false. ! Switch from anelastic to pseudo-incompressible approximation Logical :: advect_reference_state = .true. ! Set to true to advect the reference state temperature or entropy ! This has no effect for adiabatic reference states. ! Generally only do this if reference state is nonadiabatic @@ -78,7 +79,8 @@ Module Controls Namelist /Physical_Controls_Namelist/ magnetism, nonlinear, rotation, lorentz_forces, & & viscous_heating, ohmic_heating, advect_reference_state, benchmark_mode, & & benchmark_integration_interval, benchmark_report_interval, & - & momentum_advection, inertia, n_active_scalars, n_passive_scalars + & momentum_advection, inertia, n_active_scalars, n_passive_scalars, + & pseudo_incompressible !/////////////////////////////////////////////////////////////////////////// ! Temporal Controls @@ -216,6 +218,7 @@ Subroutine Restore_Physics_Defaults() lorentz_forces = .true. viscous_heating = .true. ohmic_heating = .true. + pseudo_incompressible = .false. advect_reference_state = .true. benchmark_mode = 0 benchmark_integration_interval = -1 diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 97ca8108..8d9c5262 100644 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -52,7 +52,11 @@ Module PDE_Coefficients Real*8, Allocatable :: Temperature(:) Real*8, Allocatable :: dlnT(:) + Real*8, Allocatable :: entropy(:) + Real*8, Allocatable :: exp_entropy(:) Real*8, Allocatable :: dsdr(:) + Real*8, Allocatable :: dsdr_over_cp(:) + Real*8, Allocatable :: d2s(:) Real*8, Allocatable :: heating(:) @@ -262,7 +266,11 @@ Subroutine Allocate_Reference_State Allocate(ref%dlnrho(1:N_R)) Allocate(ref%d2lnrho(1:N_R)) Allocate(ref%dlnt(1:N_R)) + Allocate(ref%entropy(1:N_R)) + Allocate(ref%exp_entropy(1:N_R)) Allocate(ref%dsdr(1:N_R)) + Allocate(ref%dsdr_over_cp(1:N_R)) + Allocate(ref%d2s(1:N_R)) Allocate(ref%Buoyancy_Coeff(1:N_R)) Allocate(ref%dpdr_w_term(1:N_R)) Allocate(ref%pressure_dwdr_term(1:N_R)) @@ -287,7 +295,11 @@ Subroutine Allocate_Reference_State ref%dlnrho(:) = Zero ref%d2lnrho(:) = Zero ref%dlnt(:) = Zero + ref%entropy(:) = Zero + ref%exp_entropy(:) = Zero ref%dsdr(:) = Zero + ref%dsdr_over_cp(:) = Zero + ref%d2s(:) = Zero ref%buoyancy_coeff(:) = Zero ref%dpdr_w_term(:) = Zero ref%pressure_dwdr_term(:) = Zero @@ -348,7 +360,11 @@ Subroutine Constant_Reference() ref%d2lnrho = 0.0d0 ref%temperature = 1.0d0 ref%dlnT = 0.0d0 + ref%entropy = 0.0d0 + ref%exp_entropy = 1.0d0 ref%dsdr = 0.0d0 + ref%dsdr_over_cp = 0.0d0 + ref%d2s = 0.0d0 amp = Rayleigh_Number/Prandtl_Number @@ -473,7 +489,12 @@ Subroutine Polytropic_ReferenceND() dtmparr = (poly_n/ref%temperature)*(2.0d0*Dissipation_Number*gravity/radius) ! (n/T)*d2Tdr2 ref%d2lnrho = ref%d2lnrho+dtmparr + ref%entropy(:) = 0.0d0 + ref%exp_entropy(:) = 1.0d0 ref%dsdr(:) = 0.0d0 + ref%dsdr_over_cp(:) = 0.0d0 + ref%d2s(:) = 0.0d0 + Call Initialize_Reference_Heating() ref%Coriolis_Coeff = 2.0d0 @@ -738,6 +759,8 @@ Subroutine Polytropic_ReferenceND_General() ref%density = zeta**poly_n ref%temperature = zeta + ref%entropy = (1.0d0/specific_heat_ratio)*(log(ref%Temperature) - (specific_heat_ratio - 1.0d0)*log(ref%density)) + ref%exp_entropy = exp(ref%entropy) ref%dlnrho = poly_n*dlnzeta ref%dlnT = dlnzeta @@ -748,6 +771,8 @@ Subroutine Polytropic_ReferenceND_General() ref%dsdr = (1.0d0/Specific_Heat_Ratio)*(ref%dlnT - (Specific_Heat_Ratio - 1.0d0) * ref%dlnrho) ! This is (1/c_p) dS/dr (where "S" is dimensional background S) ! That's fine up to a multiplicative constant which we determine below + ref%dsdr_over_cp = ref%dsdr + ref%d2s = (1.0d0/specific_heat_ratio)*(d2lnzeta - (specific_heat_ratio - 1.0d0) * ref%d2lnrho) nsquared = gravity*ref%dsdr ! N^2 (non-dimensional) up to multiplicative constant If (.not. adiabatic_polytrope) Then ! don't divide by zero! @@ -807,6 +832,9 @@ Subroutine Polytropic_ReferenceND_General() ! These are the same no matter how we non-dimensionalize time ref%dsdr = (Prandtl_Number*Buoyancy_Number_Visc/Rayleigh_Number) * nsquared/gravity + ref%dsdr_over_cp = ref%dsdr + !ref%d2s = ? BWH: It is unclear to me whether the second derivative needs to be renormalized. + ref%dpdr_w_term(:) = ref%density(:) ref%pressure_dwdr_term(:) = -ref%density(:) @@ -901,6 +929,7 @@ Subroutine Polytropic_Reference() Real*8 :: beta Real*8 :: Gravitational_Constant = 6.67d-8 ! cgs units Real*8, Allocatable :: zeta(:), gravity(:) + Real*8, Allocatable :: dlnzeta(:), d2lnzeta(:) Real*8 :: One Real*8 :: InnerRadius, OuterRadius Character*12 :: dstring @@ -949,10 +978,13 @@ Subroutine Polytropic_Reference() ! also rho_c, T_c, P_c Allocate(zeta(N_R), gravity(1:N_R)) + Allocate(dlnzeta(1:N_R), d2lnzeta(1:N_R)) d = OuterRadius - InnerRadius zeta = c0 + c1 * d / Radius + dlnzeta = -c1*d/Radius**2/zeta + d2lnzeta = 2.0d0*c1*d/Radius**3/zeta - (c1*d/Radius**2/zeta)**2 rho_c = poly_rho_i / zeta(N_R)**poly_n @@ -974,9 +1006,13 @@ Subroutine Polytropic_Reference() Ref%d2lnrho = - Ref%dlnrho*(2.0d0/Radius-c1*d/zeta/Radius**2) Ref%Temperature = T_c * zeta - Ref%dlnT = -(c1*d/Radius**2)/zeta + Ref%dlnT = dlnzeta + ref%entropy = volume_specific_heat * (log(ref%Temperature) - (specific_heat_ratio - 1.0d0) * log(ref%Temperature) - (specific_heat_ratio - 1.0d0) * log(ref%density)) + ref%exp_entropy = exp(ref%entropy/pressure_specific_heat) Ref%dsdr = volume_specific_heat * (Ref%dlnT - (Specific_Heat_Ratio - 1.0d0) * Ref%dlnrho) + ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat + ref%d2s = volume_specific_heat * (d2lnzeta - (specific_heat_ratio - 1.0d0) * ref%d2lnrho) Ref%Buoyancy_Coeff = gravity/Pressure_Specific_Heat*ref%density @@ -985,6 +1021,7 @@ Subroutine Polytropic_Reference() end do Deallocate(zeta, gravity) + Deallocate(dlnzeta, d2lnzeta) Call Initialize_Reference_Heating() @@ -1203,6 +1240,7 @@ Subroutine Get_Custom_Reference() Integer :: i, fi Character(len=2) :: ind Integer :: fi_to_check(4) = (/1, 2, 4, 6/) + Real*8 :: geofac If (my_rank .eq. 0) Then Write(6,*)'Custom reference state specified.' @@ -1254,6 +1292,14 @@ Subroutine Get_Custom_Reference() ref%ohmic_amp(:) = ra_constants(9)/(ref%density(:)*ref%temperature(:)) ref%dsdr(:) = ra_constants(11)*ra_functions(:,14) + ! Integrate dsdr to obtain a self-consistent entropy profile. Set s=0 at the upper surface. + ref%entropy(1) = 0.0 + Do i = 2, N_R + ref%entropy(i) = ref%entropy(i-1) - geofac*ref%dsdr(i)*radial_integral_weights(i)/Radius(i)**2 + Enddo + ref%exp_entropy = exp(ref%entropy/pressure_specific_heat) + ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat + Call log_deriv(ref%dsdr(:), ref%d2s(i), no_log=.true.) End Subroutine Get_Custom_Reference @@ -1652,7 +1698,11 @@ Subroutine Restore_Reference_Defaults If (allocated(ref%Temperature)) DeAllocate(ref%Temperature) If (allocated(ref%dlnT)) DeAllocate(ref%dlnT) + If (allocated(ref%entropy)) DeAllocate(ref%entropy) + If (allocated(ref%exp_entropy)) DeAllocate(ref%exp_entropy) If (allocated(ref%dsdr)) DeAllocate(ref%dsdr) + If (allocated(ref%dsdr_over_cp)) DeAllocate(ref%dsdr_over_cp) + If (allocated(ref%d2s)) DeAllocate(ref%d2s) If (allocated(ref%heating)) DeAllocate(ref%heating) From 59514cb39d67883c71eaf330c508d92e543b7461 Mon Sep 17 00:00:00 2001 From: BWHindman Date: Mon, 24 Jun 2024 14:03:40 -0600 Subject: [PATCH 2/3] Revision of the pseudo-incompressible equations. --- src/Physics/PDE_Coefficients.F90 | 61 ++++++++++++++++++--------- src/Physics/Sphere_Hybrid_Space.F90 | 16 ++++++- src/Physics/Sphere_Linear_Terms.F90 | 33 ++++++++++++++- src/Physics/Sphere_Physical_Space.F90 | 34 ++++++++++++++- 4 files changed, 118 insertions(+), 26 deletions(-) diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 8d9c5262..35b8886c 100644 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -52,11 +52,11 @@ Module PDE_Coefficients Real*8, Allocatable :: Temperature(:) Real*8, Allocatable :: dlnT(:) - Real*8, Allocatable :: entropy(:) - Real*8, Allocatable :: exp_entropy(:) + Real*8, Allocatable :: entropy(:) ! Entropy s, with s=0 on the outer boundary + Real*8, Allocatable :: exp_entropy(:) ! exp(s/c_P) Real*8, Allocatable :: dsdr(:) - Real*8, Allocatable :: dsdr_over_cp(:) - Real*8, Allocatable :: d2s(:) + Real*8, Allocatable :: dsdr_over_cp(:) ! (1/c_P) ds/dr + Real*8, Allocatable :: d2s_over_cp(:) ! (1/c_P) d2s/dr2 Real*8, Allocatable :: heating(:) @@ -270,7 +270,7 @@ Subroutine Allocate_Reference_State Allocate(ref%exp_entropy(1:N_R)) Allocate(ref%dsdr(1:N_R)) Allocate(ref%dsdr_over_cp(1:N_R)) - Allocate(ref%d2s(1:N_R)) + Allocate(ref%d2s_over_cp(1:N_R)) Allocate(ref%Buoyancy_Coeff(1:N_R)) Allocate(ref%dpdr_w_term(1:N_R)) Allocate(ref%pressure_dwdr_term(1:N_R)) @@ -299,7 +299,7 @@ Subroutine Allocate_Reference_State ref%exp_entropy(:) = Zero ref%dsdr(:) = Zero ref%dsdr_over_cp(:) = Zero - ref%d2s(:) = Zero + ref%d2s_over_cp(:) = Zero ref%buoyancy_coeff(:) = Zero ref%dpdr_w_term(:) = Zero ref%pressure_dwdr_term(:) = Zero @@ -364,7 +364,7 @@ Subroutine Constant_Reference() ref%exp_entropy = 1.0d0 ref%dsdr = 0.0d0 ref%dsdr_over_cp = 0.0d0 - ref%d2s = 0.0d0 + ref%d2s_over_cp = 0.0d0 amp = Rayleigh_Number/Prandtl_Number @@ -489,11 +489,11 @@ Subroutine Polytropic_ReferenceND() dtmparr = (poly_n/ref%temperature)*(2.0d0*Dissipation_Number*gravity/radius) ! (n/T)*d2Tdr2 ref%d2lnrho = ref%d2lnrho+dtmparr - ref%entropy(:) = 0.0d0 - ref%exp_entropy(:) = 1.0d0 - ref%dsdr(:) = 0.0d0 + ref%entropy(:) = 0.0d0 + ref%exp_entropy(:) = 1.0d0 + ref%dsdr(:) = 0.0d0 ref%dsdr_over_cp(:) = 0.0d0 - ref%d2s(:) = 0.0d0 + ref%d2s_over_cp(:) = 0.0d0 Call Initialize_Reference_Heating() @@ -760,7 +760,8 @@ Subroutine Polytropic_ReferenceND_General() ref%density = zeta**poly_n ref%temperature = zeta ref%entropy = (1.0d0/specific_heat_ratio)*(log(ref%Temperature) - (specific_heat_ratio - 1.0d0)*log(ref%density)) - ref%exp_entropy = exp(ref%entropy) + ref%entropy = ref%entropy - ref%entropy(1) ! Set the entropy to zero at the upper surface + ref%exp_entropy = exp(ref%entropy) ! Note the entropy is dimensionless: entropy = S/c_P ref%dlnrho = poly_n*dlnzeta ref%dlnT = dlnzeta @@ -772,7 +773,7 @@ Subroutine Polytropic_ReferenceND_General() ! This is (1/c_p) dS/dr (where "S" is dimensional background S) ! That's fine up to a multiplicative constant which we determine below ref%dsdr_over_cp = ref%dsdr - ref%d2s = (1.0d0/specific_heat_ratio)*(d2lnzeta - (specific_heat_ratio - 1.0d0) * ref%d2lnrho) + ref%d2s_over_cp = (1.0d0/specific_heat_ratio)*(d2lnzeta - (specific_heat_ratio - 1.0d0) * ref%d2lnrho) nsquared = gravity*ref%dsdr ! N^2 (non-dimensional) up to multiplicative constant If (.not. adiabatic_polytrope) Then ! don't divide by zero! @@ -832,8 +833,8 @@ Subroutine Polytropic_ReferenceND_General() ! These are the same no matter how we non-dimensionalize time ref%dsdr = (Prandtl_Number*Buoyancy_Number_Visc/Rayleigh_Number) * nsquared/gravity - ref%dsdr_over_cp = ref%dsdr - !ref%d2s = ? BWH: It is unclear to me whether the second derivative needs to be renormalized. + !ref%dsdr_over_cp = ref%dsdr + !ref%d2s_over_cp = ? BWH: It is unclear to me whether derivatives need to be renormalized. ref%dpdr_w_term(:) = ref%density(:) ref%pressure_dwdr_term(:) = -ref%density(:) @@ -1008,11 +1009,13 @@ Subroutine Polytropic_Reference() Ref%Temperature = T_c * zeta Ref%dlnT = dlnzeta - ref%entropy = volume_specific_heat * (log(ref%Temperature) - (specific_heat_ratio - 1.0d0) * log(ref%Temperature) - (specific_heat_ratio - 1.0d0) * log(ref%density)) - ref%exp_entropy = exp(ref%entropy/pressure_specific_heat) + ! Set the entropy to zero at the upper surface + ref%entropy = volume_specific_heat * (log(ref%Temperature) - (specific_heat_ratio - 1.0d0) * log(ref%density)) + ref%entropy = ref%entropy - ref%entropy(1) ! zeroing the entropy at the upper boundary + ref%exp_entropy = exp(ref%entropy/pressure_specific_heat) ! Used extensively in the pseudo-incompressible approximation Ref%dsdr = volume_specific_heat * (Ref%dlnT - (Specific_Heat_Ratio - 1.0d0) * Ref%dlnrho) ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat - ref%d2s = volume_specific_heat * (d2lnzeta - (specific_heat_ratio - 1.0d0) * ref%d2lnrho) + ref%d2s_over_cp = (1.0d0/specific_heat_ratio)*(d2lnzeta - (specific_heat_ratio - 1.0d0) * ref%d2lnrho) Ref%Buoyancy_Coeff = gravity/Pressure_Specific_Heat*ref%density @@ -1293,13 +1296,14 @@ Subroutine Get_Custom_Reference() ref%dsdr(:) = ra_constants(11)*ra_functions(:,14) ! Integrate dsdr to obtain a self-consistent entropy profile. Set s=0 at the upper surface. - ref%entropy(1) = 0.0 + ref%entropy(1) = 0.0d0 + geofac = (radius(1)**3 - radius(N_R)**3)/3.0d0 Do i = 2, N_R ref%entropy(i) = ref%entropy(i-1) - geofac*ref%dsdr(i)*radial_integral_weights(i)/Radius(i)**2 Enddo ref%exp_entropy = exp(ref%entropy/pressure_specific_heat) ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat - Call log_deriv(ref%dsdr(:), ref%d2s(i), no_log=.true.) + Call log_deriv(ref%dsdr_over_cp(:), ref%d2s_over_cp(i), no_log=.true.) End Subroutine Get_Custom_Reference @@ -1702,7 +1706,7 @@ Subroutine Restore_Reference_Defaults If (allocated(ref%exp_entropy)) DeAllocate(ref%exp_entropy) If (allocated(ref%dsdr)) DeAllocate(ref%dsdr) If (allocated(ref%dsdr_over_cp)) DeAllocate(ref%dsdr_over_cp) - If (allocated(ref%d2s)) DeAllocate(ref%d2s) + If (allocated(ref%d2s_over_cp)) DeAllocate(ref%d2s_over_cp) If (allocated(ref%heating)) DeAllocate(ref%heating) @@ -2029,6 +2033,7 @@ Subroutine Compute_Diffusion_Coefs() W_Diffusion_Coefs_0 = -nu*(4.0d0/3.0d0)*( dlnu*ref%dlnrho + ref%d2lnrho + ref%dlnrho/radius + & & 3.0d0*dlnu/radius ) W_Diffusion_Coefs_1 = nu*(2.0d0*dlnu-ref%dlnrho/3.0d0) + !///////////////////////////////////// ! W Coefficients for dWdr equation @@ -2059,6 +2064,20 @@ Subroutine Compute_Diffusion_Coefs() Z_Diffusion_Coefs_0 = -nu*( 2.0d0*dlnu/radius + ref%dlnrho*dlnu + & & ref%d2lnrho+2.0d0*ref%dlnrho/radius) Z_Diffusion_Coefs_1 = nu*(dlnu-ref%dlnrho) + + + If (pseudo_incompressible) Then + W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 - nu*(4.0d0/3.0d0)*ref%dsdr_over_cp*(dnu - ref%dlnrho -ref%dsdr_over_cp - 2.0d0/radius) + W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 - nu*(4.0d0/3.0d0)*ref%d2s_over_cp + W_Diffusion_Coefs_1 = W_Diffusion_Coefs_1 - nu*(7.0d0/3.0d0)*ref%dsdr_over_cp + DW_Diffusion_Coefs_0 = DW_Diffusion_Coefs_0 + nu*ref%dsdr_over_cp/3.0d0 + DW_Diffusion_Coefs_1 = DW_Diffusion_Coefs_1 - nu*ref%dsdr_over_cp*(dlnu - ref%dlnrho - ref%dsdr_over_cp) + DW_Diffusion_Coefs_1 = DW_Diffusion_Coefs_1 - nu*ref%d2s_over_cp + DW_Diffusion_Coefs_2 = DW_Diffusion_Coefs_2 - 2.0d0*nu*ref%dsdr_over_cp + Z_Diffusion_Coefs_0 = Z_Diffusion_Coefs_0 - nu*ref%dsdr_over_cp*(dlnu - ref%dlnrho - ref%dsdr_over_cp) + Z_Diffusion_Coefs_0 = Z_Diffusion_Coefs_0 - nu*ref%d2s_over_cp + Z_Diffusion_Coefs_1 = Z_Diffusion_Coefs_1 - 2.0d0*nu*ref%dsdr_over_cp + Endif !//////////////////////////////////////// ! A (vector potential) Coefficients diff --git a/src/Physics/Sphere_Hybrid_Space.F90 b/src/Physics/Sphere_Hybrid_Space.F90 index 1e84829f..0178a1df 100755 --- a/src/Physics/Sphere_Hybrid_Space.F90 +++ b/src/Physics/Sphere_Hybrid_Space.F90 @@ -52,6 +52,11 @@ Subroutine Hybrid_Init() over_rhor(r1:r2) = one_over_r(r1:r2)/ref%density(r1:r2) over_rhorsq(r1:r2) = OneOverRSquared(r1:r2)/ref%density(r1:r2) drho_term(r1:r2) = ref%dlnrho(r1:r2)+one_over_r(r1:r2) + If (pseudo_incompressible) Then + over_rhosq(r1:r2) = over_rhosq(r1:r2)/ ref%exp_entropy(r1:r2) + over_rhor(r1:r2) = over_rhor(r1:r2) / ref%exp_entropy(r1:r2) + drho_term(r1:r2) = drho_term(r1:r2) + ref%dsdr_over_cp(r1:r2) + Endif End Subroutine Hybrid_Init @@ -283,7 +288,8 @@ Subroutine Velocity_Derivatives() ! .... Small correction for density variation : - u_theta*dlnrhodr (added -u_theta/r as well here) ! Notice that there is a -u_theta/r term above. These should be combined - ! for efficiency later + ! for efficiency later. + ! The pseudo-incompressible correction has been implemented through drho_term. DO_IDX2 SBUFFA(IDX2,dvtdr) = SBUFFA(IDX2,dvtdr)- & @@ -305,6 +311,7 @@ Subroutine Velocity_Derivatives() ! .... Small correction for density variation : - u_phi*dlnrhodr ! .... moved -u_phi/r here as well + ! ... pseudo-incompressible correction has been implemented through drho_term. DO_IDX2 SBUFFA(IDX2,dvpdr) = SBUFFA(IDX2,dvpdr)- & & SBUFFA(IDX2,vphi)*drho_term(r) @@ -328,12 +335,19 @@ Subroutine Velocity_Derivatives() SBUFFA(IDX2,dvrdr) = SBUFFA(IDX2,dvrdr)- & & SBUFFA(IDX2,vr)*ref%dlnrho(r) END_DO + If (pseudo_incompressible) Then + DO_IDX2 + SBUFFA(IDX2,dvrdr) = SBUFFA(IDX2,dvrdr)- & + & SBUFFA(IDX2,vr)*ref%dsdr_over_cp(r) + END_DO + Endif Call d_by_dtheta(wsp%s2a,vr,dvrdt) ! Convert Z to ell(ell+1) Z/r^2 (i.e. omega_r) + ! Computing the radial component of the curl of the velocity DO_IDX2 SBUFFA(IDX2,zvar) = l_l_plus1(m:l_max)*SBUFFA(IDX2,zvar)*Over_RhoRSQ(r) END_DO diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90 index fb35ad02..c58ee026 100755 --- a/src/Physics/Sphere_Linear_Terms.F90 +++ b/src/Physics/Sphere_Linear_Terms.F90 @@ -282,7 +282,7 @@ Subroutine Load_Linear_Coefficients() Call add_implicit_term(chipeq(i),chipvar(i), 1, amp,lp) end do - Else + Else ! l > 0 !================================================== ! Radial Momentum Equation @@ -291,20 +291,40 @@ Subroutine Load_Linear_Coefficients() ! Temperature amp = -ref%Buoyancy_Coeff/H_Laplacian + If (pseudo_incompressible) Then + amp = amp*ref%exp_entropy + Endif Call add_implicit_term(weq, tvar, 0, amp,lp) ! Gravity ! Chi do i = 1, n_active_scalars amp = -ref%chi_buoyancy_coeff(i,:)/H_Laplacian + If (pseudo_incompressible) Then + amp = amp*ref%exp_entropy + Endif Call add_implicit_term(weq, chiavar(i), 0, amp,lp) ! Gravity end do - ! Pressure + ! Pressure Force !amp = 1.0d0/(Ek*H_Laplacian)*ref%density ! dPdr amp = ref%dpdr_W_term/H_Laplacian + If (pseudo_incompressible) Then + amp = amp*ref%exp_entropy + Endif Call add_implicit_term(weq,pvar, 1, amp,lp) + + ! Add the buoyancy term ignored under the LBR approximation + ! -(ds/dr) rho/(c_P * H_Laplacian) (P/rho) + ! amp = -rho/(c_P * H_Laplacian) (ds/dr) Non-LBR Anelastic (not implemented) + ! amp = -exp(s/c_P) rho/(c_P * H_Laplacian) (ds/dr) Pseudo-incompressible + If (pseudo_incompressible) Then + amp = -ref%exp_entropy * ref%density * ref%dsdr_over_cp / H_Laplacian + add_implicit_term(weq,pvar, 0, amp, lp) + Endif + + ! W If (inertia) Then @@ -333,6 +353,9 @@ Subroutine Load_Linear_Coefficients() ! Pressure !amp = -(1.0d0)/Ek*ref%density amp = ref%pressure_dwdr_term + If (pseudo_incompressible) Then + amp = amp*ref%exp_entropy + Endif Call add_implicit_term(peq,pvar, 0, amp,lp) ! W @@ -980,6 +1003,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) ! Else stress-free r = 1 samp = -(2.0d0/radius(r)+ref%dlnrho(r)) + If (pseudo_incompressible) Then + samp = samp - ref%dsdr_over_cp(r) + Endif Call Load_BC(lp,r,peq,wvar,one,2) Call Load_BC(lp,r,peq,wvar,samp,1) @@ -997,6 +1023,9 @@ Subroutine Set_Boundary_Conditions(mode_ind) !stress_free_bottom r = N_R samp = -(2.0d0/radius(r)+ref%dlnrho(r)) + If (pseudo_incompressible) Then + samp = samp - ref%dsdr_over_cp(r) + Endif Call Load_BC(lp,r,peq,wvar,one,2) Call Load_BC(lp,r,peq,wvar,samp,1) Call Load_BC(lp,r,zeq,zvar,one,1) diff --git a/src/Physics/Sphere_Physical_Space.F90 b/src/Physics/Sphere_Physical_Space.F90 index 968bf644..992e8a7d 100755 --- a/src/Physics/Sphere_Physical_Space.F90 +++ b/src/Physics/Sphere_Physical_Space.F90 @@ -445,6 +445,14 @@ Subroutine Momentum_Advection_Radial() !$OMP END PARALLEL DO Endif + ! Multiply by rho_*/rho = exp(s/c_P) if pseudo-incompressible + If (pseudo_incompressible) Then + !$OMP PARALLEL DO PRIVATE(t,r,k) + DO_IDX + RHSP(IDX,wvar) = RHSP(IDX,wvar)*ref%exp_entropy(r) + END_DO + !$OMP END PARALLEL DO + Endif End Subroutine Momentum_Advection_Radial @@ -561,8 +569,6 @@ Subroutine Momentum_Advection_Theta() Endif - - ! At this point, we have [u dot grad u]_theta ! Multiply by radius/sintheta so that we have r[u dot grad u]_theta/sintheta (getting ready for Z and dWdr RHS building) !$OMP PARALLEL DO PRIVATE(t,r,k) @@ -570,10 +576,22 @@ Subroutine Momentum_Advection_Theta() RHSP(IDX,pvar) = RHSP(IDX,pvar)*radius(r)*csctheta(t) END_DO !$OMP END PARALLEL DO + + + ! Multiply by rho_*/rho = exp(s/c_P) if pseudo-incompressible + If (pseudo_incompressible) Then + !$OMP PARALLEL DO PRIVATE(t,r,k) + DO_IDX + RHSP(IDX,pvar) = RHSP(IDX,pvar)*ref%exp_entropy(r) + END_DO + !$OMP END PARALLEL DO + Endif End Subroutine Momentum_Advection_Theta + + Subroutine Momentum_Advection_Phi() Implicit None Integer :: t, r, k @@ -637,7 +655,19 @@ Subroutine Momentum_Advection_Phi() RHSP(IDX,zvar) = RHSP(IDX,zvar)*radius(r)*csctheta(t) END_DO !OMP END PARALLEL DO + + ! Multiply by rho_*/rho = exp(s/c_P) if pseudo-incompressible + If (pseudo_incompressible) Then + !$OMP PARALLEL DO PRIVATE(t,r,k) + DO_IDX + RHSP(IDX,pvar) = RHSP(IDX,pvar)*ref%exp_entropy(r) + END_DO + !$OMP END PARALLEL DO + Endif + End Subroutine Momentum_Advection_Phi + + Subroutine Phi_Derivatives() Implicit None Integer :: i From 09b7b6796eb7f2b96afc6635e43b2b7ce6d6b4c1 Mon Sep 17 00:00:00 2001 From: BWHindman Date: Mon, 24 Jun 2024 15:07:40 -0600 Subject: [PATCH 3/3] Fix until compiled. --- src/Physics/Controls.F90 | 2 +- src/Physics/PDE_Coefficients.F90 | 5 +++-- src/Physics/Sphere_Hybrid_Space.F90 | 13 +++++++------ src/Physics/Sphere_Linear_Terms.F90 | 2 +- 4 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 9df1d77e..d7f1d43c 100755 --- a/src/Physics/Controls.F90 +++ b/src/Physics/Controls.F90 @@ -79,7 +79,7 @@ Module Controls Namelist /Physical_Controls_Namelist/ magnetism, nonlinear, rotation, lorentz_forces, & & viscous_heating, ohmic_heating, advect_reference_state, benchmark_mode, & & benchmark_integration_interval, benchmark_report_interval, & - & momentum_advection, inertia, n_active_scalars, n_passive_scalars, + & momentum_advection, inertia, n_active_scalars, n_passive_scalars, & & pseudo_incompressible !/////////////////////////////////////////////////////////////////////////// diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90 index 35b8886c..b89a34b0 100644 --- a/src/Physics/PDE_Coefficients.F90 +++ b/src/Physics/PDE_Coefficients.F90 @@ -1303,7 +1303,7 @@ Subroutine Get_Custom_Reference() Enddo ref%exp_entropy = exp(ref%entropy/pressure_specific_heat) ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat - Call log_deriv(ref%dsdr_over_cp(:), ref%d2s_over_cp(i), no_log=.true.) + Call log_deriv(ref%dsdr_over_cp(:), ref%d2s_over_cp(:), no_log=.true.) End Subroutine Get_Custom_Reference @@ -2067,7 +2067,8 @@ Subroutine Compute_Diffusion_Coefs() If (pseudo_incompressible) Then - W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 - nu*(4.0d0/3.0d0)*ref%dsdr_over_cp*(dnu - ref%dlnrho -ref%dsdr_over_cp - 2.0d0/radius) + W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 - nu*(4.0d0/3.0d0)*ref%dsdr_over_cp*(dlnu - ref%dlnrho -ref%dsdr_over_cp) + W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 + nu*(8.0d0/3.0d0)*ref%dsdr_over_cp/radius W_Diffusion_Coefs_0 = W_Diffusion_Coefs_0 - nu*(4.0d0/3.0d0)*ref%d2s_over_cp W_Diffusion_Coefs_1 = W_Diffusion_Coefs_1 - nu*(7.0d0/3.0d0)*ref%dsdr_over_cp DW_Diffusion_Coefs_0 = DW_Diffusion_Coefs_0 + nu*ref%dsdr_over_cp/3.0d0 diff --git a/src/Physics/Sphere_Hybrid_Space.F90 b/src/Physics/Sphere_Hybrid_Space.F90 index 0178a1df..8d24dfdc 100755 --- a/src/Physics/Sphere_Hybrid_Space.F90 +++ b/src/Physics/Sphere_Hybrid_Space.F90 @@ -49,13 +49,14 @@ Subroutine Hybrid_Init() Allocate(drho_term(my_r%min:my_r%max)) r1 = my_r%min r2 = my_r%max - over_rhor(r1:r2) = one_over_r(r1:r2)/ref%density(r1:r2) - over_rhorsq(r1:r2) = OneOverRSquared(r1:r2)/ref%density(r1:r2) - drho_term(r1:r2) = ref%dlnrho(r1:r2)+one_over_r(r1:r2) If (pseudo_incompressible) Then - over_rhosq(r1:r2) = over_rhosq(r1:r2)/ ref%exp_entropy(r1:r2) - over_rhor(r1:r2) = over_rhor(r1:r2) / ref%exp_entropy(r1:r2) - drho_term(r1:r2) = drho_term(r1:r2) + ref%dsdr_over_cp(r1:r2) + over_rhor(r1:r2) = one_over_r(r1:r2)/(ref%density(r1:r2)*ref%exp_entropy(r1:r2)) + over_rhorsq(r1:r2)= OneOverRSquared(r1:r2)/(ref%density(r1:r2)*ref%exp_entropy(r1:r2)) + drho_term(r1:r2) = ref%dlnrho(r1:r2)+one_over_r(r1:r2) + ref%dsdr_over_cp(r1:r2) + Else + over_rhor(r1:r2) = one_over_r(r1:r2)/ref%density(r1:r2) + over_rhorsq(r1:r2)= OneOverRSquared(r1:r2)/ref%density(r1:r2) + drho_term(r1:r2) = ref%dlnrho(r1:r2)+one_over_r(r1:r2) Endif End Subroutine Hybrid_Init diff --git a/src/Physics/Sphere_Linear_Terms.F90 b/src/Physics/Sphere_Linear_Terms.F90 index c58ee026..a95a05e6 100755 --- a/src/Physics/Sphere_Linear_Terms.F90 +++ b/src/Physics/Sphere_Linear_Terms.F90 @@ -321,7 +321,7 @@ Subroutine Load_Linear_Coefficients() ! amp = -exp(s/c_P) rho/(c_P * H_Laplacian) (ds/dr) Pseudo-incompressible If (pseudo_incompressible) Then amp = -ref%exp_entropy * ref%density * ref%dsdr_over_cp / H_Laplacian - add_implicit_term(weq,pvar, 0, amp, lp) + Call add_implicit_term(weq,pvar, 0, amp, lp) Endif