diff --git a/src/Physics/Controls.F90 b/src/Physics/Controls.F90 index 1b314aae..d7f1d43c 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..b89a34b0 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(:) ! 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(:) ! (1/c_P) ds/dr + Real*8, Allocatable :: d2s_over_cp(:) ! (1/c_P) d2s/dr2 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_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)) @@ -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_over_cp(:) = 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_over_cp = 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%dsdr(:) = 0.0d0 + ref%entropy(:) = 0.0d0 + ref%exp_entropy(:) = 1.0d0 + ref%dsdr(:) = 0.0d0 + ref%dsdr_over_cp(:) = 0.0d0 + ref%d2s_over_cp(:) = 0.0d0 + Call Initialize_Reference_Heating() ref%Coriolis_Coeff = 2.0d0 @@ -738,6 +759,9 @@ 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%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 @@ -748,6 +772,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_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! @@ -807,6 +833,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_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(:) @@ -901,6 +930,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 +979,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 +1007,15 @@ 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 + ! 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_over_cp = (1.0d0/specific_heat_ratio)*(d2lnzeta - (specific_heat_ratio - 1.0d0) * ref%d2lnrho) Ref%Buoyancy_Coeff = gravity/Pressure_Specific_Heat*ref%density @@ -985,6 +1024,7 @@ Subroutine Polytropic_Reference() end do Deallocate(zeta, gravity) + Deallocate(dlnzeta, d2lnzeta) Call Initialize_Reference_Heating() @@ -1203,6 +1243,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 +1295,15 @@ 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.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_over_cp(:), ref%d2s_over_cp(:), no_log=.true.) End Subroutine Get_Custom_Reference @@ -1652,7 +1702,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_over_cp)) DeAllocate(ref%d2s_over_cp) If (allocated(ref%heating)) DeAllocate(ref%heating) @@ -1979,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 @@ -2009,6 +2064,21 @@ 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*(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 + 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..8d24dfdc 100755 --- a/src/Physics/Sphere_Hybrid_Space.F90 +++ b/src/Physics/Sphere_Hybrid_Space.F90 @@ -49,9 +49,15 @@ 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_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 @@ -283,7 +289,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 +312,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 +336,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..a95a05e6 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 + Call 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