From abbdba5abf2195757d78ddcbb6b360f687b0462c Mon Sep 17 00:00:00 2001 From: Thomas Gastine Date: Wed, 12 Feb 2025 11:03:14 +0100 Subject: [PATCH 1/2] fix heat/comp diffusion at r=0 in case of full sphere MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Make use of L'Hôpital's rule to replace (2/r)*d/dr by 2*d^2/dr^2. --- src/updatePHI.f90 | 17 ++++++++++++++++- src/updateS.f90 | 44 +++++++++++++++++++++++++++++++++++++++++--- src/updateXI.f90 | 21 ++++++++++++++++++++- 3 files changed, 77 insertions(+), 5 deletions(-) diff --git a/src/updatePHI.f90 b/src/updatePHI.f90 index 69dece3f..e81ecc7b 100644 --- a/src/updatePHI.f90 +++ b/src/updatePHI.f90 @@ -17,7 +17,7 @@ module updatePhi_mod use logic, only: l_finite_diff, l_full_sphere, l_parallel_solve use parallel_mod, only: rank, chunksize, n_procs, get_openmp_blocks use radial_der, only: get_ddr, get_ddr_ghost, exch_ghosts, bulk_to_ghost - use constants, only: zero, one, two + use constants, only: zero, one, two, three use fields, only: work_LMloc use mem_alloc, only: bytes_allocated use useful, only: abortRun @@ -594,6 +594,9 @@ subroutine get_phase_rhs_imp_ghost(phig, dphidt, istage, l_calc_lin) dphidt%impl(lm,n_r,istage)= phaseDiffFac * ( work_Rloc(lm,n_r) + & & two*or1(n_r) * dphi(lm,n_r) - & & dL*or2(n_r) * phig(lm,n_r) ) + if ( l_full_sphere .and. l==0 .and. n_r==n_r_icb ) then + dphidt%impl(lm,n_r,istage)=phaseDiffFac*three*work_Rloc(lm,n_r) + end if end do end do end if @@ -925,6 +928,18 @@ subroutine get_phiMat_Rdist(tscheme,phiMat) phiMat%up(l,nR) =-tscheme%wimp_lin(1)*phaseDiffFac*( & & rscheme_oc%ddr(nR,2) + & & two*or1(nR)* rscheme_oc%dr(nR,2) ) + + if ( l==0 .and. l_full_sphere .and. nR==n_r_icb ) then + !-- Use L'Hôpital's rule to replace 2/r d/dr at the center + !-- by 2*d^2/dr^2 + phiMat%diag(l,nR)= 5.0_cp/6.0_cp*stef*pr- & + & tscheme%wimp_lin(1)*phaseDiffFac* & + & three* rscheme_oc%ddr(nR,1) + phiMat%low(l,nR)=-tscheme%wimp_lin(1)*phaseDiffFac* & + & three* rscheme_oc%ddr(nR,0) + phiMat%up(l,nR) =-tscheme%wimp_lin(1)*phaseDiffFac* & + & three* rscheme_oc%ddr(nR,2) + end if end do end do !$omp end do diff --git a/src/updateS.f90 b/src/updateS.f90 index 5f372641..88d35fab 100644 --- a/src/updateS.f90 +++ b/src/updateS.f90 @@ -23,7 +23,7 @@ module updateS_mod use radial_der, only: get_ddr, get_dr, get_dr_Rloc, get_ddr_ghost, & & exch_ghosts, bulk_to_ghost use fields, only: work_LMloc - use constants, only: zero, one, two + use constants, only: zero, one, two, three use useful, only: abortRun use time_schemes, only: type_tscheme use time_array, only: type_tarray @@ -785,6 +785,11 @@ subroutine get_entropy_rhs_imp_ghost(sg, ds, dsdt, phi, istage, l_calc_lin) & work_Rloc(lm,n_r) & & + ( beta(n_r)+two*or1(n_r)+dLkappa(n_r) ) * ds(lm,n_r) & & - dL*or2(n_r)*sg(lm,n_r) ) + if ( l_full_sphere .and. l==0 .and. n_r==n_r_icb ) then + dsdt%impl(lm,n_r,istage)= opr*hdif_S(l)* kappa(n_r) * ( & + & three * work_Rloc(lm,n_r) & + & + ( beta(n_r)+dLkappa(n_r) ) * ds(lm,n_r) ) + end if end do end do else @@ -797,6 +802,11 @@ subroutine get_entropy_rhs_imp_ghost(sg, ds, dsdt, phi, istage, l_calc_lin) & + ( beta(n_r)+dLtemp0(n_r)+two*or1(n_r)+dLkappa(n_r) ) & & * ds(lm,n_r) & & - dL*or2(n_r) * sg(lm,n_r) ) + if ( l_full_sphere .and. l==0 .and. n_r==n_r_icb ) then + dsdt%impl(lm,n_r,istage)= opr*hdif_S(l)*kappa(n_r) * ( & + & three* work_Rloc(lm,n_r) & + & + ( beta(n_r)+dLtemp0(n_r)+dLkappa(n_r) ) * ds(lm,n_r) ) + end if end do end do end if @@ -1196,7 +1206,6 @@ subroutine get_sMat_Rdist(tscheme,hdif,sMat) sMat%up(l,nR)= -tscheme%wimp_lin(1)*opr*hdif(l)* & & kappa(nR)*( rscheme_oc%ddr(nR,2) + & &( beta(nR)+two*or1(nR)+dLkappa(nR) )*rscheme_oc%dr(nR,2) ) - else sMat%diag(l,nR)=one-tscheme%wimp_lin(1)*opr*hdif(l)* & & kappa(nR)*( rscheme_oc%ddr(nR,1) + & @@ -1212,6 +1221,35 @@ subroutine get_sMat_Rdist(tscheme,hdif,sMat) & ( beta(nR)+dLtemp0(nR)+two*or1(nR)+dLkappa(nR) )* & & rscheme_oc%dr(nR,2) ) end if + + if ( l==0 .and. l_full_sphere .and. nR==n_r_icb ) then + !-- Use L'Hôpital's rule to replace 2/r d/dr at the center + !-- by 2*d^2/dr^2 + if ( l_anelastic_liquid ) then + sMat%diag(l,nR)= one - tscheme%wimp_lin(1)*opr*hdif(l)*& + & kappa(nR)*( three* rscheme_oc%ddr(nR,1) + & + & ( beta(nR)+dLkappa(nR) )*rscheme_oc%dr(nR,1)) + sMat%low(l,nR)= -tscheme%wimp_lin(1)*opr*hdif(l)* & + & kappa(nR)*( three* rscheme_oc%ddr(nR,0) + & + & ( beta(nR)+dLkappa(nR) )*rscheme_oc%dr(nR,0) ) + sMat%up(l,nR)= -tscheme%wimp_lin(1)*opr*hdif(l)* & + & kappa(nR)*( three* rscheme_oc%ddr(nR,2) + & + & ( beta(nR)+dLkappa(nR) )*rscheme_oc%dr(nR,2) ) + else + sMat%diag(l,nR)=one-tscheme%wimp_lin(1)*opr*hdif(l)* & + & kappa(nR)*( three*rscheme_oc%ddr(nR,1) + & + & ( beta(nR)+dLtemp0(nR)+dLkappa(nR) )* & + & rscheme_oc%dr(nR,1) ) + sMat%low(l,nR)= -tscheme%wimp_lin(1)*opr*hdif(l)* & + & kappa(nR)*( three*rscheme_oc%ddr(nR,0) + & + & ( beta(nR)+dLtemp0(nR)+dLkappa(nR) )* & + & rscheme_oc%dr(nR,0) ) + sMat%up(l,nR) = -tscheme%wimp_lin(1)*opr*hdif(l)* & + & kappa(nR)*( three*rscheme_oc%ddr(nR,2) + & + & ( beta(nR)+dLtemp0(nR)+dLkappa(nR) )* & + & rscheme_oc%dr(nR,2) ) + end if + end if end do end do !$omp end do @@ -1257,6 +1295,6 @@ subroutine get_sMat_Rdist(tscheme,hdif,sMat) !-- LU decomposition: call sMat%prepare_mat() - end subroutine get_Smat_Rdist + end subroutine get_sMat_Rdist !----------------------------------------------------------------------------- end module updateS_mod diff --git a/src/updateXI.f90 b/src/updateXI.f90 index 66e01904..38aa4cc2 100644 --- a/src/updateXI.f90 +++ b/src/updateXI.f90 @@ -19,7 +19,7 @@ module updateXi_mod use parallel_mod, only: rank, chunksize, n_procs, get_openmp_blocks use radial_der, only: get_ddr, get_dr, get_dr_Rloc, get_ddr_ghost, exch_ghosts,& & bulk_to_ghost - use constants, only: zero, one, two + use constants, only: zero, one, two, three use fields, only: work_LMloc use mem_alloc, only: bytes_allocated use useful, only: abortRun @@ -704,6 +704,11 @@ subroutine get_comp_rhs_imp_ghost(xig, dxidt, istage, l_calc_lin) dxidt%impl(lm,n_r,istage)= osc*hdif_Xi(l) * & & ( work_Rloc(lm,n_r)+(beta(n_r)+two*or1(n_r)) * dxi(lm,n_r) & & - dL*or2(n_r)* xig(lm,n_r) ) + + if ( l==0 .and. l_full_sphere .and. n_r==n_r_icb ) then + dxidt%impl(lm,n_r,istage)= osc*hdif_Xi(l) * & + & ( three*work_Rloc(lm,n_r)+beta(n_r)*dxi(lm,n_r) ) + end if end do end do end if @@ -1020,6 +1025,20 @@ subroutine get_xiMat_Rdist(tscheme,hdif,xiMat) xiMat%up(l,nR) =-tscheme%wimp_lin(1)*osc*hdif(l)*( & & rscheme_oc%ddr(nR,2) + & & ( beta(nR)+two*or1(nR) )* rscheme_oc%dr(nR,2) ) + + if ( l==0 .and. l_full_sphere .and. nR==n_r_icb ) then + !-- Use L'Hôpital's rule to replace 2/r d/dr at the center + !-- by 2*d^2/dr^2 + xiMat%diag(l,nR)=one-tscheme%wimp_lin(1)*osc*hdif(l)*( & + & three* rscheme_oc%ddr(nR,1) + & + & beta(nR)* rscheme_oc%dr(nR,1) ) + xiMat%low(l,nR)=-tscheme%wimp_lin(1)*osc*hdif(l)*( & + & three* rscheme_oc%ddr(nR,0) + & + & beta(nR)* rscheme_oc%dr(nR,0) ) + xiMat%up(l,nR) =-tscheme%wimp_lin(1)*osc*hdif(l)*( & + & three* rscheme_oc%ddr(nR,2) + & + & beta(nR)* rscheme_oc%dr(nR,2) ) + end if end do end do !$omp end do From 4a1ab24867c373d18acace6360c81c086276ac01 Mon Sep 17 00:00:00 2001 From: Thomas Gastine Date: Wed, 12 Feb 2025 11:05:19 +0100 Subject: [PATCH 2/2] python fix --- python/magic/coeff.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/magic/coeff.py b/python/magic/coeff.py index f2272512..6c71f926 100644 --- a/python/magic/coeff.py +++ b/python/magic/coeff.py @@ -1253,7 +1253,7 @@ def cwt(self, ell, w0=20, nfreq=256, fmin_fac=8, fmax_fac=0.5, logarithmic scale (linear by default) :type logscale: bool """ - assert version > '1.4.0' + #assert version > '1.4.0' dt = np.diff(self.time) # If the data is not regularly sampled, use splines to resample them