From 117f02008b9c54e0e14c3edaa36f5836be641aae Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 29 Jun 2023 17:33:55 +1000 Subject: [PATCH 1/3] (radiation) merge parallel sections --- src/main/radiation_implicit.f90 | 43 ++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/src/main/radiation_implicit.f90 b/src/main/radiation_implicit.f90 index 26e7d1494..22ed6229c 100644 --- a/src/main/radiation_implicit.f90 +++ b/src/main/radiation_implicit.f90 @@ -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,6 +197,7 @@ 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 @@ -202,6 +205,8 @@ subroutine do_radiation_onestep(dt,rad,xyzh,vxyzu,radprop,origEU,EU0,failed,nit, 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) @@ -367,9 +372,9 @@ subroutine fill_arrays(ncompact,ncompactlocal,npart,icompactmax,dt,xyzh,vxyzu,iv 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 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 firstprivate(dti) & !$omp private(n,i,j,k,rhoi,icompact,pmi,dvxdxi,dvxdyi,dvxdzi,dvydxi,dvydyi,dvydzi) & !$omp private(dvzdxi,dvzdyi,dvzdzi,dx,dy,dz,rij2,rij,rij1,dr,pmj,rhoj,hi,hj,hi21,hj21,hi41,hj41) & @@ -520,7 +525,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 +543,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 +579,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 +601,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 +630,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 +654,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 +732,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 +770,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 +1023,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 From 9e77708182f2088b3b2023f3688be7c6aab7bcf4 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Fri, 30 Jun 2023 15:12:56 +1000 Subject: [PATCH 2/3] (radiation) change dti to private --- src/main/radiation_implicit.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/main/radiation_implicit.f90 b/src/main/radiation_implicit.f90 index 22ed6229c..b4ec18ae7 100644 --- a/src/main/radiation_implicit.f90 +++ b/src/main/radiation_implicit.f90 @@ -371,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 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 firstprivate(dti) & - !$omp private(n,i,j,k,rhoi,icompact,pmi,dvxdxi,dvxdyi,dvxdzi,dvydxi,dvydyi,dvydzi) & + !$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) From 71551968227128285477a26c6410bb3a070e977f Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 6 Jul 2023 15:45:18 +1000 Subject: [PATCH 3/3] (radiation-implicit) fix compiler warning --- src/main/radiation_implicit.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/main/radiation_implicit.f90 b/src/main/radiation_implicit.f90 index b4ec18ae7..aa7596d81 100644 --- 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)