Skip to content

Merge prep #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ygaillard wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Merge prep #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
abbdba5
fix heat/comp diffusion at r=0 in case of full sphere
tgastine Feb 12, 2025
29c7a73
Merge branch 'master' of ssh://github.com/magic-sph/magic
tgastine Feb 12, 2025
4a1ab24
python fix
tgastine Feb 12, 2025
4189eb3
Merge branch 'magic-sph:master' into merge-prep
ygaillard Feb 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion python/magic/coeff.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
17 changes: 16 additions & 1 deletion src/updatePHI.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
44 changes: 41 additions & 3 deletions src/updateS.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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) + &
Expand All @@ -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
Expand Down Expand Up @@ -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
21 changes: 20 additions & 1 deletion src/updateXI.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down