From e9bb616d4618bce3a642fa3c5512d6dabac3be70 Mon Sep 17 00:00:00 2001 From: AnkitBarik Date: Sat, 15 Jun 2024 00:47:25 -0400 Subject: [PATCH] Eccentricity tides + cleanups --- src/Namelists.f90 | 11 +++++- src/fields.f90 | 2 +- src/init_fields.f90 | 3 +- src/phys_param.f90 | 6 ---- src/preCalculations.f90 | 45 +++++++++++++++-------- src/rIter.f90 | 6 ++-- src/special.f90 | 13 +++++++ src/startFields.f90 | 4 +-- src/updateWP.f90 | 79 +++++++++++++++++++++++++++++++++-------- src/updateWPS.f90 | 6 ++-- src/updateZ.f90 | 2 +- 11 files changed, 131 insertions(+), 46 deletions(-) diff --git a/src/Namelists.f90 b/src/Namelists.f90 index 143c4bb7..6485897f 100644 --- a/src/Namelists.f90 +++ b/src/Namelists.f90 @@ -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, & @@ -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 @@ -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" @@ -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 diff --git a/src/fields.f90 b/src/fields.f90 index 1d499ae8..b39a29b4 100644 --- a/src/fields.f90 +++ b/src/fields.f90 @@ -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, & diff --git a/src/init_fields.f90 b/src/init_fields.f90 index 7e284e8f..c15e8066 100644 --- a/src/init_fields.f90 +++ b/src/init_fields.f90 @@ -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 diff --git a/src/phys_param.f90 b/src/phys_param.f90 index 03cfdc74..3f027e62 100644 --- a/src/phys_param.f90 +++ b/src/phys_param.f90 @@ -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 diff --git a/src/preCalculations.f90 b/src/preCalculations.f90 index c342ca93..c6fbae3b 100644 --- a/src/preCalculations.f90 +++ b/src/preCalculations.f90 @@ -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 @@ -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 @@ -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) diff --git a/src/rIter.f90 b/src/rIter.f90 index 5494267f..4195d9ca 100644 --- a/src/rIter.f90 +++ b/src/rIter.f90 @@ -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 @@ -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, & diff --git a/src/special.f90 b/src/special.f90 index d7d5076d..eacbb7fa 100644 --- a/src/special.f90 +++ b/src/special.f90 @@ -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 diff --git a/src/startFields.f90 b/src/startFields.f90 index 1cc9e558..3c4e8b0a 100644 --- a/src/startFields.f90 +++ b/src/startFields.f90 @@ -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, & diff --git a/src/updateWP.f90 b/src/updateWP.f90 index fd693c83..d77dff28 100644 --- a/src/updateWP.f90 +++ b/src/updateWP.f90 @@ -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 @@ -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 @@ -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)) diff --git a/src/updateWPS.f90 b/src/updateWPS.f90 index f0719da6..e0f93840 100644 --- a/src/updateWPS.f90 +++ b/src/updateWPS.f90 @@ -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 diff --git a/src/updateZ.f90 b/src/updateZ.f90 index 34be2d19..8f44f712 100644 --- a/src/updateZ.f90 +++ b/src/updateZ.f90 @@ -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