From e9bb616d4618bce3a642fa3c5512d6dabac3be70 Mon Sep 17 00:00:00 2001 From: AnkitBarik Date: Sat, 15 Jun 2024 00:47:25 -0400 Subject: [PATCH 1/5] 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 From a23c25739f6f350eab46c236a745bc5338da4554 Mon Sep 17 00:00:00 2001 From: AnkitBarik Date: Mon, 17 Jun 2024 19:06:38 -0400 Subject: [PATCH 2/5] Fix signs + detailed comments --- src/preCalculations.f90 | 7 ++++--- src/updateWP.f90 | 7 ++++++- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/preCalculations.f90 b/src/preCalculations.f90 index c6fbae3b..1d05a0f0 100644 --- a/src/preCalculations.f90 +++ b/src/preCalculations.f90 @@ -783,9 +783,10 @@ subroutine preCalc(tscheme) 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 + tide_fac20 = amp_tide / y20_norm * r_cmb**2 ! (2,0,1) mode of Ogilvie 2014 + tide_fac22p = half * amp_tide / y22_norm / sqrt(6.0_cp) * r_cmb**2 ! Needs a half factor, (2,2,1) mode + tide_fac22n = -7.0_cp * tide_fac22p ! Half factor carried over, (2,2,3) mode, + ! has opposite sign to that of the other two (Polfliet & Smeyers, 1990) else tide_fac20 = 0.0_cp tide_fac22p = 0.0_cp diff --git a/src/updateWP.f90 b/src/updateWP.f90 index d77dff28..acc7c338 100644 --- a/src/updateWP.f90 +++ b/src/updateWP.f90 @@ -445,12 +445,17 @@ subroutine updateWP(time, s, xi, w, dw, ddw, dwdt, p, dp, dpdt, tscheme, & end if if ( amp_tide /= 0.0_cp ) then + ! tide_fac22p -> (2,2,1) + ! tide_fac22n -> (2,2,3) + ! (2,2,3) must have the same signed frequency + ! as (2,0,1) above while (2,2,1) has the opposite + ! sign in the rotating frame (Ogilvie, 2014) 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) & + & + (-tide_fac22p + tide_fac22n) & & /real(l1*(l1+1),kind=cp) & & * sin(omega_tide*(time)) end if From 8d3462df35cda70a65ddf0bf058abd94fda98cab Mon Sep 17 00:00:00 2001 From: AnkitBarik Date: Wed, 24 Jul 2024 13:03:11 -0400 Subject: [PATCH 3/5] Fix for simpson import --- python/magic/bLayers.py | 5 ++++- python/magic/butterfly.py | 5 ++++- python/magic/libmagic.py | 9 ++++++--- python/magic/thHeat.py | 6 +++++- 4 files changed, 19 insertions(+), 6 deletions(-) diff --git a/python/magic/bLayers.py b/python/magic/bLayers.py index 754ba020..c05e0ba6 100644 --- a/python/magic/bLayers.py +++ b/python/magic/bLayers.py @@ -5,7 +5,10 @@ from magic import MagicRadial, matder, intcheb, MagicSetup, scanDir, AvgField from magic.setup import labTex from scipy.signal import argrelextrema -from scipy.integrate import simps +try: + from scipy.integrate import simps +except: + from scipy.integrate import simpson as simps from scipy.interpolate import splrep, splev def getAccuratePeaks(rad, uh, uhTop, uhBot, ri, ro): diff --git a/python/magic/butterfly.py b/python/magic/butterfly.py index 4e220ef5..5931c452 100644 --- a/python/magic/butterfly.py +++ b/python/magic/butterfly.py @@ -8,8 +8,11 @@ import scipy.interpolate as S from magic.plotlib import cut from magic.setup import labTex -from scipy.integrate import simps from .npfile import * +try: + from scipy.integrate import simps +except: + from scipy.integrate import simpson as simps class Butterfly: diff --git a/python/magic/libmagic.py b/python/magic/libmagic.py index 6f9b458c..ac8f4222 100644 --- a/python/magic/libmagic.py +++ b/python/magic/libmagic.py @@ -1,10 +1,13 @@ # -*- coding: utf-8 -*- import scipy.interpolate as S -from scipy.integrate import simps import numpy as np import glob, os, re, sys import subprocess as sp from .npfile import * +try: + from scipy.integrate import simps +except: + from scipy.integrate import simpson as simps def selectField(obj, field, labTex=True, ic=False): @@ -561,7 +564,7 @@ def fd_grid(nr, a, b, fd_stretching=0.3, fd_ratio=0.1): dr_after = dr_before for i in range(1, nr): rr[i]=rr[i-1]-dr_before - + else: n_boundary_points = int( float(nr-1)/(2.*(1+ratio1)) ) ratio1 = float(nr-1)/float(2.*n_boundary_points) -1. @@ -574,7 +577,7 @@ def fd_grid(nr, a, b, fd_stretching=0.3, fd_ratio=0.1): dr_before = dr_before*dr_after dr_before = 1./(float(n_bulk_points)+ \ 2.*dr_after*((1-dr_before)/(1.-dr_after))) - + for i in range(n_boundary_points): dr_before = dr_before*dr_after diff --git a/python/magic/thHeat.py b/python/magic/thHeat.py index 4a3f97ec..eac95662 100644 --- a/python/magic/thHeat.py +++ b/python/magic/thHeat.py @@ -3,7 +3,11 @@ import numpy as np from magic import scanDir, MagicSetup, Movie, chebgrid, rderavg, AvgField import os, pickle -from scipy.integrate import simps, trapz +from scipy.integrate import trapz +try: + from scipy.integrate import simps +except: + from scipy.integrate import simpson as simps json_model = {'phys_params': [], 'time_series': { 'heat': ['topnuss', 'botnuss']}, From e0c8472469ff0e4b1363e09e25015f17bd1e90b7 Mon Sep 17 00:00:00 2001 From: Thomas Gastine Date: Thu, 25 Jul 2024 12:17:42 +0200 Subject: [PATCH 4/5] update deprecated Makefile for backup usage. Still not perfect though --- src/Makefile | 48 ++++++++++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/src/Makefile b/src/Makefile index a9ac5e41..aee0e0f9 100644 --- a/src/Makefile +++ b/src/Makefile @@ -16,7 +16,7 @@ USE_PRECOND=yes # USE_FFTLIB can be MKL, FFTW or JW USE_FFTLIB=MKL #USE_DCTLIB can be MKL, FFTW or JW -USE_DCTLIB=JW +USE_DCTLIB=MKL #USE_LAPACKLIB can be MKL, LAPACK or JW USE_LAPACKLIB=MKL #Use shtns for Legendre/Fourier transforms @@ -295,16 +295,15 @@ endif OBJS := $(addsuffix .o, $(basename $(RED_SOURCES))) OBJS += truncation.o -OBJS += c_utils.o .SUFFIXES: -ifeq ($(COMPILER),gnu) -ifeq ($(USE_MPI),yes) -all: mpi.mod - make $(OUT) -endif -endif +#ifeq ($(COMPILER),gnu) +#ifeq ($(USE_MPI),yes) +#all: mpi.mod +# make $(OUT) +#endif +#endif $(OUT): $(OBJS) $(FC) $(FFLAGS) -o $@ $^ $(LIBS) @@ -363,7 +362,7 @@ get_nl.o : truncation.o horizontal.o logic.o time_schemes.o\ phys_param.o radial.o blocking.o mem_alloc.o out_movie_file.o : truncation.o output_data.o blocking.o\ radial.o horizontal.o movie.o fields.o\ - num_param.o $(FFT_OBJS) logic.o\ + num_param.o $(FFT_OBJS) logic.o outMisc.o\ out_dtB_frame.o parallel.o communications.o storeCheckPoints.o : truncation.o phys_param.o num_param.o logic.o\ output_data.o init_fields.o dt_fieldsLast.o\ @@ -381,7 +380,8 @@ power.o : truncation.o blocking.o horizontal.o\ phys_param.o num_param.o radial.o logic.o\ output_data.o outRot.o integration.o useful.o\ mem_alloc.o mean_sd.o -integration.o : $(DCT_OBJS) precision_mod.o constants.o +integration.o : $(DCT_OBJS) precision_mod.o constants.o chebyshev.o\ + communications.o out_coeff.o : logic.o precision_mod.o parallel.o blocking.o\ truncation.o communications.o\ mem_alloc.o radial.o output_data.o phys_param.o @@ -395,9 +395,9 @@ magic.o : truncation.o num_param.o parallel.o logic.o\ step_time.o : truncation.o phys_param.o updateB.o updateZ.o\ num_param.o radial_data.o updateWP.o time_schemes.o\ logic.o output_data.o output.o movie.o updateS.o\ - timing.o courant.o signals.o updateXI.o\ + timing.o courant.o signals.o updateXI.o updatePHI.o\ radialLoop.o nonlinear_bcs.o updateWPS.o\ - LMLoop.o dt_fieldsLast.o c_utils.o useful.o probes.o + LMLoop.o dt_fieldsLast.o useful.o probes.o phys_param.o : precision_mod.o char_manip.o : precision_mod.o constants.o : precision_mod.o @@ -408,7 +408,7 @@ chebyshev_polynoms.o : constants.o logic.o precision_mod.o num_param.o rIteration.o : precision_mod.o time_schemes.o truncation.o radial_data.o rIter.o : precision_mod.o parallel.o truncation.o logic.o radial_data.o\ radial.o constants.o get_nl.o get_td.o TO.o dtB.o\ - out_graph_file.o nonlinear_bcs.o\ + out_graph_file.o nonlinear_bcs.o out_movie_file.o\ courant.o outRot.o fields.o time_schemes.o phys_param.o\ ${SHT_OBJS} outGeos.o radialLoop.o : truncation.o time_schemes.o precision_mod.o mem_alloc.o\ @@ -420,9 +420,9 @@ LMLoop.o : truncation.o blocking.o parallel.o time_array.o\ updateZ.o updateWP.o debugging.o mpi_transpose.o : precision_mod.o truncation.o radial_data.o\ mem_alloc.o parallel.o blocking.o -debugging.o : precision_mod.o +debugging.o : precision_mod.o precision_mod.o timing.o : parallel.o precision_mod.o -radial_derivatives.o : $(DCT_OBJS) constants.o precision_mod.o +radial_derivatives.o : $(DCT_OBJS) constants.o precision_mod.o chebyshev.o output.o : truncation.o blocking.o phys_param.o\ num_param.o logic.o output_data.o radial.o\ horizontal.o constants.o fields.o radial_spectra.o\ @@ -456,17 +456,20 @@ updateWPS.o : truncation.o blocking.o num_param.o time_array.o\ phys_param.o radial.o horizontal.o logic.o\ $(DCT_OBJS) radial_derivatives.o\ mem_alloc.o init_fields.o RMS.o time_schemes.o +updatePHI.o : truncation.o blocking.o num_param.o time_array.o\ + phys_param.o radial.o horizontal.o logic.o\ + $(DCT_OBJS) radial_derivatives.o init_fields.o\ + mem_alloc.o time_schemes.o get_td.o : truncation.o blocking.o horizontal.o\ phys_param.o num_param.o radial.o logic.o\ - RMS.o RMS_helpers.o fields.o\ - cutils_iface.o mem_alloc.o + RMS.o RMS_helpers.o fields.o mem_alloc.o store_movie_IC.o : truncation.o blocking.o logic.o radial_data.o\ movie.o radial.o horizontal.o\ $(SHT_OBJS) $(FFT_OBJS) out_movie_file.o phys_param.o radial_spectra.o : truncation.o blocking.o horizontal.o radial.o\ num_param.o output_data.o logic.o useful.o\ char_manip.o radial_data.o LMmapping.o constants.o -fft_fac.o : constants.o precision_mod.o +fft_fac.o : constants.o precision_mod.o mem_alloc.o cosine_transform_even.o: truncation.o fft_fac.o constants.o useful.o\ mem_alloc.o readCheckPoints.o : truncation.o blocking.o phys_param.o time_array.o\ @@ -505,7 +508,8 @@ updateB.o : truncation.o blocking.o horizontal.o fields.o\ RMS_helpers.o mem_alloc.o outGeos.o : truncation.o horizontal.o radial.o constants.o\ output_data.o logic.o radial_data.o parallel.o\ - communications.o integration.o mem_alloc.o logic.o + communications.o integration.o mem_alloc.o logic.o\ + movie.o preCalculations.o : truncation.o phys_param.o num_param.o constants.o\ radial.o horizontal.o init_fields.o time_schemes.o\ blocking.o logic.o output_data.o\ @@ -530,7 +534,7 @@ radial.o : truncation.o radial_data.o $(LAPACK_OBJS)\ chebyshev_polynoms.o radial_derivatives.o\ cosine_transform_even.o mem_alloc.o useful.o radial_scheme.o\ chebyshev.o finite_differences.o -output_data.o : precision_mod.o +output_data.o : precision_mod.o mem_alloc.o horizontal.o : truncation.o phys_param.o num_param.o precision_mod.o\ radial.o logic.o blocking.o plms.o $(FFT_OBJS)\ mem_alloc.o @@ -539,7 +543,7 @@ init_fields.o : truncation.o blocking.o radial.o horizontal.o\ constants.o logic.o $(FFT_OBJS) $(SHT_OBJS)\ useful.o phys_param.o mpi_transpose.o\ $(DCT_OBJS) mem_alloc.o parallel.o -mean_sd.o : precision_mod.o mem_alloc.o matrices.o +mean_sd.o : precision_mod.o mem_alloc.o matrices.o communications.o movie.o : truncation.o parallel.o radial_data.o\ output_data.o logic.o radial.o mem_alloc.o\ horizontal.o char_manip.o useful.o @@ -562,7 +566,7 @@ multistep_schemes.o : parallel.o precision_mod.o num_param.o constants.o\ time_array.o dt_fieldsLast.o : truncation.o blocking.o precision_mod.o time_array.o\ logic.o mem_alloc.o constants.o -fields.o : truncation.o blocking.o radial_data.o\ +fields.o : truncation.o blocking.o radial_data.o communications.o\ precision_mod.o mem_alloc.o logic.o special.o : precision_mod.o mem_alloc.o truncation.o probes.o : parallel.o precision_mod.o truncation.o radial.o num_param.o\ From 5587f18f35940233e16188e9fa9044b7cb673691 Mon Sep 17 00:00:00 2001 From: AnkitBarik Date: Thu, 25 Jul 2024 15:26:22 -0400 Subject: [PATCH 5/5] Backward compatible fix for trapz --- python/magic/compSims.py | 6 ++++-- python/magic/surf.py | 7 +++++-- python/magic/thHeat.py | 4 ++-- 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/python/magic/compSims.py b/python/magic/compSims.py index c844efd3..30d70088 100644 --- a/python/magic/compSims.py +++ b/python/magic/compSims.py @@ -1,11 +1,13 @@ # -*- coding: utf-8 -*- from magic import * -from scipy.integrate import trapz from matplotlib.ticker import ScalarFormatter import matplotlib.pyplot as plt import numpy as np import os - +try: + from scipy.integrate import trapz +except: + from scipy.integrate import trapezoid as trapz class CompSims: """ diff --git a/python/magic/surf.py b/python/magic/surf.py index 4e4a128b..0f2ac504 100644 --- a/python/magic/surf.py +++ b/python/magic/surf.py @@ -7,7 +7,10 @@ import matplotlib.pyplot as plt import os import numpy as np -from scipy.integrate import trapz +try: + from scipy.integrate import trapz +except: + from scipy.integrate import trapezoid as trapz class Surf: @@ -869,7 +872,7 @@ def avg(self, field='vphi', levels=defaultLevels, cm=None, vs = self.gr.vr * np.sin(th3D) + self.gr.vtheta * np.cos(th3D) vz = self.gr.vr * np.cos(th3D) - self.gr.vtheta * np.sin(th3D) - vortz = 1./s3D*(-phideravg(vs, self.gr.minc) + + vortz = 1./s3D*(-phideravg(vs, self.gr.minc) + sderavg(s3D*self.gr.vphi, self.gr.radius)) data = vortz * vz diff --git a/python/magic/thHeat.py b/python/magic/thHeat.py index eac95662..1f84e35a 100644 --- a/python/magic/thHeat.py +++ b/python/magic/thHeat.py @@ -3,11 +3,11 @@ import numpy as np from magic import scanDir, MagicSetup, Movie, chebgrid, rderavg, AvgField import os, pickle -from scipy.integrate import trapz try: - from scipy.integrate import simps + from scipy.integrate import simps, trapz except: from scipy.integrate import simpson as simps + from scipy.integrate import trapezoid as trapz json_model = {'phys_params': [], 'time_series': { 'heat': ['topnuss', 'botnuss']},