[wip] Simplified and minimal update to include pseudo-incompressibility

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

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

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

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)