Merge pull request #441 from danieljprice/omp-radiation

Omp radiation
danieljprice · Jul 7, 2023 · b1e1cf0 · b1e1cf0
2 parents 8d91ab0 + 7155196
commit b1e1cf0
Showing 1 changed file with 31 additions and 27 deletions.
diff --git a/src/main/radiation_implicit.f90 b/src/main/radiation_implicit.f90
@@ -79,7 +79,7 @@ subroutine do_radiation_implicit(dt,npart,rad,xyzh,vxyzu,radprop,drad,ierr)
     moresweep = .false.
     dtsub = dt/nsubsteps
     over_substeps: do i = 1,nsubsteps
-       call do_radiation_onestep(dtsub,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit,errorE,errorU,moresweep,ierr)
+       call do_radiation_onestep(dtsub,npart,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit,errorE,errorU,moresweep,ierr)
        if (failed .or. moresweep) then
           ierr = ierr_failed_to_converge
           call warning('radiation_implicit','integration failed - using U and E values anyway')
@@ -136,18 +136,19 @@ end subroutine save_radiation_energies
 !  perform single iteration
 !+
 !---------------------------------------------------------
-subroutine do_radiation_onestep(dt,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit,errorE,errorU,moresweep,ierr)
+subroutine do_radiation_onestep(dt,npart,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit,errorE,errorU,moresweep,ierr)
  use io,      only:fatal,error,iverbose,warning
  use part,    only:hfact
  use physcon, only:pi
  use kernel,  only:radkern
  real, intent(in)     :: dt,xyzh(:,:),origEU(:,:)
+ integer, intent(in)  :: npart
  real, intent(inout)  :: radprop(:,:),rad(:,:),vxyzu(:,:)
  logical, intent(out) :: failed,moresweep
  integer, intent(out) :: nit,ierr
- real, intent(out)    :: errorE,errorU,EU0(:,:)
+ real, intent(out)    :: errorE,errorU,EU0(2,npart)
  integer, allocatable :: ivar(:,:),ijvar(:)
- integer              :: ncompact,ncompactlocal,npart,icompactmax,nneigh_average,its
+ integer              :: ncompact,ncompactlocal,icompactmax,nneigh_average,its
  real, allocatable    :: vari(:,:),varij(:,:),varij2(:,:),varinew(:,:)
  real :: maxerrE2,maxerrU2,maxerrE2last,maxerrU2last
  logical :: converged
@@ -157,7 +158,6 @@ subroutine do_radiation_onestep(dt,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit,
  errorU = 0.
  ierr = 0
 
- npart = size(xyzh(1,:))
  nneigh_average = int(4./3.*pi*(radkern*hfact)**3) + 1
  icompactmax = int(1.2*nneigh_average*npart)
  allocate(ivar(3,npart),stat=ierr)
@@ -175,6 +175,8 @@ subroutine do_radiation_onestep(dt,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit,
     ierr = ierr_neighbourlist_empty
     return
  endif
+
+ !$omp parallel
  call fill_arrays(ncompact,ncompactlocal,npart,icompactmax,dt,&
                   xyzh,vxyzu,ivar,ijvar,radprop,rad,vari,varij,varij2,EU0)
 
@@ -195,13 +197,16 @@ subroutine do_radiation_onestep(dt,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit,
     maxerrU2last = maxerrU2
  enddo iterations
 
+ !$omp single
  if (converged) then
     if (iverbose >= 0) print "(1x,a,i4,a,es10.3,a,es10.3)", &
           trim(label)//': succeeded with ',its,' iterations: xi err:',maxerrE2,' u err:',maxerrU2
  else
     call warning('radiation_implicit','maximum iterations reached')
     moresweep = .true.
  endif
+ !$omp end single
+ !$omp end parallel
 
  call store_radiation_results(ncompactlocal,npart,ivar,EU0,rad,vxyzu)
 
@@ -366,18 +371,17 @@ subroutine fill_arrays(ncompact,ncompactlocal,npart,icompactmax,dt,xyzh,vxyzu,iv
          dust_kappai,dust_cooling,heatingISRi,dust_gas
 
  c_code = get_c_code()
- dti = dt
- !$omp parallel do default(none) &
- !$omp shared(EU0,radprop,rad,xyzh,vxyzu,c_code,vari,ivar,ijvar,varij,varij2,dvdx,dxbound,dybound,dzbound) &
- !$omp shared(dust_temp,ncompactlocal,ncompact,massoftype,iopacity_type,nucleation,dt,gradh,cv_type) &
- !$omp firstprivate(dti) &
- !$omp private(n,i,j,k,rhoi,icompact,pmi,dvxdxi,dvxdyi,dvxdzi,dvydxi,dvydyi,dvydzi) &
+ !$omp do &
+ !!$omp shared(EU0,radprop,rad,xyzh,vxyzu,c_code,vari,ivar,ijvar,varij,varij2,dvdx,dxbound,dybound,dzbound) &
+ !!$omp shared(dust_temp,ncompactlocal,ncompact,massoftype,iopacity_type,nucleation,dt,gradh,cv_type) &
+ !$omp private(n,i,j,k,rhoi,icompact,pmi,dvxdxi,dvxdyi,dvxdzi,dvydxi,dvydyi,dvydzi,dti) &
  !$omp private(dvzdxi,dvzdyi,dvzdzi,dx,dy,dz,rij2,rij,rij1,dr,pmj,rhoj,hi,hj,hi21,hj21,hi41,hj41) &
  !$omp private(v2i,vi,v2j,vj,dWi,dWj,dvx,dvy,dvz,rhomean,dvdotdr,dv,vmu,dvdWimj,dvdWimi,dvdWjmj) &
  !$omp private(dWidrlightrhorhom,pmjdWrijrhoi,dWjdrlightrhorhom,dWiidrlightrhorhom,cv_effective) &
  !$omp private(pmjdWrunix,pmjdWruniy,pmjdWruniz,dust_kappai,dust_cooling,heatingISRi,dust_gas)
 
  do n = 1,ncompact
+    dti = dt
     i = ivar(3,n)
     !  if (iphase(i) == 0) then
     EU0(1,i) = rad(iradxi,i)
@@ -520,7 +524,7 @@ subroutine fill_arrays(ncompact,ncompactlocal,npart,icompactmax,dt,xyzh,vxyzu,iv
     vari(2,n) = rhoi
     ! endif
  enddo
- !$omp end parallel do
+ !$omp end do
 
 end subroutine fill_arrays
 
@@ -538,8 +542,8 @@ subroutine compute_flux(ivar,ijvar,ncompact,npart,icompactmax,varij,varij2,vari,
  integer             :: i,j,k,n,icompact
  real                :: rhoi,rhoj,pmjdWrunix,pmjdWruniy,pmjdWruniz,dedxi,dedyi,dedzi,dradenij
 
- !$omp parallel do default(none)&
- !$omp shared(vari,ivar,EU0,varij2,ijvar,varij,ncompact,radprop,varinew)&
+ !$omp do schedule(runtime)&
+ !!$omp shared(vari,ivar,EU0,varij2,ijvar,varij,ncompact,radprop,varinew)&
  !$omp private(i,j,k,n,dedxi,dedyi,dedzi,rhoi,rhoj,icompact)&
  !$omp private(pmjdWrunix,pmjdWruniy,pmjdWruniz,dradenij)
 
@@ -574,7 +578,7 @@ subroutine compute_flux(ivar,ijvar,ncompact,npart,icompactmax,varij,varij2,vari,
     radprop(ifluxy,i) = dedyi
     radprop(ifluxz,i) = dedzi
  enddo
- !$omp end parallel do
+ !$omp end do
 
 end subroutine compute_flux
 
@@ -596,9 +600,9 @@ subroutine calc_lambda_and_eddington(ivar,ncompactlocal,npart,vari,EU0,radprop,i
  real                :: rhoi,gradE1i,opacity,radRi
 
  ierr = 0
- !$omp parallel do default(none)&
- !$omp shared(vari,ivar,radprop,ncompactlocal,EU0,dust_temp,nucleation,limit_radiation_flux) &
- !$omp private(i,n,rhoi,gradE1i,opacity,radRi)&
+ !$omp do &
+ !!$omp shared(vari,ivar,radprop,ncompactlocal,EU0,dust_temp,nucleation,limit_radiation_flux) &
+ !$omp private(i,n,rhoi,gradE1i,opacity,radRi) &
  !$omp reduction(max:ierr)
  do n = 1,ncompactlocal
     i = ivar(3,n)
@@ -625,7 +629,7 @@ subroutine calc_lambda_and_eddington(ivar,ncompactlocal,npart,vari,EU0,radprop,i
     radprop(ilambda,i) = (2. + radRi ) / (6. + 3.*radRi + radRi**2)  ! Levermore & Pomraning's flux limiter (e.g. eq 12, Whitehouse & Bate 2004)
     radprop(iedd,i) = radprop(ilambda,i) + radprop(ilambda,i)**2 * radRi**2  ! e.g., eq 11, Whitehouse & Bate (2004)
  enddo
- !$omp end parallel do
+ !$omp end do
 
 end subroutine calc_lambda_and_eddington
 
@@ -649,8 +653,8 @@ subroutine calc_diffusion_term(ivar,ijvar,varij,ncompact,npart,icompactmax, &
  real                 :: diffusion_numerator,diffusion_denominator,tempval1,tempval2
 
  ierr = 0
- !$omp parallel do default(none)&
- !$omp shared(vari,varij,ivar,ijvar,radprop,ncompact,EU0,varinew,dust_temp,nucleation)&
+ !$omp do &
+ !!$omp shared(vari,varij,ivar,ijvar,radprop,ncompact,EU0,varinew,dust_temp,nucleation)&
  !$omp private(i,j,k,n,rhoi,rhoj,opacityi,opacityj,bi,bj,b1,diffusion_numerator,diffusion_denominator)&
  !$omp private(dWiidrlightrhorhom,dWidrlightrhorhom,dWjdrlightrhorhom,tempval1,tempval2,icompact)&
  !$omp reduction(max:ierr)
@@ -727,7 +731,7 @@ subroutine calc_diffusion_term(ivar,ijvar,varij,ncompact,npart,icompactmax, &
     !$omp atomic
     varinew(2,i) = varinew(2,i) + diffusion_denominator
  enddo
- !$omp end parallel do
+ !$omp end do
 
 end subroutine calc_diffusion_term
 
@@ -765,10 +769,10 @@ subroutine update_gas_radiation_energy(ivar,ijvar,vari,ncompact,npart,ncompactlo
  maxerrE2 = 0.
  maxerrU2 = 0.
 
- !$omp parallel do default(none)&
- !$omp shared(vari,ivar,ijvar,radprop,rad,ncompact,ncompactlocal,EU0,varinew) &
- !$omp shared(dvdx,origEU,nucleation,dust_temp,eos_vars,implicit_radiation_store_drad,cv_type) &
- !$omp shared(moresweep,pdvvisc,metallicity,vxyzu,iopacity_type,a_code,c_code,massoftype,drad,fxyzu) &
+ !$omp do &
+ !!$omp shared(vari,ivar,ijvar,radprop,rad,ncompact,ncompactlocal,EU0,varinew) &
+ !!$omp shared(dvdx,origEU,nucleation,dust_temp,eos_vars,implicit_radiation_store_drad,cv_type) &
+ !!$omp shared(moresweep,pdvvisc,metallicity,vxyzu,iopacity_type,a_code,c_code,massoftype,drad,fxyzu) &
  !$omp private(i,j,n,rhoi,dti,diffusion_numerator,diffusion_denominator,U1i,skip_quartic,Tgas,E1i,dUcomb,dEcomb) &
  !$omp private(gradEi2,gradvPi,rpdiag,rpall,radpresdenom,stellarradiation,dust_tempi,dust_kappai,xnH2) &
  !$omp private(dust_cooling,heatingISRi,dust_gas,gas_dust_val,dustgammaval,gas_dust_cooling,cosmic_ray) &
@@ -1018,7 +1022,7 @@ subroutine update_gas_radiation_energy(ivar,ijvar,vari,ncompact,npart,ncompactlo
     endif
 
  enddo main_loop
-!$omp end parallel do
+!$omp end do
 
 end subroutine update_gas_radiation_energy