Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- This code takes a list of positions x,y,z and creates a kernal smoothed mesh using smoothing lengths hsml.
- Basically, I want to just to loop through i as fast a possible as every calculation is independent of one another. It just has to populate mesh_out which is shared. kern_result is a function representing the smoothing kernel.
- !$OMP PARALLEL DO DEFAULT(private) SHARED(mesh_out)
- do i = 1,Ntot
- xtmp = x0(i)
- ytmp = y0(i)
- ztmp = z0(i)
- hsmltmp = hsml(i)
- if (h.lt.hmin) then
- hsmltmp = hmin
- elseif (h.gt.hmax) then
- hsmltmp = hmax
- end if
- Xrel = xtmp-Xc+Lx/2;
- Yrel = ytmp-Yc+Ly/2;
- Zrel = ztmp-Zc+Lz/2;
- if (Xrel.gt.-hsmltmp.and.Yrel.gt.-hsmltmp.and.Xrel-Lx.lt.hsmltmp &
- .and.Yrel-Ly.lt.hsmltmp.and.Zrel.gt.-hsmltmp.and.Zrel-Lz.lt.hsmltmp) then
- binned_particles = binned_particles + 1
- pixelX = (FLOOR(Xrel/pixelsizeX)+.5)*pixelsizeX
- pixelY = (FLOOR(Yrel/pixelsizeY)+.5)*pixelsizeY
- pixelZ = (FLOOR(Zrel/pixelsizeZ)+.5)*pixelsizeZ
- !write(*,*)"pixelXYZ",pixelX,pixelY,pixelZ
- nx = INT(CEILING(hsmltmp/pixelsizeX))
- ny = INT(CEILING(hsmltmp/pixelsizeY))
- nz = INT(CEILING(hsmltmp/pixelsizeZ))
- !write(*,*)"Nxyz",nx,ny,nz
- normalization = 0
- do dx = -nx,nx
- do dy = -ny,ny
- do dz = -nz,nz
- distX = Xrel-(pixelX-REAL(dx)*pixelsizeX)
- distY = Yrel-(pixelY-REAL(dy)*pixelsizeY)
- distZ = Zrel-(pixelZ-REAL(dz)*pixelsizeZ)
- r2 = distX**2+distY**2+distZ**2;
- if (r2.lt.hsmltmp**2) then
- kern_result = kernel(hsmltmp,r2)
- normalization = normalization + kern_result
- xloc = CEILING(pixelX/pixelsizeX-REAL(dx));
- yloc = CEILING(pixelY/pixelsizeY-REAL(dy));
- zloc = CEILING(pixelZ/pixelsizeZ-REAL(dz));
- !write(*,*)"XYZloc",xloc,yloc,zloc
- if (xloc.gt.0.and.xloc.le.num_pixelsX.and.yloc.gt.0.and.yloc.le.num_pixelsY.and. &
- zloc.gt.0.and.zloc.le.num_pixelsZ) then
- kern_result = kernel(hsmltmp,r2)
- !$OMP ATOMIC UPDATE
- mesh_out(INT(zloc),INT(yloc),INT(xloc)) = &
- mesh_out(INT(zloc),INT(yloc),INT(xloc)) &
- + kern_result/normalization
- !write(*,*)"kern",hsmltmp,r2,kern_result
- end if
- end if
- end do
- end do
- end do
- end if
- end do
- !$OMP END PARALLEL DO
- RECURSIVE FUNCTION kernel(hdum,rdum) RESULT(kern)
- real*8 rdum
- real*8 hdum
- real*8 kern
- ! Kernal Coefficients
- real*8, parameter :: K1 = 2.546479089470
- real*8, parameter :: K2 = 15.278874536822
- real*8, parameter :: K3 = 5.092958178941
- real*8 u
- u = sqrt(rdum)/hdum;
- if (u.lt.0.5) then
- kern = K1+K2*(u-1)*u**2;
- else
- kern = K3*(1-u)**3;
- end if
- END FUNCTION kernel
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement