Skip to content

Commit

Permalink
Eccentricity tides + cleanups
Browse files Browse the repository at this point in the history
  • Loading branch information
@AnkitBarik
AnkitBarik committed Jun 15, 2024
1 parent 041c413 commit e9bb616
Show file tree
Hide file tree
Showing 11 changed files with 131 additions and 46 deletions.
11 changes: 10 additions & 1 deletion src/Namelists.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,8 @@ subroutine readNamelists(tscheme)
& omega_ma1,omegaOsz_ma1,tShift_ma1, &
& omega_ma2,omegaOsz_ma2,tShift_ma2, &
& amp_mode_ma,omega_mode_ma,m_mode_ma, &
& mode_symm_ma,ellipticity_cmb
& mode_symm_ma,ellipticity_cmb, amp_tide, &
& omega_tide

namelist/inner_core/sigma_ratio,nRotIc,rho_ratio_ic, &
& omega_ic1,omegaOsz_ic1,tShift_ic1, &
Expand Down Expand Up @@ -741,6 +742,10 @@ subroutine readNamelists(tscheme)
l_AM=l_AM .or. l_correct_AMe .or. l_correct_AMz
l_AM=l_AM .or. l_power

! Radial flow BC?
if (ellipticity_cmb /= 0.0_cp .or. ellipticity_icb /= 0.0_cp .or. &
& amp_tide /=0.0_cp ) l_radial_flow_bc=.true.

!-- Heat boundary condition
if ( impS /= 0 ) then
rad=pi/180
Expand Down Expand Up @@ -1223,6 +1228,8 @@ subroutine writeNamelists(n_out)
write(n_out,'('' m_mode_ma ='',i4,'','')') m_mode_ma
write(n_out,'('' mode_symm_ma ='',i4,'','')') mode_symm_ma
write(n_out,'('' ellipticity_cmb ='',ES14.6,'','')') ellipticity_cmb
write(n_out,'('' amp_tide ='',ES14.6,'','')') amp_tide
write(n_out,'('' omega_tide ='',ES14.6,'','')') omega_tide
write(n_out,*) "/"

write(n_out,*) "&inner_core"
Expand Down Expand Up @@ -1354,6 +1361,8 @@ subroutine defaultNamelists
radratio =0.35_cp
dilution_fac=0.0_cp ! centrifugal acceleration
ampForce =0.0_cp ! External body force amplitude
amp_tide =0.0_cp ! Amplitude of tide
omega_tide =0.0_cp ! Frequency of tide

!-- Phase field
tmelt =0.0_cp ! Melting temperature
Expand Down
2 changes: 1 addition & 1 deletion src/fields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ module fields
use precision_mod
use mem_alloc, only: bytes_allocated
use constants, only: zero
use physical_parameters, only: ampForce
use special, only: ampForce
use truncation, only: lm_max, n_r_max, lm_maxMag, n_r_maxMag, &
& n_r_ic_maxMag, fd_order, fd_order_bound
use logic, only: l_chemical_conv, l_finite_diff, l_mag, l_parallel_solve, &
Expand Down
3 changes: 2 additions & 1 deletion src/init_fields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ module init_fields
& impXi, n_impXi_max, n_impXi, phiXi, &
& thetaXi, peakXi, widthXi, osc, epscxi, &
& kbotxi, ktopxi, BuoFac, ktopp, oek, &
& epsPhase, ampForce
& epsPhase
use special, only: ampForce
use algebra, only: prepare_mat, solve_mat
use cosine_transform_odd
use dense_matrices
Expand Down
6 changes: 0 additions & 6 deletions src/phys_param.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,12 +108,6 @@ module physical_parameters
real(cp) :: penaltyFac ! Factor that enters the penalty method in the NS equations
real(cp) :: tmelt ! Melting temperature

real(cp) :: ellipticity_cmb ! Ellipticity of CMB, used for libration
real(cp) :: ellipticity_icb ! Ellipticity of ICB, used for libration
real(cp) :: ellip_fac_cmb ! d/d\phi (Y22) * d/dt exp(i\omega t) at CMB
real(cp) :: ellip_fac_icb ! d/d\phi (Y22) * d/dt exp(i\omega t) at ICB

real(cp) :: ampForce ! Amplitude of external body force
real(cp) :: gammatau_gravi ! Constant in front of the core/mantle gravitationnal torque (see Aubert et al. 2013)

end module physical_parameters
45 changes: 31 additions & 14 deletions src/preCalculations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,16 @@ module preCalculations
& ktops, kbots, interior_model, r_LCR, &
& n_r_LCR, mode, tmagcon, oek, Bn, imagcon,&
& ktopxi, kbotxi, epscxi, epscxi0, sc, osc,&
& ChemFac, raxi, Po, prec_angle, &
& ellipticity_cmb, ellip_fac_cmb, &
& ellipticity_icb, ellip_fac_icb
& ChemFac, raxi, Po, prec_angle
use horizontal_data, only: horizontal
use integration, only: rInt_R
use useful, only: logWrite, abortRun
use special, only: l_curr, fac_loop, loopRadRatio, amp_curr, Le, n_imp, l_imp
use special, only: l_curr, fac_loop, loopRadRatio, amp_curr, Le, n_imp, &
& l_imp, l_radial_flow_bc, &
& ellipticity_cmb, ellip_fac_cmb, &
& ellipticity_icb, ellip_fac_icb, &
& tide_fac20, tide_fac22p, tide_fac22n, &
& amp_tide
use time_schemes, only: type_tscheme

implicit none
Expand Down Expand Up @@ -75,6 +78,7 @@ subroutine preCalc(tscheme)
character(len=76) :: fileName
character(len=80) :: message
real(cp) :: mom(n_r_max)
real(cp) :: y20_norm, y22_norm

!-- Determine scales depending on n_tScale,n_lScale :
if ( n_tScale == 0 ) then
Expand Down Expand Up @@ -763,18 +767,31 @@ subroutine preCalc(tscheme)
call abortRun('LCR not compatible with imposed field!')
end if

if ( ellipticity_cmb /= 0.0_cp ) then
ellip_fac_cmb=-two*r_cmb**3*ellipticity_cmb*omega_ma1*omegaOsz_ma1
else
ellip_fac_cmb=0.0_cp
end if
if (l_radial_flow_bc) then
if ( ellipticity_cmb /= 0.0_cp ) then
ellip_fac_cmb=-two*r_cmb**3*ellipticity_cmb*omega_ma1*omegaOsz_ma1
else
ellip_fac_cmb=0.0_cp
end if

if ( ellipticity_icb /= 0.0_cp ) then
ellip_fac_icb=-two*r_icb**3*ellipticity_icb*omega_ic1*omegaOsz_ic1
else
ellip_fac_icb=0.0_cp
end if
if ( ellipticity_icb /= 0.0_cp ) then
ellip_fac_icb=-two*r_icb**3*ellipticity_icb*omega_ic1*omegaOsz_ic1
else
ellip_fac_icb=0.0_cp
end if

if ( amp_tide /= 0.0_cp ) then
y20_norm = 0.5_cp * sqrt(5.0_cp/pi)
y22_norm = 0.25_cp * sqrt(7.5_cp/pi)
tide_fac20 = amp_tide / y20_norm * r_cmb**2
tide_fac22p = half * amp_tide / y22_norm / sqrt(6.0_cp) * r_cmb**2 ! Needs a half factor
tide_fac22n = -7.0_cp * tide_fac22p ! Half factor carried over
else
tide_fac20 = 0.0_cp
tide_fac22p = 0.0_cp
tide_fac22n = 0.0_cp
end if
end if
end subroutine preCalc
!-------------------------------------------------------------------------------
subroutine preCalcTimes(time,tEND,dt,n_time_step,n_time_steps)
Expand Down
6 changes: 3 additions & 3 deletions src/rIter.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,12 +50,12 @@ module rIter_mod
& w_Rloc, dw_Rloc, ddw_Rloc, xi_Rloc, omega_ic,&
& omega_ma, phi_Rloc
use time_schemes, only: type_tscheme
use physical_parameters, only: ktops, kbots, n_r_LCR, ktopv, kbotv,&
& ellip_fac_cmb, ellip_fac_icb
use physical_parameters, only: ktops, kbots, n_r_LCR, ktopv, kbotv
use rIteration, only: rIter_t
use RMS, only: get_nl_RMS, transform_to_lm_RMS, compute_lm_forces, &
& transform_to_grid_RMS
use probe_mod
use special, only: ellip_fac_cmb, ellip_fac_icb, amp_tide, l_radial_flow_bc

implicit none

Expand Down Expand Up @@ -586,7 +586,7 @@ subroutine transform_to_grid_space(this, nR, nBc, lViscBcCalc, lRmsCalc, &
& this%gsa%dvtdpc, this%gsa%dvpdpc, l_R(nR))
end if
else if ( nBc == 2 ) then
if ( nR == n_r_cmb .and. ellip_fac_cmb == 0.0_cp ) then
if ( nR == n_r_cmb .and. (.not. l_radial_flow_bc) ) then
call v_rigid_boundary(nR, omega_ma, lDeriv, this%gsa%vrc, &
& this%gsa%vtc, this%gsa%vpc, this%gsa%cvrc, &
& this%gsa%dvrdtc, this%gsa%dvrdpc, &
Expand Down
13 changes: 13 additions & 0 deletions src/special.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,4 +38,17 @@ module special
real(cp), public :: omega_mode_ic, omega_mode_ma ! Frequency of forcing
integer, public :: mode_symm_ic, mode_symm_ma ! Symmetry of forcing: 1 : eq symm, 0 : eq antisymm

real(cp), public :: ellipticity_cmb ! Ellipticity of CMB, used for libration
real(cp), public :: ellipticity_icb ! Ellipticity of ICB, used for libration
real(cp), public :: ellip_fac_cmb ! d/d\phi (Y22) * d/dt exp(i\omega t) at CMB
real(cp), public :: ellip_fac_icb ! d/d\phi (Y22) * d/dt exp(i\omega t) at ICB
real(cp), public :: amp_tide ! Amplitude of tide, such that amp_tide = 1 => max ur(2,0) = 1 at boundary
real(cp), public :: omega_tide ! Frequency of tide
real(cp), public :: tide_fac20 ! Pre-factor for tidal forcing for (2,0)
real(cp), public :: tide_fac22p ! Pre-factor for tidal forcing for (2,2)
real(cp), public :: tide_fac22n ! Pre-factor for tidal forcing for (2,-2)
logical, public :: l_radial_flow_bc ! A flag that tells whether radial flow BCs are being used

real(cp), public :: ampForce ! Amplitude of external body force

end module special
4 changes: 2 additions & 2 deletions src/startFields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,9 +14,9 @@ module start_fields
& otemp1, ogrun, dentropy0, dxicond, r_icb
use physical_parameters, only: interior_model, epsS, impS, n_r_LCR, &
& ktopv, kbotv, LFfac, imagcon, ThExpNb, &
& ViscHeatFac, impXi, ampForce
& ViscHeatFac, impXi
use num_param, only: dtMax, alpha
use special, only: lGrenoble
use special, only: lGrenoble, ampForce
use output_data, only: log_file, n_log_file
use blocking, only: lo_map, llm, ulm, ulmMag, llmMag
use logic, only: l_conv, l_mag, l_cond_ic, l_heat, l_SRMA, l_SRIC, &
Expand Down
79 changes: 65 additions & 14 deletions src/updateWP.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,11 @@ module updateWP_mod
& alpha0, temp0, beta, dbeta, ogrun, l_R, &
& rscheme_oc, ddLvisc, ddbeta, orho1
use physical_parameters, only: kbotv, ktopv, ra, BuoFac, ChemFac, &
& ViscHeatFac, ThExpNb, ktopp, &
& ellipticity_cmb, ellipticity_icb, &
& ellip_fac_cmb, ellip_fac_icb
& ViscHeatFac, ThExpNb, ktopp
use special, only: ellipticity_cmb, ellipticity_icb, &
& ellip_fac_cmb, ellip_fac_icb, &
& tide_fac20, tide_fac22p, tide_fac22n, &
& omega_tide, amp_tide, l_radial_flow_bc
use num_param, only: dct_counter, solve_counter
use init_fields, only: omegaOsz_ma1, tShift_ma1, omegaOsz_ic1, tShift_ic1
use blocking, only: lo_sub_map, lo_map, st_sub_map, llm, ulm, st_map
Expand Down Expand Up @@ -413,19 +415,45 @@ subroutine updateWP(time, s, xi, w, dw, ddw, dwdt, p, dp, dpdt, tscheme, &
rhs1(n_r_max,2*lm-1,threadid)=0.0_cp
rhs1(n_r_max,2*lm,threadid) =0.0_cp

if ( l1 == 2 .and. m1 == 2 ) then
if ( ellipticity_cmb /= 0.0_cp ) then
rhs1(1,2*lm-1,threadid)=ellip_fac_cmb/real(l1*(l1+1),kind=cp) &
& *cos(omegaOsz_ma1*(time+tShift_ma1))
rhs1(1,2*lm,threadid) =ellip_fac_cmb/real(l1*(l1+1),kind=cp) &
& *sin(omegaOsz_ma1*(time+tShift_ma1))
if (l_radial_flow_bc) then

if ( l1 == 2 .and. m1 == 0 ) then
if ( amp_tide /= 0.0_cp ) then

rhs1(1,2*lm-1,threadid)=rhs1(1,2*lm-1,threadid) &
& + tide_fac20/real(l1*(l1+1),kind=cp) &
& * cos(omega_tide*(time))
rhs1(1,2*lm,threadid) =rhs1(1,2*lm,threadid) &
& + tide_fac20/real(l1*(l1+1),kind=cp) &
& * sin(omega_tide*(time))
end if
end if

if ( ellipticity_icb /= 0.0_cp ) then
rhs1(n_r_max,2*lm-1,threadid)=ellip_fac_icb/real(l1*(l1+1),kind=cp) &
& *cos(omegaOsz_ic1*(time+tShift_ic1))
rhs1(n_r_max,2*lm,threadid) =ellip_fac_icb/real(l1*(l1+1),kind=cp) &
& *sin(omegaOsz_ic1*(time+tShift_ic1))
if ( l1 == 2 .and. m1 == 2 ) then
if ( ellipticity_cmb /= 0.0_cp ) then
rhs1(1,2*lm-1,threadid)=ellip_fac_cmb/real(l1*(l1+1),kind=cp) &
& *cos(omegaOsz_ma1*(time+tShift_ma1))
rhs1(1,2*lm,threadid) =ellip_fac_cmb/real(l1*(l1+1),kind=cp) &
& *sin(omegaOsz_ma1*(time+tShift_ma1))
end if

if ( ellipticity_icb /= 0.0_cp ) then
rhs1(n_r_max,2*lm-1,threadid)=ellip_fac_icb/real(l1*(l1+1),kind=cp) &
& *cos(omegaOsz_ic1*(time+tShift_ic1))
rhs1(n_r_max,2*lm,threadid) =ellip_fac_icb/real(l1*(l1+1),kind=cp) &
& *sin(omegaOsz_ic1*(time+tShift_ic1))
end if

if ( amp_tide /= 0.0_cp ) then
rhs1(1,2*lm-1,threadid)=rhs1(1,2*lm-1,threadid) &
& + (tide_fac22p + tide_fac22n) &
& /real(l1*(l1+1),kind=cp) &
& * cos(omega_tide*(time))
rhs1(1,2*lm,threadid) =rhs1(1,2*lm,threadid) &
& + (tide_fac22p - tide_fac22n) &
& /real(l1*(l1+1),kind=cp) &
& * sin(omega_tide*(time))
end if
end if
end if

Expand All @@ -434,6 +462,29 @@ subroutine updateWP(time, s, xi, w, dw, ddw, dwdt, p, dp, dpdt, tscheme, &
rhs1(2,2*lm,threadid) =0.0_cp
rhs1(n_r_max-1,2*lm-1,threadid)=0.0_cp
rhs1(n_r_max-1,2*lm,threadid) =0.0_cp

! Special BC for free-slip and radial flow

if (l_radial_flow_bc .and. (ktopv == 1 .or. kbotv == 1)) then
if (l1 == 2 .and. m1 == 2) then
if ( (ellipticity_cmb /= 0.0_cp .or. amp_tide /= 0.0_cp) &
& .and. ktopv == 1 ) then
rhs1(2,2*lm-1,threadid) = - real(l1*(l1+1),kind=cp)*or2(1)
rhs1(2,2*lm,threadid) = - real(l1*(l1+1),kind=cp)*or2(1)
else if (ellipticity_icb /= 0.0_cp .and. kbotv == 1) then
rhs1(n_r_max-1,2*lm-1,threadid) = - real(l1*(l1+1),kind=cp)*or2(n_r_max)
rhs1(n_r_max-1,2*lm,threadid) = - real(l1*(l1+1),kind=cp)*or2(n_r_max)
end if
end if

if (l1 == 2 .and. m1 == 0) then
if (amp_tide /= 0.0_cp .and. ktopv == 1) then
rhs1(2,2*lm-1,threadid) = - real(l1*(l1+1),kind=cp)*or2(1)
end if
end if
end if


do nR=3,n_r_max-2
rhs1(nR,2*lm-1,threadid)= real(work_LMloc(lm1,nR))
rhs1(nR,2*lm,threadid) =aimag(work_LMloc(lm1,nR))
Expand Down
6 changes: 3 additions & 3 deletions src/updateWPS.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,9 @@ module updateWPS_mod
& ogrun, kappa, orho1, dentropy0, temp0, l_R
use physical_parameters, only: kbotv, ktopv, ktops, kbots, ra, opr, &
& ViscHeatFac, ThExpNb, BuoFac, &
& CorFac, ktopp, ellipticity_cmb, &
& ellipticity_icb, ellip_fac_cmb, &
& ellip_fac_icb
& CorFac, ktopp
use special, only: ellipticity_cmb, ellipticity_icb, ellip_fac_cmb, &
& ellip_fac_icb
use num_param, only: dct_counter, solve_counter
use init_fields, only: tops, bots, omegaOsz_ma1, tShift_ma1, &
& omegaOsz_ic1, tShift_ic1
Expand Down
2 changes: 1 addition & 1 deletion src/updateZ.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module updateZ_mod
use radial_data, only: n_r_cmb, n_r_icb, nRstart, nRstop
use radial_functions, only: visc, or1, or2, rscheme_oc, dLvisc, beta, &
& rho0, r_icb, r_cmb, r, beta, dbeta
use physical_parameters, only: kbotv, ktopv, prec_angle, po, oek, ampForce, &
use physical_parameters, only: kbotv, ktopv, prec_angle, po, oek, &
& gammatau_gravi
use num_param, only: AMstart, dct_counter, solve_counter
use torsional_oscillations, only: ddzASL
Expand Down

0 comments on commit e9bb616

Please sign in to comment.