Skip to content

Commit

Permalink
some cleanups/typo/vectorizations
Browse files Browse the repository at this point in the history
  • Loading branch information
@tgastine
tgastine committed Jul 30, 2024
1 parent 2c46270 commit aa90cd3
Show file tree
Hide file tree
Showing 15 changed files with 94 additions and 121 deletions.
2 changes: 1 addition & 1 deletion src/fields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ subroutine initialize_fields
z_Rloc(:,:)=zero
s_Rloc(:,:)=zero
if ( l_mag ) then
if( l_mag_par_solve ) then
if ( l_mag_par_solve ) then
allocate(b_Rloc(lm_max,nRstart:nRstop), aj_Rloc(lm_max,nRstart:nRstop))
b_Rloc(:,:) =zero
aj_Rloc(:,:)=zero
Expand Down
4 changes: 2 additions & 2 deletions src/init_fields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -541,7 +541,7 @@ subroutine initS(s,p)
call abortRun('Stop run in init')
end if
lm=lo_map%lm2(l,m)
if( (lm>=llm) .and. (lm<=ulm) ) then
if ( (lm>=llm) .and. (lm<=ulm) ) then
do n_r=1,n_r_max
c_r=s1(n_r)*amp_s1
s(lm,n_r)=s(lm,n_r)+cmplx(c_r,0.0_cp,kind=cp)
Expand Down Expand Up @@ -569,7 +569,7 @@ subroutine initS(s,p)
lm=lo_map%lm2(l,m)
s_r=amp_s2
s_i=0.0_cp
if( (lm>=llm) .and. (lm<=ulm) ) then
if ( (lm>=llm) .and. (lm<=ulm) ) then
if ( amp_s2 < 0.0_cp .and. m /= 0 ) then
!-------- Sin(phi)-mode initialized for amp_s2<0
s_r = 0.0_cp
Expand Down
8 changes: 4 additions & 4 deletions src/integration.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,8 +241,8 @@ subroutine cylmean_otc(a,v,n_s_max,n_s_otc,r,s,theta,zDensIn)
exit rbracket
end if
end do rbracket
if(n_r2 == n_r_max-1) n_r2=n_r_max-2
if(n_r2 == 1 ) n_r2=2
if (n_r2 == n_r_max-1) n_r2=n_r_max-2
if (n_r2 == 1 ) n_r2=2
n_r3=n_r2-1
n_r1=n_r2+1
n_r0=n_r2+2
Expand All @@ -253,7 +253,7 @@ subroutine cylmean_otc(a,v,n_s_max,n_s_otc,r,s,theta,zDensIn)
else
n_th1=n_theta_max
tbracket: do n_th=n_theta_max,1,-1
if( theta(n_th) <= thet) then
if ( theta(n_th) <= thet) then
n_th1=n_th
exit tbracket
end if
Expand Down Expand Up @@ -438,7 +438,7 @@ subroutine cylmean_itc(a,vn,vs,n_s_max,n_s_otc,r,s,theta, zDensIn)
else
n_th1=n_theta_max
tbracket: do n_th=n_theta_max,1,-1
if( theta(n_th) <= thet) then
if ( theta(n_th) <= thet) then
n_th1=n_th
exit tbracket
end if
Expand Down
2 changes: 1 addition & 1 deletion src/movie.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -532,7 +532,7 @@ subroutine get_movie_type
n_type=4 ! Horizontal field
file_name='Bh_'
lIC=.true.
else if( index(string,'BS') /= 0) then
else if ( index(string,'BS') /= 0) then
n_type=114
typeStr='cyl radial magnetic field'
file_name='Bs_'
Expand Down
114 changes: 49 additions & 65 deletions src/nonlinear_bcs.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,12 @@ module nonlinear_bcs
use iso_fortran_env, only: output_unit
use precision_mod
use blocking, only: lm2
use truncation, only: n_phi_max, l_max, n_theta_max, nlat_padded, lm_max
use truncation, only: n_phi_max, l_max, nlat_padded, lm_max
use radial_data, only: n_r_cmb, n_r_icb
use radial_functions, only: r_cmb, r_icb, rho0
use physical_parameters, only: sigma_ratio, conductance_ma, prmag, oek
use horizontal_data, only: cosTheta, sinTheta_E2, phi, sinTheta
use constants, only: two, y10_norm, y11_norm
use constants, only: two, y10_norm, y11_norm, zero
use useful, only: abortRun
use sht, only: spat_to_sphertor

Expand Down Expand Up @@ -47,74 +47,65 @@ subroutine get_br_v_bcs(br,vt,vp,omega,O_r_E_2,O_rho,br_vt_lm,br_vp_lm)
real(cp), intent(in) :: O_rho ! :math:`1/\tilde{\rho}` (anelastic)

!-- Output variables:
! br*vt/(sin(theta)**2*r**2)
complex(cp), intent(inout) :: br_vt_lm(lm_max)
! br*(vp/(sin(theta)**2*r**2)-omega_ma)
complex(cp), intent(inout) :: br_vp_lm(lm_max)
complex(cp), intent(inout) :: br_vt_lm(:) ! br*vt/(sin(theta)**2*r**2)
complex(cp), intent(inout) :: br_vp_lm(:) ! br*(vp/(sin(theta)**2*r**2)-omega_ma)

!-- Local variables:
integer :: n_theta, n_phi
integer :: n_phi
real(cp) :: br_vt(nlat_padded,n_phi_max), br_vp(nlat_padded,n_phi_max)
real(cp) :: fac ! 1/( r**2 sin(theta)**2 )

fac=O_r_E_2*O_rho
!$omp parallel do default(shared) &
!$omp& private(n_theta,n_phi)
!$omp parallel do default(shared)
do n_phi=1,n_phi_max
do n_theta=1,n_theta_max
br_vt(n_theta,n_phi)= fac*br(n_theta,n_phi)*vt(n_theta,n_phi)
br_vp(n_theta,n_phi)= br(n_theta,n_phi)* ( fac*vp(n_theta,n_phi)- &
& omega*sinTheta_E2(n_theta) )
end do
br_vt(:,n_phi)= fac*br(:,n_phi)*vt(:,n_phi)
br_vp(:,n_phi)= br(:,n_phi)* ( fac*vp(:,n_phi)-omega*sinTheta_E2(:) )
end do
!$omp end parallel do

call spat_to_sphertor(br_vt, br_vp, br_vt_lm, br_vp_lm, l_max)

end subroutine get_br_v_bcs
!----------------------------------------------------------------------------
subroutine get_b_nl_bcs(bc,br_vt_lm,br_vp_lm,lm_min_b,lm_max_b,b_nl_bc,aj_nl_bc)
subroutine get_b_nl_bcs(bc,br_vt_lm,br_vp_lm,b_nl_bc,aj_nl_bc)
!
! Purpose of this subroutine is to calculate the nonlinear term
! of the magnetic boundary condition for a conducting mantle in
! physical space (theta,phi), assuming that the conductance
! of the mantle is much smaller than that of the core.
! Calculation is performed for the theta block:
!
! .. code-block:: fortran
!
! n_theta_min<=n_theta<=n_theta_min+n_theta_block-1
!

!-- Input variables:
character(len=3), intent(in) :: bc ! Distinguishes 'CMB' and 'ICB'
integer, intent(in) :: lm_min_b,lm_max_b ! limits of lm-block
complex(cp), intent(in) :: br_vt_lm(lm_max) ! :math:`B_r u_\theta/(r^2\sin^2\theta)`
complex(cp), intent(in) :: br_vp_lm(lm_max) ! :math:`B_r u_\phi/(r^2\sin^2\theta)`
complex(cp), intent(in) :: br_vt_lm(:) ! :math:`B_r u_\theta/(r^2\sin^2\theta)`
complex(cp), intent(in) :: br_vp_lm(:) ! :math:`B_r u_\phi/(r^2\sin^2\theta)`

!-- Output variables:
complex(cp), intent(out) :: b_nl_bc(lm_min_b:lm_max_b) ! nonlinear bc for b
complex(cp), intent(out) :: aj_nl_bc(lm_min_b:lm_max_b) ! nonlinear bc for aj
complex(cp), intent(out) :: b_nl_bc(:) ! nonlinear bc for b
complex(cp), intent(out) :: aj_nl_bc(:) ! nonlinear bc for aj

!-- Local variables:
integer :: lm ! position of degree and order
real(cp) :: fac

if ( bc == 'CMB' ) then

b_nl_bc(1) =zero
aj_nl_bc(1)=zero
fac=conductance_ma*prmag
!$omp parallel do default(shared)
do lm=lm_min_b,lm_max_b
do lm=2,lm_max
b_nl_bc(lm) =-fac * br_vt_lm(lm)
aj_nl_bc(lm)=-fac * br_vp_lm(lm)
end do
!$omp end parallel do

else if ( bc == 'ICB' ) then

aj_nl_bc(1)=zero
fac=sigma_ratio*prmag
!$omp parallel do default(shared) private(lm)
do lm=lm_min_b,lm_max_b
do lm=2,lm_max
aj_nl_bc(lm)=-fac * br_vp_lm(lm)
end do
!$omp end parallel do
Expand Down Expand Up @@ -149,7 +140,7 @@ subroutine v_rigid_boundary(nR,omega,lDeriv,vrr,vtr,vpr,cvrr,dvrdtr, &

!-- Local variables:
real(cp) :: r2
integer :: nTheta,nPhi
integer :: nPhi

if ( nR == n_r_cmb ) then
r2=r_cmb*r_cmb
Expand All @@ -162,20 +153,18 @@ subroutine v_rigid_boundary(nR,omega,lDeriv,vrr,vtr,vpr,cvrr,dvrdtr, &
return
end if

!$omp parallel do default(shared) private(nPhi,nTheta)
!$omp parallel do default(shared)
do nPhi=1,n_phi_max
do nTheta=1,n_theta_max
vrr(nTheta,nPhi)=0.0_cp
vtr(nTheta,nPhi)=0.0_cp
vpr(nTheta,nPhi)=r2*rho0(nR)*sinTheta_E2(nTheta)*omega
if ( lDeriv ) then
cvrr(nTheta,nPhi) =r2*rho0(nR)*two*cosTheta(nTheta)*omega
dvrdtr(nTheta,nPhi)=0.0_cp
dvrdpr(nTheta,nPhi)=0.0_cp
dvtdpr(nTheta,nPhi)=0.0_cp
dvpdpr(nTheta,nPhi)=0.0_cp
end if
end do
vrr(:,nPhi)=0.0_cp
vtr(:,nPhi)=0.0_cp
vpr(:,nPhi)=r2*rho0(nR)*sinTheta_E2(:)*omega
if ( lDeriv ) then
cvrr(:,nPhi) =r2*rho0(nR)*two*cosTheta(:)*omega
dvrdtr(:,nPhi)=0.0_cp
dvrdpr(:,nPhi)=0.0_cp
dvtdpr(:,nPhi)=0.0_cp
dvpdpr(:,nPhi)=0.0_cp
end if
end do
!$omp end parallel do

Expand All @@ -190,45 +179,40 @@ subroutine v_center_sphere(ddw, vrr, vtr, vpr)
!

!-- Input variable
complex(cp), intent(in) :: ddw(lm_max)
complex(cp), intent(in) :: ddw(:)

!-- Output variables:
real(cp), intent(out) :: vrr(:,:), vtr(:,:), vpr(:,:)

!-- Local variables:
integer :: nTheta, nPhi, lm10, lm11
integer :: nPhi, lm10, lm11

lm10=lm2(1,0)
lm11=lm2(1,1)

if ( lm11 > 0 ) then ! minc = 1
!$omp parallel do default(shared) private(nPhi,nTheta)
!$omp parallel do default(shared)
do nPhi=1,n_phi_max
do nTheta=1,n_theta_max
vrr(nTheta,nPhi)=y10_norm*real(ddw(lm10))*cosTheta(nTheta) +&
& two*y11_norm*sinTheta(nTheta)*( &
& real(ddw(lm11))*cos(phi(nPhi))- &
& aimag(ddw(lm11))*sin(phi(nPhi)) )
vtr(nTheta,nPhi)=sinTheta(nTheta)*( &
& -y10_norm*real(ddw(lm10))*sinTheta(nTheta) +&
& two*y11_norm*cosTheta(nTheta)*( &
& real(ddw(lm11))*cos(phi(nPhi))- &
& aimag(ddw(lm11))*sin(phi(nPhi)) ) )
vpr(nTheta,nPhi)=-two*y11_norm*sinTheta(nTheta)*( &
& real(ddw(lm11))*sin(phi(nPhi))+ &
& aimag(ddw(lm11))*cos(phi(nPhi)) )
end do
vrr(:,nPhi)=y10_norm*real(ddw(lm10))*cosTheta(:) + &
& two*y11_norm*sinTheta(:)*( &
& real(ddw(lm11))*cos(phi(nPhi))- &
& aimag(ddw(lm11))*sin(phi(nPhi)) )
vtr(:,nPhi)=sinTheta(:)*( &
& -y10_norm*real(ddw(lm10))*sinTheta(:) + &
& two*y11_norm*cosTheta(:)*( &
& real(ddw(lm11))*cos(phi(nPhi))- &
& aimag(ddw(lm11))*sin(phi(nPhi)) ) )
vpr(:,nPhi)=-two*y11_norm*sinTheta(:)*( &
& real(ddw(lm11))*sin(phi(nPhi))+ &
& aimag(ddw(lm11))*cos(phi(nPhi)) )
end do
!$omp end parallel do
else ! minc /= 1
!$omp parallel do default(shared) private(nPhi,nTheta)
!$omp parallel do default(shared)
do nPhi=1,n_phi_max
do nTheta=1,n_theta_max
vrr(nTheta,nPhi)=y10_norm*real(ddw(lm10))*cosTheta(nTheta)
vtr(nTheta,nPhi)=-sinTheta(nTheta)*y10_norm*real(ddw(lm10))* &
& sinTheta(nTheta)
vpr(nTheta,nPhi)=0.0_cp
end do
vrr(:,nPhi)=y10_norm*real(ddw(lm10))*cosTheta(:)
vtr(:,nPhi)=-sinTheta(:)*y10_norm*real(ddw(lm10))*sinTheta(:)
vpr(:,nPhi)=0.0_cp
end do
!$omp end parallel do
end if
Expand Down
28 changes: 13 additions & 15 deletions src/outPar.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -252,21 +252,19 @@ subroutine outPar(s, ds, xi, p, dp, timePassed, timeNorm, l_stop_time, ekinR, &
end if

if ( rank == 0 ) then
do nR=1,n_r_max
! Re must be independant of the timescale
ReR(nR)=sqrt(two*ekinR(nR)*or2(nR)/(4*pi*mass)/eScale)
RoR(nR)=ReR(nR)*ekScaled
if ( dlVR(nR) /= 0.0_cp ) then
RolR(nR)=RoR(nR)/dlVR(nR)
else
RolR(nR)=RoR(nR)
end if
if ( l_mag_nl ) then
RmR(nR)=ReR(nR)*prmag*sigma(nR)*r(nR)*r(nR)
else
RmR(nR)=ReR(nR)*r(nR)*r(nR)
end if
end do
ReR(:)=sqrt(two*ekinR(:)*or2(:)/(4*pi*mass)/eScale)
RoR(:)=ReR(:)*ekScaled
where ( dlVr /= 0.0_cp )
RolR=RoR/dlVR
else where
RolR=RoR
end where

if ( l_mag_nl ) then
RmR(:)=ReR(:)*prmag*sigma(:)*r(:)*r(:)
else
RmR(:)=ReR(:)*r(:)*r(:)
end if

call dlV%compute(dlVR, n_calls, timePassed, timeNorm)
call dlVc%compute(dlVRc, n_calls, timePassed, timeNorm)
Expand Down
2 changes: 1 addition & 1 deletion src/out_TO.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -679,7 +679,7 @@ subroutine interp_theta(a,ac,rr,cyl,theta)
!-- Find indices of angular grid levels that bracket thet
n_th1=n_theta_max
tbracket: do n_th=n_theta_max,1,-1
if( theta(n_th) <= thet) then
if ( theta(n_th) <= thet) then
n_th1=n_th
exit tbracket
end if
Expand Down
12 changes: 6 additions & 6 deletions src/plms.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,9 @@ subroutine plm_theta(theta,max_degree,min_order,max_order,m0, &
end do

plm=sqrt(fac)
if( sin(theta) /= 0.0_cp ) then
if ( sin(theta) /= 0.0_cp ) then
plm=plm*sin(theta)**m
else if( m /= 0 ) then
else if ( m /= 0 ) then
plm=0.0_cp
end if

Expand Down Expand Up @@ -125,13 +125,13 @@ subroutine plm_theta(theta,max_degree,min_order,max_order,m0, &
l=m
pos=pos+1
if ( m < max_degree ) then
if( norm == 0 .OR. norm == 2 ) then
if ( norm == 0 .OR. norm == 2 ) then
dtheta_plma(pos)= l/sqrt(real(2*l+3,cp)) * plma(pos+1)
else if ( norm == 1 ) then
dtheta_plma(pos)= l/sqrt(real(2*l+1,cp)) * plma(pos+1)
end if
else
if( norm == 0 .OR. norm == 2 ) then
if ( norm == 0 .OR. norm == 2 ) then
dtheta_plma(pos)= l/sqrt(real(2*l+3,cp)) * dtheta_plma(pos)
else if ( norm == 1 ) then
dtheta_plma(pos)= l/sqrt(real(2*l+1,cp)) * dtheta_plma(pos)
Expand All @@ -142,7 +142,7 @@ subroutine plm_theta(theta,max_degree,min_order,max_order,m0, &
do l=m+1,max_degree-1

pos=pos+1
if( norm == 0 .OR. norm == 2 ) then
if ( norm == 0 .OR. norm == 2 ) then
dtheta_plma(pos)= &
& l*sqrt( real((l+m+1)*(l-m+1),cp) / &
& real((2*l+1)*(2*l+3),cp) &
Expand All @@ -166,7 +166,7 @@ subroutine plm_theta(theta,max_degree,min_order,max_order,m0, &
if ( m < max_degree ) then
l=max_degree
pos=pos+1
if( norm == 0 .OR. norm == 2 ) then
if ( norm == 0 .OR. norm == 2 ) then
dtheta_plma(pos)= &
& l*sqrt( real((l+m+1)*(l-m+1),cp) / &
& real((2*l+1)*(2*l+3),cp) &
Expand Down
8 changes: 3 additions & 5 deletions src/preCalculations.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -302,12 +302,10 @@ subroutine preCalc(tscheme)
end if

!-- Calculate auxiliary arrays containing effective Courant grid intervals:
delxh2(1) =r_cmb**2/real(l_R(1)*(l_R(1)+1),kind=cp)
delxh2(n_r_max)=r_icb**2/real(l_R(n_r_max)*(l_R(n_r_max)+1),kind=cp)
delxh2(:)=r(:)**2/real(l_R(:)*(l_R(:)+1),kind=cp)
delxr2(1) =(r(1)-r(2))**2
delxr2(n_r_max)=(r(n_r_max-1)-r(n_r_max))**2
do n_r=2,n_r_max-1
delxh2(n_r)=r(n_r)**2/real(l_R(n_r)*(l_R(n_r)+1),kind=cp)
delmin=min((r(n_r-1)-r(n_r)),(r(n_r)-r(n_r+1)))
delxr2(n_r)=delmin*delmin
end do
Expand Down Expand Up @@ -743,7 +741,7 @@ subroutine preCalc(tscheme)
do l=1,l_max
fac_loop(l)=0.0_cp
if (mod(l,2)/=0) then
if(l==1) then
if ( l==1 ) then
fac_loop(l)= one
else
fac_loop(l)= -fac_loop(l-2)*loopRadRatio**2*real(l,kind=cp)/ &
Expand All @@ -752,7 +750,7 @@ subroutine preCalc(tscheme)
end if
end do

if (l_non_rot) then
if ( l_non_rot ) then
amp_curr = Le
else
amp_curr = Le * sqrt(prmag/ek)
Expand Down
Loading

0 comments on commit aa90cd3

Please sign in to comment.