diff --git a/doxygen/pafiX_documentation b/documentation/doxygen_codeGuide/codeGuide
similarity index 99%
rename from doxygen/pafiX_documentation
rename to documentation/doxygen_codeGuide/codeGuide
index 1b569f4af512a20b2786d3e9275d8d188862571a..1eb64fbd1fa493d85db465b598df4223e4c1bf8a 100644
--- a/doxygen/pafiX_documentation
+++ b/documentation/doxygen_codeGuide/codeGuide
@@ -58,7 +58,7 @@ PROJECT_LOGO = pafiXlogo.png
# entered, it will be relative to the location where doxygen was started. If
# left blank the current directory will be used.
-OUTPUT_DIRECTORY = ../../../../projects/pafiX/documentation
+OUTPUT_DIRECTORY = .
# If the CREATE_SUBDIRS tag is set to YES then doxygen will create 4096 sub-
# directories (in 2 levels) under the output directory of each output format and
@@ -771,7 +771,7 @@ WARN_LOGFILE =
# spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
# Note: If this tag is empty the current directory is searched.
-INPUT = ../src
+INPUT = ../../src
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/doxygen/pafiXlogo.png b/documentation/doxygen_codeGuide/pafiXlogo.png
similarity index 100%
rename from doxygen/pafiXlogo.png
rename to documentation/doxygen_codeGuide/pafiXlogo.png
diff --git a/documentation/userGuide/pafiX_userGuide.pdf b/documentation/userGuide/pafiX_userGuide.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..c0791cc0c383d8fa177d5aad7770bf0d480bce36
Binary files /dev/null and b/documentation/userGuide/pafiX_userGuide.pdf differ
diff --git a/input/input.dat b/input/input_example.dat
similarity index 55%
rename from input/input.dat
rename to input/input_example.dat
index 44b08ceedae1a5c818b2288c35e420bad02c527e..0aa6a278b383b55b784d8c0bd665ca4a06e1d491 100644
--- a/input/input.dat
+++ b/input/input_example.dat
@@ -1,24 +1,28 @@
pafiX input data
----------------
dimensions in x,y,z-direction (m,m,m) =
-0.24,0.04,0.04
+0.24,0.125,0.04
+boundary conditions in x,y,z-direction [(w)all/(p)eriodic/(i)n-outlet] =
+p,p,w
number of cells in x,y,z-direction (-,-,-) =
-60,60,60
+20,20,20
+grid in x,y,z-direction [(u)niform/(s)tretch] =
+u,u,s
start time-step (-), time-steps to compute (-) =
-1,10000
+1,10
interval of writing result files, restart files (-,-) =
10,10000
CFL number (-) =
-0.3
-initial bulk flow (m/s), pressure gradient (N/m**3) =
-3.65,-5.8
+0.4
+in/initial flow (m/s), pressure gradient (N/m**3) =
+3.65,-2.9
fluid density (kg/m**3), fluid kinematic viscosity (m**2/s) =
1.2,1.46E-5
-particle number density (-/m**3), seet at nt (-) =
-1E8,10000
+particle number density (-/m**3) or rate (-/s), start seed at nt (-) =
+1E7,1
particle radius (m), material density (kg/m**3) =
50E-6,1200
-particle charge (C) =
-1.E-13
+particle charge: initial, equilibrium (C,C) =
+0.,1E-12
gravity vector in x,y,z-dir. (m/s**2,m/s**2,m/s**2) =
9.81,0,0
diff --git a/monitor.plt b/monitor.plt
index 95e01264a1a5518fe9a61779863f146eb3fd9128..063b67b21ba26e804611c24678beeb617a29a425 100644
--- a/monitor.plt
+++ b/monitor.plt
@@ -1,37 +1,17 @@
-#set term postscript enhanced eps 10 color linewidth 1 'Times-Roman'
-#set term dumb
-#set output "post.eps"
-#set xrange [0:]
-set multiplot layout 3,2 title 'pafiX -- by H. Grosshans'
+set multiplot layout 4,4 title 'pafiX -- by H. Grosshans'
f = 'output/monitor.dat'
+header=system("head -3 ".f." | tail -1")
set key box opaque width 2 t l
set xlabel offset 0,0.5
set ylabel offset 2
set linetype 1 ps 0.5 lc 3
-set linetype 2 ps 0.5 lc 4
-set xlabel 's'
+set xlabel word(header,3) #header columns are shifted by 1 bc of #
-set ylabel 'm/s'
-p f u 2:4 t 'u_{cl}', f u 2:5 t 'u_{av}'
-
-set ylabel '-'
-p f u 2:8 t 'Re_c', f u 2:9 t 'Re_b'
-
-set ylabel 'm/s'
-p f u 2:6 t 'v_{rms}', f u 2:7 t 'w_{rms}'
-
-set ylabel '-'
-p f u 2:($6/$11) t 'v_{rms}/u_{tau}', f u 2:($7/$11) t 'w_{rms}/u_{tau}'
-
-#set ylabel '1/s'
-#p f u 2:14 t 'gradu'
-
-set ylabel 'N/m**2'
-p f u 2:10 t 'tau_w'
-
-set ylabel '-'
-p f u 2:12 t 'Re_{tau}'
+do for [i=4:16] {
+set ylabel word(header,i+1)
+p f u 2:i notitle ls 1
+}
pause 20
reread
diff --git a/run.sh b/run.sh
index a49e0e732d54eb3fe2615bce5f38596d9c4c1e33..aa8cf587cbcd12f195bd78dfd4e87712d4088797 100644
--- a/run.sh
+++ b/run.sh
@@ -1 +1,2 @@
-mpiexec -np 4 src/pafiX
+mkdir -p results restart output
+mpiexec -np 2 --oversubscribe src/pafiX
diff --git a/src/bc.f90 b/src/bc.f90
index d9c22d6d29d0b21225d4ee4b89119d1f41717c0b..a211c927a04087700e29ff5c7b9fc3c29629aa9a 100644
--- a/src/bc.f90
+++ b/src/bc.f90
@@ -1,42 +1,128 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief boundary condition for velocity field
- subroutine bc
+!> @brief boundary condition for velocity and pressure field
+ subroutine bcUVWP
use var
real(kind=pr) :: &
ue,uw,ve,vw,we,ww, &
uba,ufr,vba,vfr,wba,wfr, &
- un,us,vn,vs,wn,ws
- integer :: m
+ un,us,vn,vs,wn,ws, &
+ alpha,dist
+ integer :: i,j,l,m,mm
-! ue=0._pr; uw=0._pr; ve=0._pr; vw=0._pr; ww=0._pr; we=0._pr
+ if (bcx.eq.'w') then
+ ue=0._pr; uw=0._pr; ve=0._pr; vw=0._pr; ww=0._pr; we=0._pr
+ endif
+ if (bcy.eq.'w') then
uba=0._pr; ufr=0._pr; vba=0._pr; vfr=0._pr; wba=0._pr; wfr=0._pr
+ endif
+ if (bcz.eq.'w') then
un =0._pr; us =0._pr; vs =0._pr; vn =0._pr; wn =0._pr; ws =0._pr
+ endif
+
+
+ if (bcx.eq.'w') then
+ if (myid.eq.0) then
+ do 1 m=1,gc
+ u(imin-m,:,:)= uw
+ v(imin-m,:,:)= v(imin,:,:)-2._pr*(v(imin,:,:)-vw)
+ w(imin-m,:,:)= w(imin,:,:)-2._pr*(w(imin,:,:)-ww)
+1 enddo
+ endif
+ if (myid.eq.nrprocs-1) then
+ do 2 m=1,gc
+ u(imax-1+m,:,:)=ue
+ v(imax+m,:,:)= v(imax,:,:)-2._pr*(v(imax,:,:)-ve)
+ w(imax+m,:,:)= w(imax,:,:)-2._pr*(w(imax,:,:)-we)
+2 enddo
+ endif
+ elseif (bcx.eq.'p') then
+ ! done through the sync routines
+ elseif (bcx.eq.'i') then
+ ! see inflow routine
+ if (myid.eq.nrprocs-1) then
+ do 8 m=1,gc
+ u(imax+m-1,:,:)=u(imax-1,:,:)
+ v(imax+m,:,:)= v(imax,:,:)
+ w(imax+m,:,:)= w(imax,:,:)
+ p(imax+m,:,:)= p(imax,:,:)
+8 enddo
+ endif
+ endif
+
+ if (bcy.eq.'w') then
+ do 3 m=1,gc
+ u(:,jmin-m,:)= u(:,jmin,:)-2._pr*(u(:,jmin,:)-ufr)
+ u(:,jmax+m,:)= u(:,jmax,:)-2._pr*(u(:,jmax,:)-uba)
+ v(:,jmin-m,:)= vfr
+ v(:,jmax-1+m,:)=vba
+ w(:,jmin-m,:)= w(:,jmin,:)-2._pr*(w(:,jmin,:)-wfr)
+ w(:,jmax+m,:)= w(:,jmax,:)-2._pr*(w(:,jmax,:)-wba)
+3 enddo
+ elseif (bcy.eq.'p') then
+ do 4 m=1,gc
+ u(:,jmin-m,:)= u(:,jmax+1-m,:)
+ u(:,jmax+m,:)= u(:,jmin-1+m,:)
+ v(:,jmin-m,:)= v(:,jmax+1-m,:)
+ v(:,jmax+m,:)= v(:,jmin-1+m,:)
+ w(:,jmin-m,:)= w(:,jmax+1-m,:)
+ w(:,jmax+m,:)= w(:,jmin-1+m,:)
+ p(:,jmin-m,:)= p(:,jmax+1-m,:)
+ p(:,jmax+m,:)= p(:,jmin-1+m,:)
+4 enddo
+ endif
+
+ if (bcz.eq.'w') then
+ do 5 m=1,gc
+ u(:,:,lmin-m)= u(:,:,lmin)-2._pr*(u(:,:,lmin)-us)
+ u(:,:,lmax+m)= u(:,:,lmax)-2._pr*(u(:,:,lmax)-un)
+ v(:,:,lmin-m)= v(:,:,lmin)-2._pr*(v(:,:,lmin)-vs)
+ v(:,:,lmax+m)= v(:,:,lmax)-2._pr*(v(:,:,lmax)-vn)
+ w(:,:,lmin-m)= ws
+ w(:,:,lmax-1+m)=wn
+5 enddo
+ elseif (bcz.eq.'p') then
+ do 6 m=1,gc
+ u(:,:,lmin-m)= u(:,:,lmax+1-m)
+ u(:,:,lmax+m)= u(:,:,lmin-1+m)
+ v(:,:,lmin-m)= v(:,:,lmax+1-m)
+ v(:,:,lmax+m)= v(:,:,lmin-1+m)
+ w(:,:,lmin-m)= w(:,:,lmax+1-m)
+ w(:,:,lmax+m)= w(:,:,lmin-1+m)
+ p(:,:,lmin-m)= p(:,:,lmax+1-m)
+ p(:,:,lmax+m)= p(:,:,lmin-1+m)
+6 enddo
+ endif
+
+
+ return
+ end
+
+!####################################################################
+!> @author Holger Grosshans
+!> @brief inflow boundary condition
+ subroutine inflow
+ use var
+ real(kind=pr),dimension(dimj*diml) :: random
+ real(kind=pr) :: &
+ alpha,dist
+ integer :: i,j,l,m,mm
-! commented because of periodic bc in x
-! u(imin-1,j,l)=uw
-! u(imax,j,l) =ue
-! v(imin-1,j,l)=v(imin)-2._pr*(v(imin,j,l)-vw)
-! v(imax+1,j,l)=v(imax)-2._pr*(v(imax,j,l)-ve)
-! w(imin-1,j,l)=w(imin)-2._pr*(w(imin,j,l)-ww)
-! w(imax+1,j,l)=w(imax)-2._pr*(w(imax,j,l)-we)
-
- do 1 m=1,wc
- u(:,jmin-m,:) =u(:,jmin,:)-2._pr*(u(:,jmin,:)-ufr)
- u(:,jmax+m,:) =u(:,jmax,:)-2._pr*(u(:,jmax,:)-uba)
- v(:,jmin-m,:) =vfr
- v(:,jmax+1-m,:)=vba
- w(:,jmin-m,:) =w(:,jmin,:)-2._pr*(w(:,jmin,:)-wfr)
- w(:,jmax+m,:) =w(:,jmax,:)-2._pr*(w(:,jmax,:)-wba)
-
- u(:,:,lmin-m) =u(:,:,lmin)-2._pr*(u(:,:,lmin)-us)
- u(:,:,lmax+m) =u(:,:,lmax)-2._pr*(u(:,:,lmax)-un)
- v(:,:,lmin-m) =v(:,:,lmin)-2._pr*(v(:,:,lmin)-vs)
- v(:,:,lmax+m) =v(:,:,lmax)-2._pr*(v(:,:,lmax)-vn)
- w(:,:,lmin-m) =ws
- w(:,:,lmax+1-m)=wn
-1 enddo
+ call random_number(random)
+ alpha=.1_pr
+ mm=0
+ do j=jmin,jmax; do l=lmin,lmax
+ mm=mm+1
+ if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l)))
+ if (bcy.eq.'w'.and.bcz.eq.'p') dist=delta-abs(yc(j))
+ if (bcy.eq.'p'.and.bcz.eq.'w') dist=delta-abs(zc(l))
+ do m=1,gc
+ u(imin-m,j,l)=max(0._pr,ubulk0*dist/delta + (random(mm)-.5_pr)*alpha*ubulk0)
+ v(imin-m,j,l)=0._pr
+ w(imin-m,j,l)=0._pr
+ enddo
+ enddo; enddo
return
end
@@ -44,75 +130,216 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief boundary condition for electric potential
- subroutine bcPhi_el
+!> @brief boundary condition for electric field
+ subroutine bcE_el
use var
integer :: m
- do 1 m=1,wc
- phi_el(:,jmin-m,:)= 0._pr !zero on a conductive grounded surface
- phi_el(:,jmax+m,:)= 0._pr
-! phi_el(:,jmin-m,:)= phi_el(:,jmin,:) !zero gradient on a conductive grounded surface
-! phi_el(:,jmax+m,:)= phi_el(:,jmax,:)
+ if (bcx.eq.'w') then
+ if (myid.eq.0) then
+ do 1 m=1,gc
+ Ex_el(imin-m,:,:)= Ex_el(imin,:,:)
+ Ey_el(imin-m,:,:)= Ey_el(imin,:,:)
+ Ez_el(imin-m,:,:)= Ez_el(imin,:,:)
+1 enddo
+ endif
+ if (myid.eq.nrprocs-1) then
+ do 2 m=1,gc
+ Ex_el(imax+m,:,:)= Ex_el(imax,:,:)
+ Ey_el(imax+m,:,:)= Ey_el(imax,:,:)
+ Ez_el(imax+m,:,:)= Ez_el(imax,:,:)
+2 enddo
+ endif
+ elseif (bcx.eq.'p') then
+ ! done through the sync routines
+ elseif (bcx.eq.'i') then
+ if (myid.eq.0) then
+ do 7 m=1,gc
+ Ex_el(imin-m,:,:)= Ex_el(imin,:,:)
+ Ey_el(imin-m,:,:)= Ey_el(imin,:,:)
+ Ez_el(imin-m,:,:)= Ez_el(imin,:,:)
+7 enddo
+ endif
+ if (myid.eq.nrprocs-1) then
+ do 8 m=1,gc
+ Ex_el(imax+m,:,:)= Ex_el(imax,:,:)
+ Ey_el(imax+m,:,:)= Ey_el(imax,:,:)
+ Ez_el(imax+m,:,:)= Ez_el(imax,:,:)
+8 enddo
+ endif
+ endif
- phi_el(:,:,lmin-m)= 0._pr
- phi_el(:,:,lmax+m)= 0._pr
-! phi_el(:,:,lmin-m)= phi_el(:,:,lmin) !zero gradient on a conductive grounded surface
-! phi_el(:,:,lmax+m)= phi_el(:,:,lmax)
-1 enddo
+ if (bcy.eq.'w') then
+ do 3 m=1,gc
+ Ex_el(:,jmin-m,:)= Ex_el(:,jmin,:)
+ Ex_el(:,jmax+m,:)= Ex_el(:,jmax,:)
+ Ey_el(:,jmin-m,:)= Ey_el(:,jmin,:)
+ Ey_el(:,jmax+m,:)= Ey_el(:,jmax,:)
+ Ez_el(:,jmin-m,:)= Ez_el(:,jmin,:)
+ Ez_el(:,jmax+m,:)= Ez_el(:,jmax,:)
+3 enddo
+ elseif (bcy.eq.'p') then
+ do 4 m=1,gc
+ Ex_el(:,jmin-m,:)= Ex_el(:,jmax+1-m,:)
+ Ex_el(:,jmax+m,:)= Ex_el(:,jmin-1+m,:)
+ Ey_el(:,jmin-m,:)= Ey_el(:,jmax+1-m,:)
+ Ey_el(:,jmax+m,:)= Ey_el(:,jmin-1+m,:)
+ Ez_el(:,jmin-m,:)= Ez_el(:,jmax+1-m,:)
+ Ez_el(:,jmax+m,:)= Ez_el(:,jmin-1+m,:)
+4 enddo
+ endif
+ if (bcz.eq.'w') then
+ do 5 m=1,gc
+ Ex_el(:,:,lmin-m)= Ex_el(:,:,lmin)
+ Ex_el(:,:,lmax+m)= Ex_el(:,:,lmax)
+ Ey_el(:,:,lmin-m)= Ey_el(:,:,lmin)
+ Ey_el(:,:,lmax+m)= Ey_el(:,:,lmax)
+ Ez_el(:,:,lmin-m)= Ez_el(:,:,lmin)
+ Ez_el(:,:,lmax+m)= Ez_el(:,:,lmax)
+5 enddo
+ elseif (bcz.eq.'p') then
+ do 6 m=1,gc
+ Ex_el(:,:,lmin-m)= Ex_el(:,:,lmax+1-m)
+ Ex_el(:,:,lmax+m)= Ex_el(:,:,lmin-1+m)
+ Ey_el(:,:,lmin-m)= Ey_el(:,:,lmax+1-m)
+ Ey_el(:,:,lmax+m)= Ey_el(:,:,lmin-1+m)
+ Ez_el(:,:,lmin-m)= Ez_el(:,:,lmax+1-m)
+ Ez_el(:,:,lmax+m)= Ez_el(:,:,lmin-1+m)
+6 enddo
+ endif
-! phi_el(i,j,l) = 2*phi_el(i,j+1,l)-phi_el(i,j+2,l) !constant gradient on a non-dielectric insulater, i.e. a surface which does not influence the field
-! phi_el(i,j,l) = 1.25*phi_el(i,j,l+1)-0.25*phi_el(i,j,l+2)
return
end
+
!####################################################################
!> @author Holger Grosshans
-!> @brief boundary condition for electric field
- subroutine bcE_el
+!> @brief boundary condition for electric potential
+ subroutine bcPhi_el
use var
+ real(kind=pr) :: wval
integer :: m
+ wval= 0._pr !zero on a conductive grounded surface
+
+ if (bcx.eq.'w') then
+ if (myid.eq.0) then
+ do 1 m=1,gc
+ phi_el(imin-m,:,:)= -phi_el(imin,:,:)+2._pr*wval
+1 enddo
+ endif
+ if (myid.eq.nrprocs-1) then
+ do 2 m=1,gc
+ phi_el(imax+m,:,:)= -phi_el(imax,:,:)+2._pr*wval
+2 enddo
+ endif
+ elseif (bcx.eq.'p') then
+ ! done through the sync routines
+ elseif (bcx.eq.'i') then
+ if (myid.eq.0) then
+ do 7 m=1,gc
+ phi_el(imin-m,:,:)= phi_el(imin,:,:)
+7 enddo
+ endif
+ if (myid.eq.nrprocs-1) then
+ do 8 m=1,gc
+ phi_el(imax+m,:,:)= phi_el(imax,:,:)
+8 enddo
+ endif
+ endif
- do 1 m=1,wc
- Ex_el(:,jmin-m,:)= Ex_el(:,jmin,:)
- Ex_el(:,jmax+m,:)= Ex_el(:,jmax,:)
- Ey_el(:,jmin-m,:)= Ey_el(:,jmin,:)
- Ey_el(:,jmax+m,:)= Ey_el(:,jmax,:)
- Ez_el(:,jmin-m,:)= Ez_el(:,jmin,:)
- Ez_el(:,jmax+m,:)= Ez_el(:,jmax,:)
+ if (bcy.eq.'w') then
+ do 3 m=1,gc
+ phi_el(:,jmin-m,:)= -phi_el(:,jmin,:)+2._pr*wval
+ phi_el(:,jmax+m,:)= -phi_el(:,jmax,:)+2._pr*wval
+3 enddo
+ elseif (bcy.eq.'p') then
+ do 4 m=1,gc
+ phi_el(:,jmin-m,:)= phi_el(:,jmax+1-m,:)
+ phi_el(:,jmax+m,:)= phi_el(:,jmin-1+m,:)
+4 enddo
+ endif
- Ex_el(:,:,lmin-m)= Ex_el(:,:,lmin)
- Ex_el(:,:,lmax+m)= Ex_el(:,:,lmax)
- Ey_el(:,:,lmin-m)= Ey_el(:,:,lmin)
- Ey_el(:,:,lmax+m)= Ey_el(:,:,lmax)
- Ez_el(:,:,lmin-m)= Ez_el(:,:,lmin)
- Ez_el(:,:,lmax+m)= Ez_el(:,:,lmax)
-1 enddo
+ if (bcz.eq.'w') then
+ do 5 m=1,gc
+ phi_el(:,:,lmin-m)= -phi_el(:,:,lmin)+2._pr*wval
+ phi_el(:,:,lmax+m)= -phi_el(:,:,lmax)+2._pr*wval
+5 enddo
+ elseif (bcz.eq.'p') then
+ do 6 m=1,gc
+ phi_el(:,:,lmin-m)= phi_el(:,:,lmax+1-m)
+ phi_el(:,:,lmax+m)= phi_el(:,:,lmin-1+m)
+6 enddo
+ endif
return
end
+
!####################################################################
!> @author Holger Grosshans
!> @brief boundary condition for rho_el
subroutine bcRho_el
use var
+ real(kind=pr) :: wval
integer :: m
+ wval= 0._pr !no charge on surface
- do 1 m=1,wc
- rho_el(:,jmin-m,:)= 0._pr !no charge on surface
- rho_el(:,jmax+m,:)= 0._pr
- rho_el(:,:,lmin-m)= 0._pr
- rho_el(:,:,lmax+m)= 0._pr
-1 enddo
+ if (bcx.eq.'w') then
+ if (myid.eq.0) then
+ do 1 m=1,gc
+ rho_el(imin-m,:,:)= -rho_el(imin,:,:)+2._pr*wval
+1 enddo
+ endif
+ if (myid.eq.nrprocs-1) then
+ do 2 m=1,gc
+ rho_el(imax+m,:,:)= -rho_el(imax,:,:)+2._pr*wval
+2 enddo
+ endif
+ elseif (bcx.eq.'p') then
+ ! done through the sync routines
+ elseif (bcx.eq.'i') then
+ if (myid.eq.0) then
+ do 7 m=1,gc
+ rho_el(imin-m,:,:)= rho_el(imin,:,:)
+7 enddo
+ endif
+ if (myid.eq.nrprocs-1) then
+ do 8 m=1,gc
+ rho_el(imax+m,:,:)= rho_el(imax,:,:)
+8 enddo
+ endif
+ endif
+
+ if (bcy.eq.'w') then
+ do 3 m=1,gc
+ rho_el(:,jmin-m,:)= -rho_el(:,jmin,:)+2._pr*wval
+ rho_el(:,jmax+m,:)= -rho_el(:,jmax,:)+2._pr*wval
+3 enddo
+ elseif (bcy.eq.'p') then
+ do 4 m=1,gc
+ rho_el(:,jmin-m,:)= rho_el(:,jmax+1-m,:)
+ rho_el(:,jmax+m,:)= rho_el(:,jmin-1+m,:)
+4 enddo
+ endif
+
+ if (bcz.eq.'w') then
+ do 5 m=1,gc
+ rho_el(:,:,lmin-m)= -rho_el(:,:,lmin)+2._pr*wval
+ rho_el(:,:,lmax+m)= -rho_el(:,:,lmax)+2._pr*wval
+5 enddo
+ elseif (bcz.eq.'p') then
+ do 6 m=1,gc
+ rho_el(:,:,lmin-m)= rho_el(:,:,lmax+1-m)
+ rho_el(:,:,lmax+m)= rho_el(:,:,lmin-1+m)
+6 enddo
+ endif
return
end
-
diff --git a/src/electrostatics.f90 b/src/electrostatics.f90
index aae8c7faedcffb6946631eb1ea9bce34105faa4c..142a19f5d5be6bf879758148c9ad9facb1c61a9b 100644
--- a/src/electrostatics.f90
+++ b/src/electrostatics.f90
@@ -3,11 +3,11 @@
!> @brief solve E field
subroutine solveElectrostatics
use var
- real(kind=pr) err,err00,err2
+ real(kind=pr) err,err00,err2,syncMax
integer it
- if (pnd.ne.0._pr.and.qelp.ne.0._pr) then
+ if (pnd.ne.0._pr.and.syncMax(maxval(abs(q_el))).ne.0._pr) then
do it=1,itmax
call calcPhi_el(err)
@@ -34,7 +34,7 @@
subroutine calcPhi_el(err)
use var
real(kind=pr),dimension(ii,jj,ll) :: phi_el1
- real(kind=pr) :: diag,dex,dey,dez,err,dphi_el
+ real(kind=pr) :: diag,dex,dey,dez,err,dphi_el,syncSum
integer :: i,j,l
! phi_el: electrostatic potential (V)
@@ -44,12 +44,9 @@
err=0._pr
phi_el1=0._pr
+ do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax
- do 1 i=imin,imax,1
- do 1 j=jmin,jmax,1
- do 1 l=lmin,lmax,1
-
- if (celltype(i,j,l).ne.active) cycle
+ if (celltype(i,j,l).ne.active) cycle
! solve the poisson equation:
diag= (2._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
@@ -70,13 +67,13 @@
dphi_el=phi_el1(i,j,l)-phi_el(i,j,l)
err=err+dphi_el*dphi_el
-1 enddo
+ enddo; enddo; enddo
phi_el=phi_el1
call sync(phi_el)
call bcPhi_el
- call syncSum(err); err=(err/dimgptot)**(0.5_pr)
+ err=(syncSum(err)/dimgptot)**(0.5_pr)
return
end
@@ -87,7 +84,7 @@
!> @param err L2 error norm
subroutine calcE_el(err)
use var
- real(kind=pr) :: dexel,deyel,dezel,temp,err
+ real(kind=pr) :: dexel,deyel,dezel,temp,err,syncSum
integer :: i,j,l
! phi_el: electrostatic potential (V)
@@ -97,11 +94,9 @@
err=0._pr
- do 1 i=imin,imax,1
- do 1 j=jmin,jmax,1
- do 1 l=lmin,lmax,1
+ do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax
- if (celltype(i,j,l).ne.active) cycle
+ if (celltype(i,j,l).ne.active) cycle
! calculate electrostatic field
temp= Ex_el(i,j,l)
@@ -118,11 +113,11 @@
err=err+dexel*dexel+deyel*deyel+dezel*dezel
-1 enddo
+ enddo; enddo; enddo
call sync(Ex_el); call sync(Ey_el); call sync(Ez_el)
call bcE_el
- call syncSum(err); err=(err/dimgptot/3._pr)**(0.5_pr)
+ err=(syncSum(err)/dimgptot/3._pr)**(0.5_pr)
! write(*,'(x,a,es9.2e2,a,3(es9.2e2))') &
! 'res. E =',err
@@ -149,14 +144,12 @@
call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0)
- do 5 i=ibeg,iend
- do 5 j=jbeg,jend
- do 5 l=lbeg,lend
+ do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend
iw=i+1-ibeg
jw=j+1-jbeg
lw=l+1-lbeg
rho_el(i,j,l) = rho_el(i,j,l) + q_el(n)*partn(n)*weight(iw,jw,lw)/volE
-5 enddo
+ enddo; enddo; enddo
1 enddo
@@ -186,16 +179,14 @@
call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0)
- do 5 i=ibeg,iend
- do 5 j=jbeg,jend
- do 5 l=lbeg,lend
+ do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend
iw=i+1-ibeg
jw=j+1-jbeg
lw=l+1-lbeg
Extot=Extot+Ex_el(i,j,l)*weight(iw,jw,lw)
Eytot=Eytot+Ey_el(i,j,l)*weight(iw,jw,lw)
Eztot=Eztot+Ez_el(i,j,l)*weight(iw,jw,lw)
-5 enddo
+ enddo; enddo; enddo
partmass=4._pr/3._pr*pi*rhop*partn(n)*radp(n)**3
fx_el(n)= q_el(n)*Extot/partmass
@@ -215,51 +206,72 @@
use var
real(kind=pr) :: dis,dis3,disx,disy,disz, &
Extot,Eytot,Eztot,partmass
- integer :: n1,n2,ip,jp,lp
+ integer :: n1,n2,ip,jp,lp,m1,m2,i,j,l
- if (elforceScheme.eq.2) then
+ select case (elforceScheme)
+ case (2)
fx_el=0._pr
fy_el=0._pr
fz_el=0._pr
- endif
- do 1 n1=2,np
- do 2 n2=1,(n1-1)
-
- if (n1.eq.n2) cycle
- if ((elforceScheme.eq.3).and. &
- ((ip(n1).lt.ip(n2)-1.or.ip(n1).gt.ip(n2)+1).or. &
- (jp(n1).lt.jp(n2)-1.or.jp(n1).gt.jp(n2)+1).or. &
- (lp(n1).lt.lp(n2)-1.or.lp(n1).gt.lp(n2)+1))) cycle
-
-! force of n2 acting on n1
- disx=xp(n1)-xp(n2)
- disy=yp(n1)-yp(n2)
- disz=zp(n1)-zp(n2)
- dis=(disx**2+disy**2+disz**2)**0.5
- dis3=dis**3
-
- partmass=4._pr/3._pr*pi*rhop*partn(n1)*radp(n1)**3
+ do 1 n1=2,np
+ do 2 n2=1,(n1-1)
+ if (n1.eq.n2) cycle
+ call forcesCoulombN1N2(n1,n2)
+2 enddo
+1 enddo
+
+ case (3)
+ do i=1,ii; do j=1,jj; do l=1,ll
+
+ do 4 m1=2,npic(i,j,l)
+ do 5 m2=1,(m1-1)
+ n1=nic(i,j,l,m1)
+ n2=nic(i,j,l,m2)
+ call forcesCoulombN1N2(n1,n2)
+5 enddo
+4 enddo
+ enddo; enddo; enddo
- fx_el(n1)= fx_el(n1) &
- + q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
- fy_el(n1)= fy_el(n1) &
- + q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
- fz_el(n1)= fz_el(n1) &
- + q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
+ end select
- partmass=4._pr/3._pr*pi*rhop*partn(n2)*radp(n2)**3
+ return
+ end
- fx_el(n2)= fx_el(n2) &
- - q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
- fy_el(n2)= fy_el(n2) &
- - q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
- fz_el(n2)= fz_el(n2) &
- - q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
-2 enddo
-1 enddo
+!####################################################################
+!> @author Holger Grosshans
+!> @brief compute Coulomb force of n2 acting on n1
+ subroutine forcesCoulombN1N2(n1,n2)
+ use var
+ real(kind=pr) :: dis,dis3,disx,disy,disz,partmass
+ integer :: n1,n2,ip,jp,lp
+
+
+ disx=xp(n1)-xp(n2)
+ disy=yp(n1)-yp(n2)
+ disz=zp(n1)-zp(n2)
+ dis=(disx**2+disy**2+disz**2)**0.5
+ dis3=dis**3
+
+ partmass=4._pr/3._pr*pi*rhop*partn(n1)*radp(n1)**3
+
+ fx_el(n1)= fx_el(n1) &
+ + q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
+ fy_el(n1)= fy_el(n1) &
+ + q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
+ fz_el(n1)= fz_el(n1) &
+ + q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
+
+ partmass=4._pr/3._pr*pi*rhop*partn(n2)*radp(n2)**3
+
+ fx_el(n2)= fx_el(n2) &
+ - q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
+ fy_el(n2)= fy_el(n2) &
+ - q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
+ fz_el(n2)= fz_el(n2) &
+ - q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass
return
@@ -271,17 +283,19 @@
!> @brief compute particle-wall charging
subroutine chargeParticleWall(n,direction)
use var
- real(kind=pr) ::ut2,un2,alpha1,dtcontact,alpha2, &
- C,tau,A,qc,q_el_temp,qmax
+ real(kind=pr) :: ut2,un2,alpha1,AoAtot,dqp,random(2),random_normal
integer :: n,i,direction
+ character(70) :: filename
- wcollnum(n)=wcollnum(n)+1
! pipe:
! ut2=up(n)**2
! un2=vp(n)**2+wp(n)**2
-
- if (direction.eq.2) then
+!
+ if (direction.eq.1) then
+ ut2=vp(n)**2+wp(n)**2
+ un2=up(n)**2
+ elseif (direction.eq.2) then
ut2=up(n)**2+wp(n)**2
un2=vp(n)**2
elseif (direction.eq.3) then
@@ -289,88 +303,63 @@
un2=wp(n)**2
endif
-
-!! this is the expression from Kolniak (1989)
-! alpha1=0.864_pr*2._pr*radp(n)*(un2*(ut2+un2))**(0.2_pr)* &
-! (rhop*((1._pr-nyw**(2._pr))/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr)
-!
-!! this is the expression from John (1979)
-! alpha2=(0.625_pr*pi*rhop*(1._pr+restRatio)*(un2+ut2)* &
+! John (1979)
+! alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*un2* &
! ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr)
-!
-! alpha1=alpha2
-!
-! A=pi*radp(n)**2*alpha1
-!
-! dtcontact=2.94_pr*alpha1*radp(n)/ut2**(0.5_pr)
-!
-! C=eps_el*A/h0
-!
-! tau=epsw*eps_el*r_ep
-!! tau=eps_el*r_ep
-!
-!! qmax=4.1038e-9*pi*(2*radp(n))**(1.7)
-!
-! qc=C*Vc*(1._pr-exp(-dtcontact/tau))
-! q_el_temp=qc+0.25_pr*alpha1*q_el(n)
-!! switch off wall charging
-! q_el_temp = 0._pr
-!
-! q_el(n)=-q_el_temp+q_el(n)
+! AoAtot=alpha1/4._pr
+! dqp=Qaccfactor*AoAtot*(qpmax-q_el(n))
- open(access='append',unit=10,file='results/partWall.dat')
- write(10,77) sqrt(ut2),sqrt(un2)
-77 format(2(es11.2e2,x))
+ call random_number(random)
+ random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method
+ dqp=Qaccfactor*(qpmax-q_el(n))*(random_normal+1._pr) ! ...(mu*random_normal+sigma)
+ q_el(n)=max(min(q_el(n)+dqp,qpmax),qp0)
+
+
+ write(filename,'(a,i3.3,a)') 'results/particlesWall_p',myid,'.dat'
+ open(access='append',unit=10,file=filename)
+ write(10,77) nt,t,sqrt(ut2),sqrt(un2),dqp,q_el(n)
+77 format(i7,5(x,es11.3e2))
close(10)
return
end
-!!####################################################################
+!####################################################################
!> @author Holger Grosshans
!> @brief compute particle-particle charging
-! subroutine chargeParticleParticle(n1,n2)
-! use var
-! real(kind=pr) uax2,urad2,alpha1,dtcontact,alpha2, &
-! c1,c2,tau,a,qc,q_el_temp,qmax
-! integer :: n1,n2,i
-!
-! uax2=(up(n1)-up(n2))**2
-! urad2=(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
+ subroutine chargeParticleParticle(n1,n2)
+ use var
+ real(kind=pr) urel2,alpha1,dqp,AoAtot,random(2),random_normal
+ integer :: n1,n2
+ character(70) :: filename
+
+! urel2=(up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2
!
-! if(urad2.eq.0.) then
-! qc=0.5*(q_el(n2)-q_el(n1))
-! goto 10
+! if (urel2.eq.0._pr) then
+! dqp=0.5_pr*(q_el(n2)-q_el(n1))
+! else
+!! Soo (1971)
+! alpha1=radp(n1)*radp(n2)*(0.625_pr*pi*rhop*(1._pr+restRatio)*urel2 &
+! *(radp(n1)+radp(n2))**(0.5_pr)/(radp(n1)**3+radp(n2)**3) &
+! *(1._pr-nyp**2)/Ep)**(0.4_pr)
+! dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1))
! endif
-!
-!! this is the expression from Soo (1971)
-! alpha1=radp(n1)*radp(n2)*(0.625*pi*rhop*(1+restratio)*(uax2+urad2) &
-! *(radp(n1)+radp(n2))**(0.5)/(radp(n1)**3+radp(n2)**3) &
-! *(1-nyp**(2.0))/ep)**(0.4)
-!
-! a=pi*radp(n1)*radp(n2)*alpha1/(radp(n1)+radp(n2))
-!
-! dtcontact=2.94*alpha1/urad2**(0.5)
-! c1=4*pi*eps_el*radp(n1)
-! c2=4*pi*eps_el*radp(n2)
-!
-! tau=c1*c2/(c1+c2)*(radp(n1)+radp(n2))/a*r_ep
-!
-!! qmax=4.1038e-9*pi*(2*radp(n))**(1.7)
-!
-! qc= c1*c2/(c1+c2)*(q_el(n2)/c2-q_el(n1)/c1) &
-! *(1.0-exp(-dtcontact/tau))
-!
-!
-! open(access='append',unit=10,file='results/partPart.dat')
-! write(10,77) q_el(n1),q_el(n2),qc,a,dtcontact
-!77 format(5(es11.2e2,x))
-! close(10)
-!
-!10 q_el(n1)=q_el(n1)+qc
-! q_el(n2)=q_el(n2)-qc
-!
-!
-! return
-! end
+! dqp=Qaccfactor*AoAtot*(q_el(n2)-q_el(n1))
+
+ call random_number(random)
+ random_normal=cos(2._pr*pi*random(1))*sqrt(-2._pr*log(random(2)+1.e-12_pr)) ! Box-Muller method
+ dqp=Qaccfactor*(q_el(n2)-q_el(n1))*(random_normal+1._pr)
+ q_el(n1)=max(min(q_el(n1)+dqp,qpmax),qp0)
+ q_el(n2)=max(min(q_el(n2)-dqp,qpmax),qp0)
+
+
+ write(filename,'(a,i3.3,a)') 'results/particlesPart_p',myid,'.dat'
+ open(access='append',unit=10,file=filename)
+ write(10,77) nt,t,q_el(n1),q_el(n2),dqp
+77 format(i7,4(x,es11.3e2))
+ close(10)
+
+
+ return
+ end
diff --git a/src/fluid.f90 b/src/fluid.f90
index efd9b2a079f634ad6e1be162b82e4b67f0e67e29..b96194da4ec6bb371ce7436e417f2dcec4f4f04b 100644
--- a/src/fluid.f90
+++ b/src/fluid.f90
@@ -9,6 +9,7 @@
call ddt
call deferredCorrection
+ if ((bcx.eq.'i').and.(myid.eq.0)) call inflow
do 1 it=1,itmax
call momentum(err1)
@@ -41,15 +42,18 @@
!> @brief calculate the adaptive time-step size
subroutine timestep
use var
+ real(kind=pr) syncMax,umax
- dt02=dt01
- dt01=dt
- dt= cfl*minval(hx)/maxval(u)
- call syncMax(dt)
+ umax= syncMax(maxval(abs(u)))
- if (nt.eq.1) then
+ if ((nt.eq.1).or.(umax.eq.0._pr)) then
+ dt= cfl*dimx/dimi/10._pr
dt01=dt
dt02=dt01
+ else
+ dt02=dt01
+ dt01=dt
+ dt= cfl*dimx/dimi/umax
endif
! fixed time-step:
diff --git a/src/main.f90 b/src/main.f90
index 47f26e6b2938d2819d8ff28d6c8bbce1a01196a9..91feab50c39d2725fe4808f8a67d6f2f93c9fa28 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -3,11 +3,11 @@
program pafiX
use var
use mpi
+ real(kind=pr) :: syncAv
call initParallel
if (myid.eq.0) write(*,'(/a/a)') repeat('#',72),version
-
if (myid.eq.0) write(*,'(a/a)') repeat('#',72),'Pre-processing'
call preProcessing
@@ -20,16 +20,16 @@
call timestep
call solveFluid
-
call solveElectrostatics
call solveParticles
+
call postProcessing
- call syncAv(timecom)
+ timecom(1)=syncAv(timecom)
if (myid.eq.0) write(*,'(3(a,es8.2e2))') &
'time (s) |physical= ',t, &
' |comput. = ',(mpi_wtime()-timebeg), &
- ' |MPI com = ',timecom
+ ' |MPI com = ',timecom(1)
10 enddo
diff --git a/src/makefile b/src/makefile
index 0d713b579b4d72f9e76ab1b7c9310120be59fa39..6ed0367b799943482043167ac24eec15e7df1821 100644
--- a/src/makefile
+++ b/src/makefile
@@ -3,29 +3,24 @@
OBJ = var.o main.o pre.o bc.o restart.o fluid.o \
momentum.o post.o schemes.o \
pressure.o writevtk_fluid_xz.o \
- writevtk_fluid_yz.o writedat_fluid.o \
parallel.o writevtk_grid.o particles.o \
mom1.o mom5.o particlesTransport.o \
mass.o electrostatics.o \
- writevtk_particles.o writedat_particles.o
+ writevtk_particles.o writedat_particles.o \
+ writevtk_fluid_xy.o writevtk_fluid_xyz.o \
+ writevtk_fluid_yz.o writedat_fluid_xz.o \
+ writedat_fluid_xy.o
INC =
-#O0: debugging, O3: optimization
-#dynamic linking (currently also for release):
# gcc:
CMP = mpifort
-# intel:
-#CMP = mpif90
-#FLAGS = -O3 -mcmodel=medium
+#O3: optimization
+FLAGS = -O3 -mcmodel=medium
#debugging:
-FLAGS = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
+#FLAGS = -g3 -O0 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -fno-tree-vectorize
FFLAGS = -c $(FLAGS)
-#static linking:
-#CMP = /home/holger/fortran/static/bin/mpifort
-#FLAGS = -O3 -fbounds-check -mcmodel=medium -fimplicit-none -fcheck=all -fbacktrace -floop-nest-optimize -ffpe-trap=invalid,zero,overflow -Wconversion -static
-#FFLAGS = -c $(FLAGS)
newfu: $(OBJ)
$(CMP) $(FLAGS) -o pafiX $(OBJ)
@@ -49,12 +44,18 @@ pressure.o: pressure.f90 $(INC)
$(CMP) $(FFLAGS) pressure.f90
mass.o: mass.f90 $(INC)
$(CMP) $(FFLAGS) mass.f90
+writevtk_fluid_xy.o: writevtk_fluid_xy.f90 $(INC)
+ $(CMP) $(FFLAGS) writevtk_fluid_xy.f90
writevtk_fluid_xz.o: writevtk_fluid_xz.f90 $(INC)
$(CMP) $(FFLAGS) writevtk_fluid_xz.f90
writevtk_fluid_yz.o: writevtk_fluid_yz.f90 $(INC)
$(CMP) $(FFLAGS) writevtk_fluid_yz.f90
-writedat_fluid.o: writedat_fluid.f90 $(INC)
- $(CMP) $(FFLAGS) writedat_fluid.f90
+writevtk_fluid_xyz.o: writevtk_fluid_xyz.f90 $(INC)
+ $(CMP) $(FFLAGS) writevtk_fluid_xyz.f90
+writedat_fluid_xy.o: writedat_fluid_xy.f90 $(INC)
+ $(CMP) $(FFLAGS) writedat_fluid_xy.f90
+writedat_fluid_xz.o: writedat_fluid_xz.f90 $(INC)
+ $(CMP) $(FFLAGS) writedat_fluid_xz.f90
writevtk_particles.o: writevtk_particles.f90 $(INC)
$(CMP) $(FFLAGS) writevtk_particles.f90
writedat_particles.o: writedat_particles.f90 $(INC)
diff --git a/src/momentum.f90 b/src/momentum.f90
index 56e611e7d41eae39f18df8a4068420fdf5f40f90..b1e95c167ea07ca90f369ab586d0d4fbacb54c6d 100644
--- a/src/momentum.f90
+++ b/src/momentum.f90
@@ -6,61 +6,57 @@
use var
real(kind=pr),dimension(ii,jj,ll) :: u1,v1,w1
- real(kind=pr) :: err
+ real(kind=pr) :: err,syncSum
real(kind=pr) :: dum,dvm,dwm
real(kind=pr) :: du,dv,dw,dua,dva,dwa
real(kind=pr) :: tu,tv,tw,qu,qv,qw
integer :: i,j,l
err=0._pr
- u1=0._pr; v1=0._pr; w1=0._pr
+ u1=u; v1=v; w1=w
dum=0._pr; dvm=0._pr; dwm=0._pr
- do 1 i=imin,imax,1
- do 1 j=jmin,jmax,1
- do 1 l=lmin,lmax,1
-
- if (celltype(i,j,l).ne.active) cycle
-
- dua=0._pr; dva=0._pr; dwa=0._pr
-
-! if (i.ne.imax) then
- call momx1(tu,qu,i,j,l)
- du=(defect_u(i,j,l)+tu)/qu*urfu
- u1(i,j,l)=u(i,j,l)+du
- dua=abs(du)
- if (dum.lt.dua) dum=dua
-! endif
-
- if (j.ne.jmax) then
- call momy1(tv,qv,i,j,l)
- dv=(defect_v(i,j,l)+tv)/qv*urfv
- v1(i,j,l)=v(i,j,l)+dv
- dva=abs(dv)
- if (dvm.lt.dva) dvm=dva
- endif
-
-
- if (l.ne.lmax) then
- call momz1(tw,qw,i,j,l)
- dw=(defect_w(i,j,l)+tw)/qw*urfw
- w1(i,j,l)=w(i,j,l)+dw
- dwa=abs(dw)
- if (dwm.lt.dwa) dwm=dwa
- endif
-
-
-12 err=err+dua*dua+dva*dva+dwa*dwa
-
-1 enddo
-
+ do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax
+
+ if (celltype(i,j,l).ne.active) cycle
+
+ dua=0._pr; dva=0._pr; dwa=0._pr
+
+ if (celltype(i+1,j,l).ne.wall) then
+ call momx1(tu,qu,i,j,l)
+ du=(defect_u(i,j,l)+tu)/qu*urfu
+ u1(i,j,l)=u(i,j,l)+du
+ dua=abs(du)
+ if (dum.lt.dua) dum=dua
+ endif
+
+ if (celltype(i,j+1,l).ne.wall) then
+ call momy1(tv,qv,i,j,l)
+ dv=(defect_v(i,j,l)+tv)/qv*urfv
+ v1(i,j,l)=v(i,j,l)+dv
+ dva=abs(dv)
+ if (dvm.lt.dva) dvm=dva
+ endif
+
+ if (celltype(i,j,l+1).ne.wall) then
+ call momz1(tw,qw,i,j,l)
+ dw=(defect_w(i,j,l)+tw)/qw*urfw
+ w1(i,j,l)=w(i,j,l)+dw
+ dwa=abs(dw)
+ if (dwm.lt.dwa) dwm=dwa
+ endif
+
+ err=err+dua*dua+dva*dva+dwa*dwa
+
+ enddo; enddo; enddo
u=u1; v=v1; w=w1
call sync(u); call sync(v); call sync(w)
- call bc
- call syncSum(err); err=(err/dimgptot/3._pr)**(0.5_pr)
+ call bcUVWP
+
+ err=(syncSum(err)/dimgptot/3._pr)**(0.5_pr)
! write(*,'(x,a,es9.2e2,a,3(es9.2e2))') &
! 'res. inner it. mom =',err,', max du,dv,dw =',dum,dvm,dwm
@@ -96,57 +92,52 @@
tu,tv,tw,qu,qv,qw,ra,raHO,tuHO,tvHO,twHO,quHO,qvHO,qwHO
integer :: i,j,l
-! velocity components are defined on ranges:
-! u(imin:imax, jmin:jmax, lmin:lmax)
-! v(imin:imax, jmin:jmax-1, lmin:lmax)
-! w(imin:imax, jmin:jmax, lmin:lmax-1)
-
- do 1 i=imin,imax
- do 1 j=jmin+1,jmax-1
- do 1 l=lmin+1,lmax-1
-
- if (celltype(i,j,l).ne.active) cycle
-
- if ((j-2.ge.jmin-1.and.j+1.le.jmax).and. &
- (l-2.ge.lmin-1.and.l+1.le.lmax)) then
- call mass2(ra,i,j,l)
- call mass4(raHO,i,j,l)
- defect_c(i,j,l)=raHO-ra
- else
- defect_c(i,j,l)=0._pr
- endif
+
+ do i=imin,imax+1; do j=jmin,jmax+1; do l=lmin,lmax+1
+
+ if ((celltype(i+1,j,l).ne.wall).and.(celltype(i-1,j,l).ne.wall).and. &
+ (celltype(i,j+1,l).ne.wall).and.(celltype(i,j-1,l).ne.wall).and. &
+ (celltype(i,j,l+1).ne.wall).and.(celltype(i,j,l-1).ne.wall)) then
+ call mass2(ra,i,j,l)
+ call mass4(raHO,i,j,l)
+ defect_c(i,j,l)=raHO-ra
+ else
+ defect_c(i,j,l)=0._pr
+ endif
- if ((j-3.ge.jmin.and.j+3.le.jmax).and. &
- (l-3.ge.lmin.and.l+3.le.lmax)) then
- call momx1(tu,qu,i,j,l)
-! call momx3(tuHO,quHO,i,j,l)
- call momx5(tuHO,quHO,i,j,l)
- defect_u(i,j,l)=tuHO-tu
- else
- defect_u(i,j,l)=0._pr
- endif
-
- if ((j-3.ge.jmin.and.j+3.le.jmax-1).and. &
- (l-3.ge.lmin.and.l+3.le.lmax)) then
- call momy1(tv,qv,i,j,l)
-! call momy3(tvHO,qvHO,i,j,l)
- call momy5(tvHO,qvHO,i,j,l)
- defect_v(i,j,l)=tvHO-tv
- else
- defect_v(i,j,l)=0._pr
- endif
-
- if ((j-3.ge.jmin.and.j+3.le.jmax).and. &
- (l-3.ge.lmin.and.l+3.le.lmax-1)) then
- call momz1(tw,qw,i,j,l)
-! call momz3(twHO,qwHO,i,j,l)
- call momz5(twHO,qwHO,i,j,l)
- defect_w(i,j,l)=twHO-tw
- else
- defect_w(i,j,l)=0._pr
- endif
-
-1 enddo
+ if (celltype(i,j,l).ne.active) cycle ! xmax+1 cells are only needed for defect_c
+
+ if ((celltype(i+3,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. &
+ (celltype(i,j+2,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. &
+ (celltype(i,j,l+2).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then
+ call momx1(tu,qu,i,j,l)
+ call momx5(tuHO,quHO,i,j,l)
+ defect_u(i,j,l)=tuHO-tu
+ else
+ defect_u(i,j,l)=0._pr
+ endif
+
+ if ((celltype(i+2,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. &
+ (celltype(i,j+3,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. &
+ (celltype(i,j,l+2).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then
+ call momy1(tv,qv,i,j,l)
+ call momy5(tvHO,qvHO,i,j,l)
+ defect_v(i,j,l)=tvHO-tv
+ else
+ defect_v(i,j,l)=0._pr
+ endif
+
+ if ((celltype(i+2,j,l).ne.wall).and.(celltype(i-2,j,l).ne.wall).and. &
+ (celltype(i,j+2,l).ne.wall).and.(celltype(i,j-2,l).ne.wall).and. &
+ (celltype(i,j,l+3).ne.wall).and.(celltype(i,j,l-2).ne.wall)) then
+ call momz1(tw,qw,i,j,l)
+ call momz5(twHO,qwHO,i,j,l)
+ defect_w(i,j,l)=twHO-tw
+ else
+ defect_w(i,j,l)=0._pr
+ endif
+
+ enddo; enddo; enddo
return
diff --git a/src/parallel.f90 b/src/parallel.f90
index 18867bb8b7c1dfc62b03f8fb51db33e02f1c749a..89947601574ea71f66fd3182662ad81e084e8d12 100644
--- a/src/parallel.f90
+++ b/src/parallel.f90
@@ -13,16 +13,8 @@
timecom=0._pr
timenow=real(mpi_wtime(),kind=pr)
timebeg=real(mpi_wtime(),kind=pr)
-
-
- next=myid+1
- prev=myid-1
- if(prev.eq.-1) prev=nrprocs-1 ! periodic bc in x-dir
- if(next.eq.nrprocs) next=0 ! periodic bc in x-dir
-! write(*,*) myid,'/',nrprocs,' Prev. ',prev,' Next ',next
-
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
@@ -30,7 +22,7 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief syncronize particles between processors
+!> @brief synchronize particles between processors
subroutine syncPart(np_sendl,np_sendr,np_recvl,np_recvr,myvar)
use var
use mpi
@@ -68,7 +60,8 @@
call mpi_wait(rrleft,mpistatus,mpierr)
call mpi_wait(rsright,mpistatus,mpierr)
call mpi_wait(rrright,mpistatus,mpierr)
-
+ if (next.eq.mpi_proc_null) n_recvr=0
+ if (prev.eq.mpi_proc_null) n_recvl=0
if (n_sendl.gt.0) call mpi_isend&
(p_sendl,n_sendl,mpi_pr,prev,2,mpi_comm_world,rsleft,mpierr)
@@ -96,7 +89,7 @@
3 enddo
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
@@ -141,6 +134,8 @@
call mpi_wait(rrleft,mpistatus,mpierr)
call mpi_wait(rsright,mpistatus,mpierr)
call mpi_wait(rrright,mpistatus,mpierr)
+ if (next.eq.mpi_proc_null) n_recvr=0
+ if (prev.eq.mpi_proc_null) n_recvl=0
if (n_sendl.gt.0) call mpi_isend&
@@ -169,7 +164,7 @@
3 enddo
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
@@ -177,7 +172,7 @@
!####################################################################
!> @author Holger Grosshans
-!> @brief syncronize fluid between processors
+!> @brief synchronize fluid between processors
subroutine sync(myvar)
use var
use mpi
@@ -190,18 +185,18 @@
timenow=real(mpi_wtime(),kind=pr)
n=0
- do 1 i = imin,imin+gc-1
- do 1 l = 1,ll; do 1 j = 1,jj
+ do i= imin,imin+gc-1
+ do l= 1,ll; do j= 1,jj
n=n+1
- sendleft(n) = myvar(i,j,l)
-1 enddo
+ sendleft(n)= myvar(i,j,l)
+ enddo; enddo; enddo
n=0
- do 2 i = ii-gc,ii-gc-gc+1,-1
- do 2 l = 1,ll; do 2 j = 1,jj
+ do i= imax,imax-gc+1,-1
+ do l= 1,ll; do j= 1,jj
n=n+1
- sendright(n) = myvar(i,j,l)
-2 enddo
+ sendright(n)= myvar(i,j,l)
+ enddo; enddo; enddo
call mpi_isend(sendleft,gc*jj*ll,mpi_pr,prev,1,mpi_comm_world,rsleft,mpierr)
call mpi_isend(sendright,gc*jj*ll,mpi_pr,next,2,mpi_comm_world,rsright,mpierr)
@@ -212,26 +207,26 @@
call mpi_wait(rrright,mpistatus,mpierr)
call mpi_wait(rrleft,mpistatus,mpierr)
- if(next.ne.mpi_proc_null) then
- n=0
- do 3 i = ii-gc+1,ii
- do 3 l = 1,ll; do 3 j = 1,jj
- n=n+1
- myvar(i,j,l) = recvright(n)
-3 enddo
+ if (next.ne.mpi_proc_null) then
+ n=0
+ do i= imax+1,ii
+ do l= 1,ll; do j= 1,jj
+ n=n+1
+ myvar(i,j,l)= recvright(n)
+ enddo; enddo; enddo
endif
- if(prev.ne.mpi_proc_null) then
- n=0
- do 4 i = gc,1,-1
- do 4 l = 1,ll; do 4 j = 1,jj
- n=n+1
- myvar(i,j,l) = recvleft(n)
-4 enddo
+ if (prev.ne.mpi_proc_null) then
+ n=0
+ do i= gc,1,-1
+ do l= 1,ll; do j= 1,jj
+ n=n+1
+ myvar(i,j,l)= recvleft(n)
+ enddo; enddo; enddo
endif
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
@@ -240,7 +235,7 @@
!> @author Holger Grosshans
!> @brief compute the max of a scalar over all processors
!> @param myvar scalar
- subroutine syncMax(myvar)
+ real(kind=pr) function syncMax(myvar)
use var
use mpi
real(kind=pr) :: myvar,myvar2
@@ -249,40 +244,11 @@
timenow=real(mpi_wtime(),kind=pr)
! mpi_allreduce does not work with intel compiler
-! call mpi_allreduce&
-! (myvar,myvar,1,mpi_pr,mpi_max,mpi_comm_world,mpierr)
-!
-! myvar=myvar/nrprocs
-
- if (myid.ne.0) then
- call mpi_isend(myvar,1,mpi_pr,0,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
- endif
-
- if (myid.eq.0) then
- do 1 proc=1,nrprocs-1
- call mpi_irecv(myvar2,1,mpi_pr,proc,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- myvar=max(myvar,myvar2)
-1 enddo
- endif
-
- if (myid.eq.0) then
- do 2 proc=1,nrprocs-1
- call mpi_isend(myvar,1,mpi_pr,proc,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
-2 enddo
- endif
-
- if (myid.ne.0) then
- call mpi_irecv(myvar,1,mpi_pr,0,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- endif
-
- myvar=myvar
-
+ call mpi_allreduce&
+ (myvar,syncMax,1,mpi_pr,mpi_max,mpi_comm_world,mpierr)
+
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
@@ -292,55 +258,31 @@
!> @author Holger Grosshans
!> @brief compute the sum of a scalar over all processors
!> @param myvar real scalar
- subroutine syncSum(myvar)
+ real(kind=pr) function syncSum(myvar)
use var
use mpi
- real(kind=pr) :: myvar,myvar2
+ real(kind=pr),intent(in) :: myvar
+ real(kind=pr) :: myvar2,sumvar
integer :: proc,rs,rr
timenow=real(mpi_wtime(),kind=pr)
! mpi_allreduce does not work with intel compiler
-! call mpi_allreduce&
-! (myvar,myvar,1,mpi_pr,mpi_sum,mpi_comm_world,mpierr)
-!
-! myvar=myvar/nrprocs
-
- if (myid.ne.0) then
- call mpi_isend(myvar,1,mpi_pr,0,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
- endif
-
- if (myid.eq.0) then
- do 1 proc=1,nrprocs-1
- call mpi_irecv(myvar2,1,mpi_pr,proc,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- myvar=myvar+myvar2
-1 enddo
- endif
-
- if (myid.eq.0) then
- do 2 proc=1,nrprocs-1
- call mpi_isend(myvar,1,mpi_pr,proc,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
-2 enddo
- endif
-
- if (myid.ne.0) then
- call mpi_irecv(myvar,1,mpi_pr,0,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- endif
+ call mpi_allreduce&
+ (myvar,syncSum,1,mpi_pr,mpi_sum,mpi_comm_world,mpierr)
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
- end
+ end function
+
+
!####################################################################
!> @author Holger Grosshans
!> @brief compute the sum of a scalar over all processors
!> @param myvar integer scalar
- subroutine syncSumI(myvar)
+ integer function syncSumI(myvar)
use var
use mpi
integer :: myvar,myvar2
@@ -349,38 +291,11 @@
timenow=real(mpi_wtime(),kind=pr)
! mpi_allreduce does not work with intel compiler
-! call mpi_allreduce&
-! (myvar,myvar,1,mpi_integer,mpi_sum,mpi_comm_world,mpierr)
-!
-! myvar=myvar/nrprocs
-
- if (myid.ne.0) then
- call mpi_isend(myvar,1,mpi_integer,0,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
- endif
-
- if (myid.eq.0) then
- do 1 proc=1,nrprocs-1
- call mpi_irecv(myvar2,1,mpi_integer,proc,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- myvar=myvar+myvar2
-1 enddo
- endif
-
- if (myid.eq.0) then
- do 2 proc=1,nrprocs-1
- call mpi_isend(myvar,1,mpi_integer,proc,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
-2 enddo
- endif
-
- if (myid.ne.0) then
- call mpi_irecv(myvar,1,mpi_integer,0,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- endif
+ call mpi_allreduce&
+ (myvar,syncSumI,1,mpi_integer,mpi_sum,mpi_comm_world,mpierr)
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
@@ -389,7 +304,7 @@
!> @author Holger Grosshans
!> @brief compute the average of a scalar over all processors
!> @param myvar scalar
- subroutine syncAv(myvar)
+ real(kind=pr) function syncAv(myvar)
use var
use mpi
real(kind=pr) :: myvar,myvar2
@@ -398,47 +313,20 @@
timenow=real(mpi_wtime(),kind=pr)
! mpi_allreduce does not work with intel compiler
-! call mpi_allreduce&
-! (myvar,myvar,1,mpi_pr,mpi_???,mpi_comm_world,mpierr)
-!
-! myvar=myvar/nrprocs
-
- if (myid.ne.0) then
- call mpi_isend(myvar,1,mpi_pr,0,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
- endif
-
- if (myid.eq.0) then
- do 1 proc=1,nrprocs-1
- call mpi_irecv(myvar2,1,mpi_pr,proc,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- myvar=myvar+myvar2
-1 enddo
- endif
-
- if (myid.eq.0) then
- do 2 proc=1,nrprocs-1
- call mpi_isend(myvar,1,mpi_pr,proc,0,mpi_comm_world,rs,mpierr)
- call mpi_wait(rs,mpistatus,mpierr)
-2 enddo
- endif
-
- if (myid.ne.0) then
- call mpi_irecv(myvar,1,mpi_pr,0,0,mpi_comm_world,rr,mpierr)
- call mpi_wait(rr,mpistatus,mpierr)
- endif
-
- myvar=myvar/nrprocs
+ call mpi_allreduce&
+ (myvar,myvar,1,mpi_pr,mpi_sum,mpi_comm_world,mpierr)
+
+ syncAv=myvar/nrprocs
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
!####################################################################
!> @author Holger Grosshans
-!> @brief checks if processors run syncronous
+!> @brief checks if processors run synchronous
subroutine syncCheck
use var
use mpi
@@ -451,7 +339,7 @@
call mpi_irecv(ntproc,1,mpi_integer,proc,0,mpi_comm_world,rr,mpierr)
call mpi_wait(rr,mpistatus,mpierr)
if (ntproc.ne.nt) then
- write(*,'(a)') 'Processors not syncronous!'
+ write(*,'(a)') 'Processors not synchronous!'
stop
endif
1 enddo
@@ -463,7 +351,7 @@
endif
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
@@ -472,7 +360,7 @@
!> @author Holger Grosshans
!> @brief generate m independent random numbers for each processor
!> @param m number of elements required per processor
-!> @return random random numbers
+!> @return random random numbers between 0 and 1
subroutine syncRandom(m,random)
use var
use mpi
@@ -496,7 +384,7 @@
endif
timeend=real(mpi_wtime(),kind=pr)
- timecom=timecom+timeend-timenow
+ timecom(1)=timecom(1)+timeend-timenow
return
end
diff --git a/src/particles.f90 b/src/particles.f90
index d9a1204dfcf23c867b21c29f242742fbcd5e1923..94637091ae48a66fd70195326cc111d4a58c433d 100644
--- a/src/particles.f90
+++ b/src/particles.f90
@@ -3,18 +3,19 @@
!> @brief solve particulate phase
subroutine solveParticles
use var
+ integer ip,jp,lp
- if (pnd.ne.0.or.np.ne.0) then
+ if (pnd.ne.0.or.npTot.ne.0) then
if (myid.eq.0) write(*,'(a)',advance='no') 'particle'
- if (nt.eq.ntinj) call partSeed
+ if (((bcx.eq.'i').and.(nt.ge.ntseed)).or.(nt.eq.ntseed)) call particlesSeed
call fluidVelocity
if (elForceScheme.ne.2) call forcesGauss
if (elForceScheme.ne.1) call forcesCoulomb
- call partVelocity
- call partCollide
- call partTransport !particles from neighbour gc added
+ call particlesVelocity
+ call particlesCollide
+ call particlesTransport !particles from neighbour gc added
call chargeDensity
- call momentumCoupling !particles from neighbour gc removed
+ call momentumCoupling !particles from neighbour gc removed
if (myid.eq.0) write(*,*)
endif
@@ -47,6 +48,7 @@
partn(m)=partn(n)
wcollnum(m)=wcollnum(n)
ppcollnum(m)=ppcollnum(n)
+ nGlob(m)=nGlob(n)
return
end
@@ -166,7 +168,7 @@
lend=lpa
endif
-! volume on the Eulerien grid
+! volume on the Eulerian grid
if (direction.eq.0) volE=(xf(iend)-xf(ibeg-1))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1))
if (direction.eq.1) volE=(xc(iend+1)-xc(ibeg))*(yf(jend)-yf(jbeg-1))*(zf(lend)-zf(lbeg-1))
if (direction.eq.2) volE=(xf(iend)-xf(ibeg-1))*(yc(jend+1)-yc(jbeg))*(zf(lend)-zf(lbeg-1))
@@ -186,9 +188,7 @@
deltaf2=(deltaf**2)/6._pr
!compute weights
- do 5 i=ibeg,iend
- do 5 j=jbeg,jend
- do 5 l=lbeg,lend
+ do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend
disx=xp(n)-xc(i)
if (direction.eq.1) disx=xp(n)-xf(i)
disy=yp(n)-yc(j)
@@ -198,7 +198,7 @@
dis2=disx*disx+disy*disy+disz*disz
weight(i+1-ibeg,j+1-jbeg,l+1-lbeg)=exp(-dis2/deltaf2)
-5 enddo
+ enddo; enddo; enddo
weight=weight/sum(weight) ! sum normalized
@@ -236,7 +236,7 @@
jbeg=jpa; jend=jpa
lbeg=lpa; lend=lpa
-! volume on the Eulerien grid
+! volume on the Eulerian grid
if (direction.eq.0) volE=(xf(ipa)-xf(ipa-1))*(yf(jpa)-yf(jpa-1))*(zf(lpa)-zf(lpa-1))
if (direction.eq.1) volE=(xc(ipa+1)-xc(ipa))*(yf(jpa)-yf(jpa-1))*(zf(lpa)-zf(lpa-1))
if (direction.eq.2) volE=(xf(ipa)-xf(ipa-1))*(yc(jpa+1)-yc(jpa))*(zf(lpa)-zf(lpa-1))
@@ -244,5 +244,115 @@
weight=1._pr
+ return
+ end
+
+!####################################################################
+!> @author Holger Grosshans
+!> @brief allocate particle arrays
+ subroutine allocateParticleArrays
+ use var
+ real(kind=pr), allocatable, dimension(:) :: &
+ tup,tvp,twp,tuf,tvf,twf,tuf01,tvf01,twf01,txp,typ,tzp,tradp, &
+ tq_el,tfx_el,tfy_el,tfz_el,tpartn,tnGlob
+ integer, allocatable, dimension(:) :: &
+ twcollnum,tppcollnum
+ integer :: syncSumI,nppTot
+
+ npTot=syncSumI(np)
+ nppTot=syncSumI(npp)
+ maxnp=npTot*(nrprocs)*3 !if nrprocs=1, one particle can be in domain, gc and sync storage
+
+ call allocateParticleArray(nppTot,up)
+ call allocateParticleArray(nppTot,vp)
+ call allocateParticleArray(nppTot,wp)
+ call allocateParticleArray(nppTot,uf)
+ call allocateParticleArray(nppTot,vf)
+ call allocateParticleArray(nppTot,wf)
+ call allocateParticleArray(nppTot,uf01)
+ call allocateParticleArray(nppTot,vf01)
+ call allocateParticleArray(nppTot,wf01)
+ call allocateParticleArray(nppTot,xp)
+ call allocateParticleArray(nppTot,yp)
+ call allocateParticleArray(nppTot,zp)
+ call allocateParticleArray(nppTot,radp)
+ call allocateParticleArray(nppTot,q_el)
+ call allocateParticleArray(nppTot,fx_el)
+ call allocateParticleArray(nppTot,fy_el)
+ call allocateParticleArray(nppTot,fz_el)
+ call allocateParticleArray(nppTot,partn)
+
+ call allocateParticleArrayI(nppTot,wcollnum)
+ call allocateParticleArrayI(nppTot,ppcollnum)
+ call allocateParticleArrayI(nppTot,nGlob)
+
+ contains
+
+ subroutine allocateParticleArray(nppTot,myvar)
+ use var
+ real(kind=pr), allocatable, dimension(:) :: &
+ myvar,tmyvar
+ integer :: nppTot
+
+ if (nppTot.ge.1) then
+ allocate(tmyvar(maxnp))
+ tmyvar(1:npp)=myvar(1:npp)
+ endif
+ if (allocated(myvar)) deallocate(myvar)
+ allocate(myvar(maxnp))
+ myvar=0._pr
+ if (nppTot.ge.1) then
+ myvar(1:npp)=tmyvar(1:npp)
+ endif
+
+ end
+
+ subroutine allocateParticleArrayI(nppTot,myvar)
+ use var
+ integer, allocatable, dimension(:) :: &
+ myvar,tmyvar
+ integer :: nppTot
+
+ if (nppTot.ge.1) then
+ allocate(tmyvar(maxnp))
+ tmyvar(1:npp)=myvar(1:npp)
+ endif
+ if (allocated(myvar)) deallocate(myvar)
+ allocate(myvar(maxnp))
+ myvar=0
+ if (nppTot.ge.1) then
+ myvar(1:npp)=tmyvar(1:npp)
+ endif
+
+ end
+
+ end
+
+
+!####################################################################
+!> @author Holger Grosshans
+!> @brief locate particles on Eulerian grid
+ subroutine particlesToCells
+ use var
+ integer :: i,j,l,n,ip,jp,lp
+
+ npic(:,:,:)=0
+
+ do 1 n=1,np
+ i=ip(n); j=jp(n); l=lp(n)
+ npic(i,j,l)=npic(i,j,l)+1
+1 enddo
+
+ if (allocated(nic)) deallocate(nic)
+ allocate(nic(ii,jj,ll,maxval(npic)))
+
+ nic=0
+
+ do 2 n=1,np
+ i=ip(n); j=jp(n); l=lp(n)
+ nic(i,j,l,count(nic(i,j,l,:).ne.0)+1)=n
+2 enddo
+
+
return
end
diff --git a/src/particlesTransport.f90 b/src/particlesTransport.f90
index 8562b17cfaec1e13f4e2367e83d795543dc510f3..ee4045f745c57ea0a1177b6e30eef1b14a728764 100644
--- a/src/particlesTransport.f90
+++ b/src/particlesTransport.f90
@@ -12,27 +12,18 @@
uf=0._pr; vf=0._pr; wf=0._pr
do 1 n=1,np
-
if (celltype(ip(n),jp(n),lp(n)).ne.active) cycle
-
do 4 direction=1,3
-
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction)
-
- do 5 i=ibeg,iend
- do 5 j=jbeg,jend
- do 5 l=lbeg,lend
- iw=i+1-ibeg
- jw=j+1-jbeg
- lw=l+1-lbeg
-
- if (direction.eq.1) uf(n)=uf(n)+weight(iw,jw,lw)*u(i,j,l)
- if (direction.eq.2) vf(n)=vf(n)+weight(iw,jw,lw)*v(i,j,l)
- if (direction.eq.3) wf(n)=wf(n)+weight(iw,jw,lw)*w(i,j,l)
-
-5 enddo
+ call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction)
+ do i=ibeg,iend; do j=jbeg,jend; do l=lbeg,lend
+ iw=i+1-ibeg
+ jw=j+1-jbeg
+ lw=l+1-lbeg
+ if (direction.eq.1) uf(n)=uf(n)+weight(iw,jw,lw)*u(i,j,l)
+ if (direction.eq.2) vf(n)=vf(n)+weight(iw,jw,lw)*v(i,j,l)
+ if (direction.eq.3) wf(n)=wf(n)+weight(iw,jw,lw)*w(i,j,l)
+ enddo; enddo; enddo
4 enddo
-
1 enddo
return
@@ -47,6 +38,7 @@
real(kind=pr) :: urel,Reyp,cd
integer :: n
+
urel= sqrt(( uf01(n)-up(n))*(uf01(n)-up(n)) &
+(vf01(n)-vp(n))*(vf01(n)-vp(n)) &
+(wf01(n)-wp(n))*(wf01(n)-wp(n)))
@@ -58,7 +50,8 @@
cd= 24._pr/Reyp * (1._pr + 1._pr/6._pr*Reyp**(2._pr/3._pr))
endif
drag= 0.375_pr*rhof*urel*cd/(rhop*radp(n))
-! drag=0._pr
+ !drag=0._pr
+
return
end
@@ -66,7 +59,7 @@
!####################################################################
!> @author Holger Grosshans
!> @brief explicit first order integration of particle velocities
- subroutine partVelocity
+ subroutine particlesVelocity
use var
real(kind=pr),dimension(3) :: f_fl,f_g
real(kind=pr) :: fact,dragdt,drag
@@ -103,17 +96,21 @@
!####################################################################
!> @author Holger Grosshans
!> @brief move particles
- subroutine partTransport
+ subroutine particlesTransport
use var
+ use mpi
real(kind=pr) :: xpold,ypold,zpold
integer, dimension (maxnp) :: &
np_sendr,np_sendl,np_recvl,np_recvr
integer :: n,m,i,j,l
+ integer :: syncSumI
+ character(70) :: filename
+ logical :: remove(maxnp)
+ remove=.false.
- m=np
- call syncSumI(m)
- if (myid.eq.0) write(*,'(a,i9)',advance='no') ' |transp. =',m
+ npTot=syncSumI(np)
+ if (myid.eq.0) write(*,'(a,i9)',advance='no') ' |transp. =',npTot
np_sendl=0
np_sendr=0
@@ -127,41 +124,75 @@
yp(n)=ypold+vp(n)*dt
zp(n)=zpold+wp(n)*dt
-! reflective bc's at side walls
- if ((yp(n).lt.ymin+radp(n)).or.(yp(n).gt.ymax-radp(n))) then
- vp(n)=-vp(n)*restRatio
- yp(n)=ypold+vp(n)*dt
- call chargeParticleWall(n,2)
+! boundary conditions
+ if (bcx.eq.'w') then
+ if (((myid.eq.0.).and.(xp(n).lt.xmin+radp(n))).or. &
+ ((myid.eq.nrprocs-1).and.(xp(n).gt.xmax-radp(n)))) then
+ up(n)=-up(n)*restRatio
+ xp(n)=xpold+up(n)*dt
+ if (qpmax.ne.qp0) call chargeParticleWall(n,1)
+ endif
+ elseif (bcx.eq.'i') then
+ if (myid.eq.nrprocs-1) then
+ if (xp(n).gt.xmax) then
+ remove(n)= .true.
+ else
+ remove(n)= .false.
+ endif
+ endif
+ endif
+
+ if (bcy.eq.'p') then
+ if (yp(n).lt.ymin) yp(n)=yp(n)+dimy
+ if (yp(n).gt.ymax) yp(n)=yp(n)-dimy
+ elseif (bcy.eq.'w') then
+ if ((yp(n).lt.ymin+radp(n)).or.(yp(n).gt.ymax-radp(n))) then
+ vp(n)=-vp(n)*restRatio
+ yp(n)=ypold+vp(n)*dt
+ wcollnum(n)=wcollnum(n)+1
+ if (qpmax.ne.qp0) call chargeParticleWall(n,2)
+ endif
endif
- if ((zp(n).lt.zmin+radp(n)).or.(zp(n).gt.zmax-radp(n))) then
- wp(n)=-wp(n)*restRatio
- zp(n)=zpold+wp(n)*dt
- call chargeParticleWall(n,3)
+ if (bcz.eq.'p') then
+ if (zp(n).lt.zmin) zp(n)=zp(n)+dimz
+ if (zp(n).gt.zmax) zp(n)=zp(n)-dimz
+ elseif (bcz.eq.'w') then
+ if ((zp(n).lt.zmin+radp(n)).or.(zp(n).gt.zmax-radp(n))) then
+ wp(n)=-wp(n)*restRatio
+ zp(n)=zpold+wp(n)*dt
+ wcollnum(n)=wcollnum(n)+1
+ if (qpmax.ne.qp0) call chargeParticleWall(n,3)
+ endif
endif
! particles are send (=1) if they are either on the active ghost
! cells or crossed the border to the next processor
- if ((xp(n).gt.xmin).and.(xp(n).lt.xf(gc+gc)).or. &
- (xpold.gt.xmin).and.(xp(n).le.xmin)) then
- np_sendl(n)=1
- else
- np_sendl(n)=0
- endif
- if ((xp(n).gt.xf(imax-gc)).and.(xp(n).lt.xmax).or. &
- (xpold.lt.xmax).and.(xp(n).ge.xmax)) then
- np_sendr(n)=1
- else
- np_sendr(n)=0
- endif
+ if (prev.ne.mpi_proc_null) then
+ if ((xp(n).gt.xmin).and.(xp(n).lt.xf(gc+gc)).or. &
+ (xpold.gt.xmin).and.(xp(n).le.xmin)) then
+ np_sendl(n)=1
+ else
+ np_sendl(n)=0
+ endif
+ endif
+ if (next.ne.mpi_proc_null) then
+ if ((xp(n).gt.xf(imax-gc)).and.(xp(n).lt.xmax).or. &
+ (xpold.lt.xmax).and.(xp(n).ge.xmax)) then
+ np_sendr(n)=1
+ else
+ np_sendr(n)=0
+ endif
+ endif
1 enddo
+
+
npp=np
np_recvl=0
np_recvr=0
-! synchronize
call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,xp)
call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,yp)
call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,zp)
@@ -173,17 +204,45 @@
call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,uf)
call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,vf)
call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,wf)
+ call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,uf01)
+ call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,vf01)
+ call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,wf01)
call syncPart(np_sendl,np_sendr,np_recvl,np_recvr,partn)
call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,wcollnum)
call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,ppcollnum)
+ call syncPartI(np_sendl,np_sendr,np_recvl,np_recvr,nGlob)
+
+
+ if ((bcx.eq.'p').and.(myid.eq.nrprocs-1)) then
+ do 2 n=npp+1,np
+ if (np_recvr(n).eq.1) xp(n)= xp(n) + dimxtot
+2 enddo
+ endif
+ if ((bcx.eq.'p').and.(myid.eq.0)) then
+ do 3 n=npp+1,np
+ if (np_recvl(n).eq.1) xp(n)= xp(n) - dimxtot
+3 enddo
+ endif
+
+ if ((bcx.eq.'i').and.(myid.eq.nrprocs-1)) then
+ m=0
+ do 4 n=1,np
+ if (remove(n)) then
+ write(filename,'(a,i3.3,a)') 'results/particlesOut_p',myid,'.dat'
+ open(access='append',unit=10,file=filename)
+ write(10,77) nt,t,q_el(n)
+77 format(i7,3(x,es11.2e2))
+ close(10)
+ else
+ m=m+1
+ call partN2M(n,m)
+ endif
+4 enddo
+ np=m
+ endif
+
+ call particlesToCells
-! add new particles
- do 2 n=npp+1,np
- if ((myid+1.eq.nrprocs).and.(np_recvr(n).eq.1)) &
- xp(n)= xp(n) + real(nrprocs,kind=pr)*dimx
- if ((myid.eq.0).and.(np_recvl(n).eq.1)) &
- xp(n)= xp(n) - real(nrprocs,kind=pr)*dimx
-2 enddo
return
end
@@ -209,7 +268,10 @@
do 1 n=1,np
- if (ip(n).lt.(imin-1).or.ip(n).gt.(imax+1)) cycle
+ if ((ip(n).lt.(imin-1)).or.(ip(n).gt.(imax+1)).or. &
+ (jp(n).lt.(jmin-1)).or.(jp(n).gt.(jmax+1)).or. &
+ (lp(n).lt.(lmin-1)).or.(lp(n).gt.(lmax+1))) cycle
+
! only drag forces
f_d(1)= drag(n)*(uf(n)-up(n))
@@ -219,31 +281,26 @@
volp= 4._pr/3._pr*pi*partn(n)*radp(n)**3
do 2 direction=1,3
-
- call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction)
+ call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,direction)
! compute source terms for the parcel
- do 5 l=lbeg,lend
- do 5 i=ibeg,iend
- do 5 j=jbeg,jend
- iw=i+1-ibeg
- jw=j+1-jbeg
- lw=l+1-lbeg
-
- if (direction.eq.1) then
- source= rhop/rhof * volp/volE * f_d(1)
- rui(i,j,l)= rui(i,j,l) + source*weight(iw,jw,lw)
- elseif (direction.eq.2) then
- source= rhop/rhof * volp/volE * f_d(2)
- rvi(i,j,l)= rvi(i,j,l) + source*weight(iw,jw,lw)
- elseif (direction.eq.3) then
- source= rhop/rhof * volp/volE * f_d(3)
- rwi(i,j,l)= rwi(i,j,l) + source*weight(iw,jw,lw)
- endif
+ do l=lbeg,lend; do i=ibeg,iend; do j=jbeg,jend
+ iw=i+1-ibeg
+ jw=j+1-jbeg
+ lw=l+1-lbeg
+ if (direction.eq.1) then
+ source= rhop/rhof * volp/volE * f_d(1)
+ rui(i,j,l)= rui(i,j,l) + source*weight(iw,jw,lw)
+ elseif (direction.eq.2) then
+ source= rhop/rhof * volp/volE * f_d(2)
+ rvi(i,j,l)= rvi(i,j,l) + source*weight(iw,jw,lw)
+ elseif (direction.eq.3) then
+ source= rhop/rhof * volp/volE * f_d(3)
+ rwi(i,j,l)= rwi(i,j,l) + source*weight(iw,jw,lw)
+ endif
+ enddo; enddo; enddo
-5 enddo
2 enddo
-
1 enddo
! hg not sure if I want this
@@ -259,14 +316,14 @@
! rwi_o(i,j,l)=rwi(io)
!6 enddo
- do 6 i=imin,imax; do 6 j=jmin,jmax; do 6 l=lmin,lmax
+ do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax
fsx(i,j,l)= -rui(i,j,l)
fsy(i,j,l)= -rvi(i,j,l)
fsz(i,j,l)= -rwi(i,j,l)
! fsx(i,j,l)= 0._pr
! fsy(i,j,l)= 0._pr
! fsz(i,j,l)= 0._pr
-6 enddo
+ enddo; enddo; enddo
! remove particles of another processor, they were only needed for source term
m=0
@@ -286,7 +343,7 @@
!####################################################################
!> @author Holger Grosshans
!> @brief compute particles collisions using ray casting
- subroutine partCollide
+ subroutine particlesCollide
use var
real(kind=pr),dimension(maxnp) :: xpc,ypc,zpc
real(kind=pr),dimension(3) :: relvelo,midlink,ncontact
@@ -295,18 +352,19 @@
r3n1,r3n2,rxx,projvtox,srel,projvton, &
closestrelvelo,fracdt, &
vex,vey,vez,vnx,vny,vnz
- integer :: n1,n2,numcol,ip,jp,lp,m
+ integer :: n1,n2,numcol,ip,jp,lp,m1,m2,i,j,l,syncSumI
numcol=0
- do 1 n1=2,np
- do 2 n2=1,(n1-1)
+ do i=1,ii; do j=1,jj; do l=1,ll
+! condition I only particles in same cell:
+ do 1 m1=2,npic(i,j,l)
+ do 2 m2=1,(m1-1)
+
+ n1=nic(i,j,l,m1)
+ n2=nic(i,j,l,m2)
-! condition I to check collision (same or neighbour cell):
- if ((ip(n1).lt.ip(n2)-1.or.ip(n1).gt.ip(n2)+1).or. &
- (jp(n1).lt.jp(n2)-1.or.jp(n1).gt.jp(n2)+1).or. &
- (lp(n1).lt.lp(n2)-1.or.lp(n1).gt.lp(n2)+1)) cycle
! setting the 2 particles in the rest frame of particle n2
relvelo(1)=up(n1)-up(n2)
relvelo(2)=vp(n1)-vp(n2)
@@ -324,13 +382,17 @@
! condition IIb to check collision (no rel velocity, no collision):
if (abs(srel).lt.1.e-19_pr) cycle
- partdist=dot_product(midlink,midlink)
+ partdist=sqrt(dot_product(midlink,midlink))
sumrad=radp(n1)+radp(n2)
! condition III to check collision (contact expected):
- if ((partdist-(projvtox/srel)**2).gt.sumrad**2) cycle
+ if ((partdist**2-(projvtox/srel)**2).gt.sumrad**2) then
+ cycle
+ elseif (partdist.le.sumrad) then ! if overlapping after seeding
+ cycle
+ endif
- closestrelvelo=sqrt(sumrad**2-(partdist-(projvtox/srel)**2))
+ closestrelvelo=sqrt(sumrad**2-(partdist**2-(projvtox/srel)**2))
fracdt=(projvtox/srel-closestrelvelo)/srel
! condition IV to check collision (contact during next timestep):
@@ -340,7 +402,7 @@
ppcollnum(n1)=ppcollnum(n1)+1
ppcollnum(n2)=ppcollnum(n2)+1
numcol=numcol+1
-! if (q_el(n1).ne.q_el(n2)) call chargeParticleParticle(n1,n2)
+ if (q_el(n1).ne.q_el(n2)) call chargeParticleParticle(n1,n2)
! fictitious contact point
xpc(n1)=xp(n1)+fracdt*up(n1)*dt
@@ -411,7 +473,9 @@
2 enddo
1 enddo
+ enddo; enddo; enddo
+ numcol=syncSumI(numcol)
if (myid.eq.0) write(*,'(a,i9)',advance='no') &
' |collide =',numcol
@@ -421,43 +485,62 @@
!####################################################################
!> @author Holger Grosshans
!> @brief seed particles
- subroutine partSeed
+ subroutine particlesSeed
use var
real(kind=pr),allocatable,dimension(:) :: randomx,randomy,randomz
- integer :: n,i,j,l,ip,jp,lp
+ integer :: n,i,j,l,ip,jp,lp,npseed,npseedTot, &
+ syncSumI
+ character(70) :: filename
-!> number of particles to be injected
- np=int(pnd*dimx*dimy*dimz)
-! np=2
-
- allocate(randomx(np),randomy(np),randomz(np))
+!> number of particles to be seeded
+ npp=np
+ if (bcx.eq.'i') then
+ if (myid.eq.0) then
+ npseed=int(dt*pnd)
+ else
+ npseed=0
+ endif
+ else
+ npseed=int(pnd*dimx*dimy*dimz)
+ ! npseed=2
+ endif
+ npseedTot=syncSumI(npseed)
+ np=npp+npseed
+ call allocateParticleArrays
if (myid.eq.0) write(*,'(a,i9)',advance='no') &
- ' |inject =',np*nrprocs
-
- open(status='replace',unit=10,file="results/injparticels.dat")
+ ' |seed =',npseedTot
- call syncRandom(np,randomx)
- call syncRandom(np,randomy)
- call syncRandom(np,randomz)
+ if (npseed.eq.0) return
- do 1 n=1,np
-!> seed not directly on wall or an active gc:
- xp(n)=randomx(n)*(xf(ii-gc)-xf(gc)-2._pr*prad)+xf(gc)+prad
- yp(n)=(2._pr*randomy(n)-1._pr)*(dimy/2._pr-2._pr*prad)
- zp(n)=(2._pr*randomz(n)-1._pr)*(dimz/2._pr-2._pr*prad)
+ write(filename,'(a,i3.3,a)') 'results/particlesSeed_p',myid,'.dat'
+ open(access='append',unit=10,file=filename)
-! if (mod(n,4).eq.0) yp(n)=-1.99e-2_pr
-! if (mod(n,4).eq.1) yp(n)=1.99e-2_pr
-! yp(n)=0_pr
-! if (mod(n,4).eq.2) zp(n)=-1.99e-2_pr
-! if (mod(n,4).eq.3) zp(n)=1.99e-2_pr
+ allocate(randomx(npseed),randomy(npseed),randomz(npseed))
+ if (bcx.eq.'i') then
+ if (myid.eq.0) then
+ call random_number(randomx)
+ call random_number(randomy)
+ call random_number(randomz)
+ endif
+ else
+ call syncRandom(npseed,randomx)
+ call syncRandom(npseed,randomy)
+ call syncRandom(npseed,randomz)
+ endif
-! yp(1)=-0.019_pr;zp(1)=-0.019_pr
-! yp(2)=0.019_pr;zp(2)=0.019_pr
-! vp(1)=-1._pr;vp(2)=1._pr;wp(1)=-1._pr;wp(2)=1._pr
+ do 1 n=npp+1,np
+!> seed not directly on wall or an active gc:
+ if (bcx.eq.'i') then
+ xp(n)=randomx(n-npp)*ubulk*dt+xf(imin-1)
+ else
+ xp(n)=randomx(n-npp)*(dimx-2._pr*prad)+xf(imin-1)+prad
+ endif
+ yp(n)=randomy(n-npp)*(dimy-2._pr*prad)+yf(jmin-1)+prad
+ zp(n)=randomz(n-npp)*(dimz-2._pr*prad)+zf(lmin-1)+prad
-! if(myid.eq.0) print*,xp(n),yp(n),zp(n),up(n),vp(n),wp(n)
+ nseedLoc=nseedLoc+1
+ nGlob(n)=nseedLoc*1000+myid
wcollnum(n)=0
ppcollnum(n)=0
@@ -465,20 +548,29 @@
i=ip(n)
j=jp(n)
l=lp(n)
+
up(n)=u(i,j,l)
if (xp(n).lt.xc(i)) up(n)=u(i-1,j,l)
vp(n)=v(i,j,l)
if (yp(n).lt.yc(j)) vp(n)=v(i,j-1,l)
wp(n)=w(i,j,l)
if (zp(n).lt.zc(l)) wp(n)=w(i,j,l-1)
+
+
+! xp(1)=0.02;yp(1)=0.;zp(1)=0.
+! xp(2)=0.03;yp(2)=0.;zp(2)=0.
+! up(1)=1.;up(2)=-1;vp(1)=0.;vp(2)=0.;wp(1)=0.;wp(2)=0.
+ uf(n)=up(n); vf(n)=vp(n); wf(n)=wp(n);
radp(n)=prad
partn(n)=1._pr
- q_el(n)=qelp
+ q_el(n)=qp0
write(10,'(3(es10.2e2,x))') radp(n),q_el(n),partn(n)
1 enddo
+
+ call particlesToCells
close(10)
diff --git a/src/post.f90 b/src/post.f90
index 7a25d037a87395fc8981d8512bb746e66cd26dad..515f085ff35d92529f748bf77858cc3941d8c478 100644
--- a/src/post.f90
+++ b/src/post.f90
@@ -6,13 +6,17 @@
if (mod(nt,int_results).eq.0.or.(nt).eq.1) then
call monitor
+ call writevtk_fluid_xy
call writevtk_fluid_xz
call writevtk_fluid_yz
- call writedat_ufluid_xz
- call writedat_vfluid_xz
- call writedat_wfluid_xz
+ call writedat_ufluid_xy
+ call writedat_vfluid_xy
+ call writedat_wfluid_xy
+! call writedat_ufluid_xz
+! call writedat_vfluid_xz
+! call writedat_wfluid_xz
! call writedat_coulomb
-! call writevtk_3d_fluid
+! call writevtk_fluid_xyz
call writevtk_particles
call writedat_particles
! call writevtk_fluid_av
@@ -31,11 +35,12 @@
!> @author Holger Grosshans
subroutine monitor
use var
- real(kind=pr) :: ucl,avu,tke,gradu, &
- Rec,Reb,avvp,avwp,avyp,avzp,avvf,avwf,avxp
- integer :: i,j,l,m
+ real(kind=pr) :: ucl,avu,avv,avw,tke,gradu,graduy,graduz, &
+ Rec,Reb,avvp,avwp,avyp,avzp,avvf,avwf,avxp,avqp, &
+ syncAv,syncSum,C01,A01
+ integer :: i,j,l,m,N01,syncSumI
logical :: file_ex
- character*70 :: filename2
+ character(70) :: filename2
@@ -43,7 +48,11 @@
ucl= sum(u(imin:imax,j,l))/(imax-imin+1)
- avu= sum(abs(u(imin:imax,jmin:jmax,lmin:lmax)) &
+ avu= sum(u(imin:imax,jmin:jmax,lmin:lmax) &
+ *vol(imin:imax,jmin:jmax,lmin:lmax)) / volTot
+ avv= sum(v(imin:imax,jmin:jmax,lmin:lmax) &
+ *vol(imin:imax,jmin:jmax,lmin:lmax)) / volTot
+ avw= sum(w(imin:imax,jmin:jmax,lmin:lmax) &
*vol(imin:imax,jmin:jmax,lmin:lmax)) / volTot
rmsv = sum((v(imin:imax,jmin:jmax,lmin:lmax))**2 &
@@ -53,54 +62,67 @@
*vol(imin:imax,jmin:jmax,lmin:lmax)) / volTot
rmsw=sqrt(rmsw)
- gradu= &
- (sum(u(imin:imax,jmin:jmax,lmin)*vol(imin:imax,jmin:jmax,lmin)) &
- +sum(u(imin:imax,jmin:jmax,lmax)*vol(imin:imax,jmin:jmax,lmax)) &
- +sum(u(imin:imax,jmin,lmin:lmax)*vol(imin:imax,jmin,lmin:lmax)) &
+ graduy= &
+ (sum(u(imin:imax,jmin,lmin:lmax)*vol(imin:imax,jmin,lmin:lmax)) &
+sum(u(imin:imax,jmax,lmin:lmax)*vol(imin:imax,jmax,lmin:lmax)))&
/ &
- (sum(vol(imin:imax,jmin:jmax,lmin)) &
- +sum(vol(imin:imax,jmin:jmax,lmax)) &
- +sum(vol(imin:imax,jmin,lmin:lmax)) &
+ (sum(vol(imin:imax,jmin,lmin:lmax)) &
+sum(vol(imin:imax,jmax,lmin:lmax)))&
/ (ymax-yc(jmax))
-
-
- if (np.gt.0) then
-! avvp=sum(vp(1:np))/np
-! avwp=sum(wp(1:np))/np
-! avxp=sum(xp(1:np))/np
- avyp=sum(yp(1:np))/np
- avzp=sum(zp(1:np))/np
-! avvf=sum(vf(1:np))/np
-! avwf=sum(wf(1:np))/np
+ graduz= &
+ (sum(u(imin:imax,jmin:jmax,lmin)*vol(imin:imax,jmin:jmax,lmin)) &
+ +sum(u(imin:imax,jmin:jmax,lmax)*vol(imin:imax,jmin:jmax,lmax))) &
+ / &
+ (sum(vol(imin:imax,jmin:jmax,lmin)) &
+ +sum(vol(imin:imax,jmin:jmax,lmax))) &
+ / (zmax-zc(lmax))
+
+ if (bcy.eq.'w'.and.bcz.eq.'w') gradu=(graduy+graduz)/2._pr
+ if (bcy.eq.'w'.and.bcz.eq.'p') gradu=graduy
+ if (bcy.eq.'p'.and.bcz.eq.'w') gradu=graduz
+
+ if (npTot.gt.0) then
+ avyp=syncSum(sum(yp(1:np)))/npTot
+ avzp=syncSum(sum(zp(1:np)))/npTot
+ avqp=syncSum(sum(q_el(1:np)))/npTot
+ if (bcy.eq.'w'.and.bcz.eq.'w') then
+ N01=count((min((yp(1:np)-ymin),(ymax-yp(1:np))).le.(dimy*0.01_pr)).or. &
+ (min((zp(1:np)-zmin),(zmax-zp(1:np))).le.(dimz*0.01_pr)))
+ A01=(dimy*dimz)-(dimy*0.98_pr)*(dimz*0.98_pr)
+ elseif (bcy.eq.'w'.and.bcz.eq.'p') then
+ N01=count(min((yp(1:np)-ymin),(ymax-yp(1:np))).le.(dimy*0.01_pr))
+ A01=(dimy*dimz)-(dimy*0.98_pr)*dimz
+ elseif (bcy.eq.'p'.and.bcz.eq.'w') then
+ N01=count(min((zp(1:np)-zmin),(zmax-zp(1:np))).le.(dimz*0.01_pr))
+ A01=(dimy*dimz)-dimy*(dimz*0.98_pr)
+ endif
+ C01=syncSumI(N01)/(A01*dimxTot)/pnd
else
- avyp=0
- avzp=0
+ npTot=0
+ avyp=0._pr
+ avzp=0._pr
+ avqp=0._pr
+ C01=0._pr
endif
-! call syncAv(avvp)
-! call syncAv(avwp)
- call syncAv(avyp)
- call syncAv(avzp)
-! call syncAv(avvf)
-! call syncAv(avwf)
- call syncAv(ucl)
- call syncAv(avu)
+ ucl=syncAv(ucl)
+ avu=syncAv(avu)
+ avv=syncAv(avv)
+ avw=syncAv(avw)
ubulk=avu
- call syncAv(rmsv)
- call syncAv(rmsw)
- call syncAv(gradu)
+ rmsv=syncAv(rmsv)
+ rmsw=syncAv(rmsw)
+ gradu=syncAv(gradu)
- if(myid.eq.0) then
+ if (myid.eq.0) then
Rec = ucl*delta/nuf
Reb = avu*delta/nuf
tau_w = muf*abs(gradu)
u_tau = sqrt(tau_w/rhof)
Re_tau = u_tau*delta/nuf
- delta_v= delta/Re_tau
+ if (Re_tau.gt.0._pr) delta_v= delta/Re_tau
filename2='output/monitor.dat'
@@ -110,16 +132,16 @@
else
open(10,file=filename2)
write(10,'(a,a)') '# ',version
- write(10,'(a,18(i10))') '#',(m,m=1,17)
- write(10,'(a,18(a10))') &
- '#','nt','t','dt','ucl','av(u)','rms(v)','rms(w)', &
- 'Rec','Reb','tau_w','u_tau','Re_tau', &
- 'delta_v','gradu','dpdx','av(yp)','av(zp)'
+ write(10,'(a,i10,i14,15(i10))') '#',(m,m=1,16)
+ write(10,'(a,a10,a14,14(a10))') &
+ '#','nt','t','dt','ucl','av(u)','av(v)','av(w)','rms(v)','rms(w)', &
+ 'Rec','Reb','tau_w','u_{tau}','Re_{tau}', &
+ 'C_{01}','av(qp)'
endif
- write(10,'(x,i10,21(es10.2e2))') &
- nt,t,dt,ucl,avu,rmsv,rmsw,Rec,Reb, &
- tau_w,u_tau,Re_tau,delta_v,gradu,dpdx,avyp,avzp
+ write(10,'(x,i10,es14.6e2,20(es10.2e2))') &
+ nt,t,dt,ucl,avu,avv,avw,rmsv,rmsw,Rec,Reb, &
+ tau_w,u_tau,Re_tau,C01,avqp
close(10)
diff --git a/src/pre.f90 b/src/pre.f90
index fb738e8df6ee0a103e38f291e3e3c437509d4c97..99c0b01c333ab965266224376ec35f7f5f7302e5 100644
--- a/src/pre.f90
+++ b/src/pre.f90
@@ -7,15 +7,14 @@
nt=0
call flowConditions
- call setParameter
call initVariables
call gridGeneration
call writevtk_grid
if (ntstart.gt.1) then
call readField
- do n=1,np
- q_el(n)=qelp !> possible variation of particle charge
- enddo
+! do n=1,np
+! q_el(n)=qelp !> possible variation of particle charge
+! enddo
! q_el(:)=qelp !> possible variation of particle charge
elseif (ntstart.eq.1) then
call initFlow
@@ -30,9 +29,10 @@
!> @brief read input data, compute flow conditions and write out
subroutine flowConditions
use var
- integer :: dimitotTemp,npTemp
- character*8 comp_date
- character*10 comp_time
+ use mpi
+ integer :: dimitotTemp
+ character(8) :: comp_date
+ character(10) :: comp_time
! read input file
@@ -43,28 +43,45 @@
read(10,*)
read(10,*) dimxtot,dimy,dimz
read(10,*)
+ read(10,*) bcx,bcy,bcz
+ read(10,*)
read(10,*) dimitotTemp,dimj,diml
read(10,*)
+ read(10,*) gridx,gridy,gridz
+ read(10,*)
read(10,*) ntstart,ntime
read(10,*)
read(10,*) int_results,int_restart
read(10,*)
read(10,*) cfl
read(10,*)
- read(10,*) ubulk,dpdx
+ read(10,*) ubulk0,dpdx
read(10,*)
read(10,*) rhof,nuf
read(10,*)
- read(10,*) pnd,ntinj
+ read(10,*) pnd,ntseed
read(10,*)
read(10,*) prad,rhop
read(10,*)
- read(10,*) qelp
+ read(10,*) qp0,qpmax
read(10,*)
read(10,*) g(1),g(2),g(3)
close(10)
+! connect processors
+ next=myid+1
+ prev=myid-1
+
+ select case (bcx)
+ case ('p')
+ if(prev.eq.-1) prev=nrprocs-1
+ if(next.eq.nrprocs) next=0
+ case ('w','i')
+ if(prev.eq.-1) prev=mpi_proc_null
+ if(next.eq.nrprocs) next=mpi_proc_null
+ end select
+
! compute flow conditions
ntend=ntstart+ntime-1
@@ -74,10 +91,14 @@
if (myid.eq.0.and.dimitot.ne.dimitotTemp) &
write(*,'(a,i8)') 'dimi set to ',dimitot
- npTemp=int(pnd*dimxtot*dimy*dimz)
+ npTot=int(pnd*dimxtot*dimy*dimz)
muf=nuf*rhof
+ ubulk=ubulk0
+
+ if (bcy.eq.'w'.and.bcz.eq.'w') delta= min(dimy/2._pr,dimz/2._pr)
+ if (bcy.eq.'w'.and.bcz.eq.'p') delta= dimy/2._pr
+ if (bcy.eq.'p'.and.bcz.eq.'w') delta= dimz/2._pr
- delta= dimy/2._pr
Re = ubulk*delta/nuf
ftt = dimxtot/ubulk
@@ -96,27 +117,34 @@
write(20,'(a,3(es11.2e2,x))') &
'dimensions in x,y,z-direction (m,m,m) = ', &
dimxtot,dimy,dimz
+ write(20,'(a,3x,3(a,x))') &
+ 'BCs in x,y,z-direction [(w)all/(p)eriodic] = ',bcx,bcy,bcz
write(20,'(a,3(i5),a,x,i8)') &
'Nr. of cells in x,y,z-direction (-,-,-) = ', &
dimitot,dimj,diml,', total =',dimitot*dimj*diml
write(20,'(a,3(i5),a,x,i8)') &
'Nr. of cells in x,y,z-direction per proc. (-,-,-) = ', &
dimi,dimj,diml,', total =',dimi*dimj*diml
+ write(20,'(a,3x,3(a,x))') &
+ 'Grid in x,y,z-direction [(u)niform/(s)tretch] = ',gridx,gridy,gridz
write(20,'(a,4(i8))') &
'nt: start, compute, results, restart (-,-,-,-) = ', &
ntstart,ntime,int_results,int_restart
write(20,'(a,es11.2e2)') &
'CFL number (-) = ', cfl
write(20,'(a,2(es11.2e2))') &
- 'initial bulk flow (m/s), pressure gradient (N/m**3) = ', ubulk,dpdx
+ 'in/initial flow (m/s), pressure gradient (N/m**3) = ', ubulk,dpdx
write(20,'(a,2(es11.2e2))') &
'fluid density (kg/m**3), kinematic visc. (m**2/s) = ', &
rhof,nuf
write(20,'(a,i8,es11.2,i8)') &
- 'particle number (-), density (-/m**3), inject nt (-) = ', npTemp,pnd,ntinj
+ 'particle number (-), (-/m**3) or (-/s), seed nt (-) = ', npTot,pnd,ntseed
write(20,'(a,3(es11.2e2))') &
- 'particle rad. (m), mat. density (kg/m**3), charge (C) = ', &
- prad,rhop,qelp
+ 'particle rad. (m), mat. density (kg/m**3) = ', &
+ prad,rhop
+ write(20,'(a,3(es11.2e2))') &
+ 'particle charge initial, maximum (C,C) = ', &
+ qp0,qpmax
write(20,'(a,3(es11.2e2))') &
'gravity vector in x,y,z-dir. (m/s**2,m/s**2,m/s**2) = ', &
g(1),g(2),g(3)
@@ -131,28 +159,6 @@
return
end
-!####################################################################
-!> @author Holger Grosshans
-!> @brief hardcode parameters which the user cannot change
- subroutine setParameter
- use var
-
- itmax=100
- tol= 1.e-2_pr
- urfu= 0.25_pr
- urfv= urfu
- urfw= urfv
- urfp= .5_pr
- gc= 3 !ghost cells between processors
- wc= 1 !wall boundary cells
- restRatio=1._pr
-! scheme to compute the electrostatic forces on particles
-! (see Grosshans, 2017): 1=Gauss, 2=Coulomb, 3=Hybrid
- elforceScheme=3
-
-
- return
- end
!####################################################################
!> @brief allocate and initialize all variables
@@ -160,12 +166,16 @@
use var
use mpi
- ii=dimi+2*gc; jj=dimj+2*wc; ll=diml+2*wc
- maxnp=int(pnd*dimx*dimy*dimz)*(nrprocs)*3 !if nrprocs=1, one particle can be in domain, gc and sync storage
+
+ ii=dimi+2*gc; jj=dimj+2*gc; ll=diml+2*gc
+ maxnp=0
+
+! parallel
+ allocate(mpistatus(mpi_status_size))
! grid
allocate( &
- hx(ii),hy(jj),hz(ll),celltype(ii,jj,ll), &
+ celltype(ii,jj,ll), &
xc(ii),yc(jj),zc(ll),xf(ii),yf(jj),zf(ll), &
dxcdi(ii),dycdj(jj),dzcdl(ll),dxfdi(ii),dyfdj(jj),dzfdl(ll), &
vol(ii,jj,ll))
@@ -182,25 +192,18 @@
dudt(ii,jj,ll),dvdt(ii,jj,ll),dwdt(ii,jj,ll))
! particles
- allocate( &
- up(maxnp),vp(maxnp),wp(maxnp), &
- uf(maxnp),vf(maxnp),wf(maxnp), &
- uf01(maxnp),vf01(maxnp),wf01(maxnp), &
- xp(maxnp),yp(maxnp),zp(maxnp), &
- radp(maxnp),q_el(maxnp))
- allocate(wcollnum(maxnp),ppcollnum(maxnp), &
- partn(maxnp), &
- fx_el(maxnp),fy_el(maxnp),fz_el(maxnp))
-
-! parallel
- allocate(mpistatus(mpi_status_size))
+ allocate(npic(ii,jj,ll))
+ np=0; npp=0; nseedLoc=0;
+ call allocateParticleArrays
+! init
u=0._pr; v=0._pr; w=0._pr; p=0._pr
dudt=0._pr; dvdt=0._pr; dwdt=0._pr
u01=u; v01=v; w01= w
fsx=0._pr; fsy=0._pr; fsz=0._pr
defect_c=0._pr; defect_u=0._pr; defect_v=0._pr; defect_w=0._pr
celltype=active
+ dxcdi=0._pr; dycdj=0._pr; dzcdl=0._pr; dxfdi=0._pr; dyfdj=0._pr; dzfdl=0._pr
t=0._pr
! electrostatics
@@ -221,14 +224,14 @@
!> @brief generate grid for gaseous phase
subroutine gridGeneration
use var
- real(kind=pr) :: scaley,scalez,fluxplu,fluxmin, &
+ real(kind=pr) :: hy_temp,hz_temp,fluxplu,fluxmin, &
dimi_pr,dimj_pr,diml_pr,gc_pr,jj_pr,ll_pr
- integer :: i,j,l,stretch,m
+ integer :: i,j,l,m
- imin=gc+1; imax=ii-gc
- jmin=wc+1; jmax=jj-wc
- lmin=wc+1; lmax=ll-wc
+ imin=1+gc; imax=ii-gc
+ jmin=1+gc; jmax=jj-gc
+ lmin=1+gc; lmax=ll-gc
gp=ii*jj*ll; dimgp=dimi*dimj*diml; dimgptot=(dimi*nrprocs)*dimj*diml
dimi_pr=real(dimi,kind=pr)
@@ -238,169 +241,174 @@
jj_pr=real(jj,kind=pr)
ll_pr=real(ll,kind=pr)
- stretch=1 !0:uniform grid, 1:stretched grid
-
- if (stretch.eq.0) then
+!> celltypes, default=active
+ do 1 i=1,ii
+ select case (bcx)
+ case ('p')
+ if ((i.lt.imin).or.(i.gt.imax)) celltype(i,:,:)= passive
+ case ('w','i')
+ if ((i.lt.imin).and.(myid.eq.0)) then
+ celltype(i,:,:)= wall
+ elseif ((i.lt.imin).and.(myid.ne.0)) then
+ celltype(i,:,:)= passive
+ elseif ((i.gt.imax).and.(myid.eq.nrprocs-1)) then
+ celltype(i,:,:)= wall
+ elseif ((i.gt.imax).and.(myid.ne.nrprocs-1)) then
+ celltype(i,:,:)= passive
+ endif
+ end select
+1 enddo
+ do 2 j=1,jj
+ select case (bcy)
+ case ('p')
+ if ((j.lt.jmin).or.(j.gt.jmax)) celltype(:,j,:)= passive
+ case ('w','i')
+ if ((j.lt.jmin).or.(j.gt.jmax)) celltype(:,j,:)= wall
+ end select
+2 enddo
+ do 3 l=1,ll
+ select case (bcz)
+ case ('p')
+ if ((l.lt.lmin).or.(l.gt.lmax)) celltype(:,:,l)= passive
+ case ('w','i')
+ if ((l.lt.lmin).or.(l.gt.lmax)) celltype(:,:,l)= wall
+ end select
+3 enddo
- hx(:)=dimx/dimi_pr
- hy(:)=dimy/dimj_pr
- hz(:)=dimz/diml_pr
+!> grid
do 4 i=1,ii
- if ((i.lt.imin).or.(i.gt.imax)) celltype(i,:,:)=ghostcell
- xc(i)=(real(i,kind=pr)-real(gc_pr)-0.5_pr)*hx(1)
- xf(i)=(real(i,kind=pr)-real(gc_pr))*hx(1)
+ if (gridx.eq.'u') then
+ xf(i)=real(i-gc,kind=pr)*dimx/dimi_pr
+ elseif (gridx.eq.'s') then
+ !not implemented yet
+ endif
+
+ if (i.eq.1) then
+ xc(i)= xf(1)-(dimx/dimi_pr)/2._pr
+ else
+ xc(i)= (xf(i)+xf(i-1))/2._pr
+ endif
4 enddo
+
do 5 j=1,jj
- yc(j)=(real(j,kind=pr)-real(jj_pr)/2._pr-0.5_pr)*hy(1)
- yf(j)=(real(j,kind=pr)-real(jj_pr)/2._pr)*hy(1)
- if (j.lt.jmin) then
- celltype(:,j,:)=wall
- elseif (j.gt.jmax) then
- celltype(:,j,:)=wall
- endif
+ if (gridy.eq.'u') then
+ hy_temp=dimy/dimj_pr
+ yf(j)=(real(j,kind=pr)-real(jj_pr)/2._pr)*hy_temp
+ elseif (gridy.eq.'s') then
+ hy_temp= (1._pr-cos(pi/dimj_pr))*dimy/2._pr !width of first active cell
+ if (j.lt.jmin) then
+ yf(j)= real(j+1-jmin,kind=pr)*hy_temp-dimy/2._pr
+ elseif (j.gt.jmax) then
+ yf(j)= real(j-jmax,kind=pr)*hy_temp+dimy/2._pr
+ else
+ yf(j)=-(cos(real(j-gc,kind=pr)*pi/dimj_pr))*dimy/2._pr
+ endif
+ endif
+
+ if (j.eq.1) then
+ yc(j)= yf(1)-hy_temp/2._pr
+ else
+ yc(j)= (yf(j)+yf(j-1))/2._pr
+ endif
5 enddo
- do 6 l=1,ll
- zc(l)=(real(l,kind=pr)-real(ll_pr)/2._pr-0.5_pr)*hz(1)
- zf(l)=(real(l,kind=pr)-real(ll_pr)/2._pr)*hz(1)
- if (l.lt.lmin) then
- celltype(:,:,l)=wall
- elseif (l.gt.lmax) then
- celltype(:,:,l)=wall
- endif
-6 enddo
- elseif (stretch.eq.1) then
+ do 6 l=1,ll
+ if (gridz.eq.'u') then
+ hz_temp=dimz/diml_pr
+ zf(l)=(real(l,kind=pr)-real(ll_pr)/2._pr)*hz_temp
+ elseif (gridz.eq.'s') then
+ hz_temp= (1._pr-cos(pi/diml_pr))*dimz/2._pr !width of first active cell
+ if (l.lt.lmin) then
+ zf(l)= real(l+1-lmin,kind=pr)*hz_temp-dimz/2._pr
+ elseif (l.gt.lmax) then
+ zf(l)= real(l-lmax,kind=pr)*hz_temp+dimz/2._pr
+ else
+ zf(l)=-(cos(real(l-gc,kind=pr)*pi/diml_pr))*dimz/2._pr
+ endif
+ endif
- hx(:)=dimx/dimi_pr
- scaley = dimy/2._pr
- scalez = dimz/2._pr
- do 1 i=1,ii
- if ((i.lt.imin).or.(i.gt.imax)) celltype(i,:,:)=ghostcell
- xc(i)=(real(i,kind=pr)-gc_pr-.5_pr)*hx(1)
- xf(i)=(real(i,kind=pr)-gc_pr)*hx(1)
-1 enddo
- do 2 j=1,jj
- if (j.lt.jmin) then
- celltype(:,j,:)=wall
- hy(j)= (cos(0._pr*pi/dimj_pr)-cos((1._pr)*pi/dimj_pr))*scaley
- yc(j)=(-(cos(0._pr*pi/dimj_pr)))*scaley+(real(j-wc,kind=pr)-.5_pr)*hy(j)
- yf(j)=(-(cos(0._pr*pi/dimj_pr)))*scaley+real(j-wc,kind=pr)*hy(j)
- elseif (j.gt.jmax) then
- celltype(:,j,:)=wall
- hy(j)= (cos(0._pr*pi/dimj_pr)-cos((1._pr)*pi/dimj_pr))*scaley
- yc(j)=(-(cos(real(jj-wc-wc-1,kind=pr)*pi/dimj_pr)))*scaley+(real(j-jj+wc,kind=pr)+.5_pr)*hy(j)
- yf(j)=(-(cos(real(jj-wc-wc-1,kind=pr)*pi/dimj_pr)))*scaley+(real(j-jj+wc,kind=pr)+1._pr)*hy(j)
- else
- hy(j)= (cos(real(j-wc-1,kind=pr)*pi/dimj_pr)-cos(real(j+1-wc-1,kind=pr)*pi/dimj_pr))*scaley
- yc(j)=(-(cos(real(j-wc-1,kind=pr)*pi/dimj_pr)))*scaley+.5_pr*hy(j)
- yf(j)=(-(cos(real(j-wc-1,kind=pr)*pi/dimj_pr)))*scaley+1._pr*hy(j)
- endif
-2 enddo
- do 3 l=1,ll
- if (l.lt.lmin) then
- celltype(:,:,l)=wall
- hz(l)= (cos((0._pr)*pi/diml_pr)-cos((1._pr)*pi/real(diml,kind=pr)))*scalez
- zc(l)=(-(cos((0._pr)*pi/diml_pr)))*scalez+(real(l-wc,kind=pr)-.5_pr)*hz(l)
- zf(l)=(-(cos((0._pr)*pi/diml_pr)))*scalez+real(l-wc,kind=pr)*hz(l)
- elseif (l.gt.lmax) then
- celltype(:,:,l)=wall
- hz(l)= (cos((0._pr)*pi/diml_pr)-cos((1._pr)*pi/diml_pr))*scalez
- zc(l)=(-(cos(real(ll-wc-wc-1,kind=pr)*pi/diml_pr)))*scalez+(real(l-ll+wc,kind=pr)+.5_pr)*hz(l)
- zf(l)=(-(cos(real(ll-wc-wc-1,kind=pr)*pi/diml_pr)))*scalez+(real(l-ll+wc,kind=pr)+1._pr)*hz(l)
- else
- hz(l)= (cos(real(l-wc-1,kind=pr)*pi/diml_pr)-cos(real(l+1-wc-1,kind=pr)*pi/diml_pr))*scalez
- zc(l)=(-(cos(real(l-wc-1,kind=pr)*pi/diml_pr)))*scalez+.5_pr*hz(l)
- zf(l)=(-(cos(real(l-wc-1,kind=pr)*pi/diml_pr)))*scalez+1._pr*hz(l)
- endif
-3 enddo
+ if (l.eq.1) then
+ zc(l)= zf(1)-hz_temp/2._pr
+ else
+ zc(l)= (zf(l)+zf(l-1))/2._pr
+ endif
+6 enddo
- endif !stretch
deltax=xf(ii-gc)*real(myid,kind=pr)
xc=xc+deltax
xf=xf+deltax
xmin=xf(gc); xmax=xf(ii-gc)
- ymin=yf(wc); ymax=yf(jj-wc)
- zmin=zf(wc); zmax=zf(ll-wc)
+ ymin=yf(gc); ymax=yf(jj-gc)
+ zmin=zf(gc); zmax=zf(ll-gc)
- do 7 i=1,ii
- do 7 j=1,jj
- do 7 l=1,ll
- vol(i,j,l)=hx(i)*hy(j)*hz(l)
-7 enddo
+ do i=2,ii; do j=2,jj; do l=2,ll
+ vol(i,j,l)=(xf(i)-xf(i-1))*(yf(j)-yf(j-1))*(zf(l)-zf(l-1))
+ enddo; enddo; enddo
volTot= sum(vol(imin:imax,jmin:jmax,lmin:lmax))
! compute derivatives of grid spacing for mapping
- do 8 i=2,ii-1
- if (i.ge.4.and.i.le.ii-3) then
- if (xc(i).gt.0._pr) then
- call weno(fluxplu,xc(i-2),xc(i-1),xc(i),xc(i+1),xc(i+2))
- call weno(fluxmin,xc(i-3),xc(i-2),xc(i-1),xc(i),xc(i+1))
- else
- call weno(fluxplu,xc(i+3),xc(i+2),xc(i+1),xc(i),xc(i-1))
- call weno(fluxmin,xc(i+2),xc(i+1),xc(i),xc(i-1),xc(i-2))
- endif
- dxcdi(i)=fluxplu-fluxmin
- if (xf(i).gt.0._pr) then
- call weno(fluxplu,xf(i-2),xf(i-1),xf(i),xf(i+1),xf(i+2))
- call weno(fluxmin,xf(i-3),xf(i-2),xf(i-1),xf(i),xf(i+1))
- else
- call weno(fluxplu,xf(i+3),xf(i+2),xf(i+1),xf(i),xf(i-1))
- call weno(fluxmin,xf(i+2),xf(i+1),xf(i),xf(i-1),xf(i-2))
+ do 8 i=imin,imax
+ if ((celltype(i+2,jmin,lmin).ne.wall).and.(celltype(i-2,jmin,lmin).ne.wall)) then
+ if (xc(i).gt.0._pr) then
+ call weno(fluxplu,xc(i-2),xc(i-1),xc(i),xc(i+1),xc(i+2))
+ call weno(fluxmin,xc(i-3),xc(i-2),xc(i-1),xc(i),xc(i+1))
+ else
+ call weno(fluxplu,xc(i+3),xc(i+2),xc(i+1),xc(i),xc(i-1))
+ call weno(fluxmin,xc(i+2),xc(i+1),xc(i),xc(i-1),xc(i-2))
+ endif
+ dxcdi(i)=fluxplu-fluxmin
+ if (xf(i).gt.0._pr) then
+ call weno(fluxplu,xf(i-2),xf(i-1),xf(i),xf(i+1),xf(i+2))
+ call weno(fluxmin,xf(i-3),xf(i-2),xf(i-1),xf(i),xf(i+1))
+ else
+ call weno(fluxplu,xf(i+3),xf(i+2),xf(i+1),xf(i),xf(i-1))
+ call weno(fluxmin,xf(i+2),xf(i+1),xf(i),xf(i-1),xf(i-2))
+ endif
+ dxfdi(i)=fluxplu-fluxmin
endif
- dxfdi(i)=fluxplu-fluxmin
- else
- dxcdi(i)=xf(i)-xf(i-1)
- dxfdi(i)=xc(i+1)-xc(i)
- endif
8 enddo
- do 9 j=2,jj-1
- if (j.ge.4.and.j.le.jj-3) then
- if (yc(j).gt.0._pr) then
- call weno(fluxplu,yc(j-2),yc(j-1),yc(j),yc(j+1),yc(j+2))
- call weno(fluxmin,yc(j-3),yc(j-2),yc(j-1),yc(j),yc(j+1))
- else
- call weno(fluxplu,yc(j+3),yc(j+2),yc(j+1),yc(j),yc(j-1))
- call weno(fluxmin,yc(j+2),yc(j+1),yc(j),yc(j-1),yc(j-2))
- endif
- dycdj(j)=fluxplu-fluxmin
- if (yf(j).gt.0._pr) then
- call weno(fluxplu,yf(j-2),yf(j-1),yf(j),yf(j+1),yf(j+2))
- call weno(fluxmin,yf(j-3),yf(j-2),yf(j-1),yf(j),yf(j+1))
- else
- call weno(fluxplu,yf(j+3),yf(j+2),yf(j+1),yf(j),yf(j-1))
- call weno(fluxmin,yf(j+2),yf(j+1),yf(j),yf(j-1),yf(j-2))
+ do 9 j=jmin,jmax
+ if ((celltype(imin,j+2,lmin).ne.wall).and.(celltype(imin,j-2,lmin).ne.wall)) then
+ if (yc(j).gt.0._pr) then
+ call weno(fluxplu,yc(j-2),yc(j-1),yc(j),yc(j+1),yc(j+2))
+ call weno(fluxmin,yc(j-3),yc(j-2),yc(j-1),yc(j),yc(j+1))
+ else
+ call weno(fluxplu,yc(j+3),yc(j+2),yc(j+1),yc(j),yc(j-1))
+ call weno(fluxmin,yc(j+2),yc(j+1),yc(j),yc(j-1),yc(j-2))
+ endif
+ dycdj(j)=fluxplu-fluxmin
+ if (yf(j).gt.0._pr) then
+ call weno(fluxplu,yf(j-2),yf(j-1),yf(j),yf(j+1),yf(j+2))
+ call weno(fluxmin,yf(j-3),yf(j-2),yf(j-1),yf(j),yf(j+1))
+ else
+ call weno(fluxplu,yf(j+3),yf(j+2),yf(j+1),yf(j),yf(j-1))
+ call weno(fluxmin,yf(j+2),yf(j+1),yf(j),yf(j-1),yf(j-2))
+ endif
+ dyfdj(j)=fluxplu-fluxmin
endif
- dyfdj(j)=fluxplu-fluxmin
- else
- dycdj(j)=yf(j)-yf(j-1)
- dyfdj(j)=yc(j+1)-yc(j)
- endif
9 enddo
- do 10 l=2,ll-1
- if (l.ge.4.and.l.le.ll-3) then
- if (zc(l).gt.0._pr) then
- call weno(fluxplu,zc(l-2),zc(l-1),zc(l),zc(l+1),zc(l+2))
- call weno(fluxmin,zc(l-3),zc(l-2),zc(l-1),zc(l),zc(l+1))
- else
- call weno(fluxplu,zc(l+3),zc(l+2),zc(l+1),zc(l),zc(l-1))
- call weno(fluxmin,zc(l+2),zc(l+1),zc(l),zc(l-1),zc(l-2))
+ do 10 l=lmin,lmax
+ if ((celltype(imin,jmin,l+2).ne.wall).and.(celltype(imin,jmin,l-2).ne.wall)) then
+ if (zc(l).gt.0._pr) then
+ call weno(fluxplu,zc(l-2),zc(l-1),zc(l),zc(l+1),zc(l+2))
+ call weno(fluxmin,zc(l-3),zc(l-2),zc(l-1),zc(l),zc(l+1))
+ else
+ call weno(fluxplu,zc(l+3),zc(l+2),zc(l+1),zc(l),zc(l-1))
+ call weno(fluxmin,zc(l+2),zc(l+1),zc(l),zc(l-1),zc(l-2))
+ endif
+ dzcdl(l)=fluxplu-fluxmin
+ if (zf(l).gt.0._pr) then
+ call weno(fluxplu,zf(l-2),zf(l-1),zf(l),zf(l+1),zf(l+2))
+ call weno(fluxmin,zf(l-3),zf(l-2),zf(l-1),zf(l),zf(l+1))
+ else
+ call weno(fluxplu,zf(l+3),zf(l+2),zf(l+1),zf(l),zf(l-1))
+ call weno(fluxmin,zf(l+2),zf(l+1),zf(l),zf(l-1),zf(l-2))
+ endif
+ dzfdl(l)=fluxplu-fluxmin
endif
- dzcdl(l)=fluxplu-fluxmin
- if (zf(l).gt.0._pr) then
- call weno(fluxplu,zf(l-2),zf(l-1),zf(l),zf(l+1),zf(l+2))
- call weno(fluxmin,zf(l-3),zf(l-2),zf(l-1),zf(l),zf(l+1))
- else
- call weno(fluxplu,zf(l+3),zf(l+2),zf(l+1),zf(l),zf(l-1))
- call weno(fluxmin,zf(l+2),zf(l+1),zf(l),zf(l-1),zf(l-2))
- endif
- dzfdl(l)=fluxplu-fluxmin
- else
- dzcdl(l)=zf(l)-zf(l-1)
- dzfdl(l)=zc(l+1)-zc(l)
- endif
10 enddo
@@ -419,15 +427,13 @@
call syncRandom(gp,random)
- alpha=0.7_pr
+ alpha=0.5_pr
- delta_v=delta/150._pr
- utau=0.05_pr
+! delta_v=delta/150._pr
+! utau=0.05_pr
m=0
- do 1 i=imin,imax
- do 1 j=jmin,jmax
- do 1 l=lmin,lmax
+ do i=imin,imax; do j=jmin,jmax; do l=lmin,lmax
m=m+1
! law of the wall/log law:
@@ -446,15 +452,16 @@
! w(i,j,l)=wplus*utau
! linear:
- dist=delta-max(abs(yc(j)),abs(zc(l)))
- u(i,j,l)=3._pr*ubulk*dist/delta + (random(m)-.5_pr)*alpha*ubulk
-! u(i,j,l)=3._pr*ubulk*dist/delta
-1 enddo
+ if (bcy.eq.'w'.and.bcz.eq.'w') dist=delta-max(abs(yc(j)),abs(zc(l)))
+ if (bcy.eq.'w'.and.bcz.eq.'p') dist=delta-abs(yc(j))
+ if (bcy.eq.'p'.and.bcz.eq.'w') dist=delta-abs(zc(l))
+ u(i,j,l)=ubulk*dist/delta + (random(m)-.5_pr)*alpha*ubulk
+ enddo; enddo; enddo
v=0._pr
w=0._pr
u01=u; v01=v; w01=w
- call bc
+ call bcUVWP
call sync(u); call sync(v); call sync(w)
return
diff --git a/src/pressure.f90 b/src/pressure.f90
index 8a161979140b842f4487f54b11e87ccd3e372318..3224269d319321cbe2a3bc42407fddfa700afa40 100644
--- a/src/pressure.f90
+++ b/src/pressure.f90
@@ -1,16 +1,13 @@
!> @author Holger Grosshans
!> @brief one sweep to relax mass conservation using distributive
!> Jacobi method based on distributive Gauss-Seidel by Brandt (1972)
-!> @param direction 1: sweep in increasing index order
-!> -1: sweep in decreasing index order
!> @param err L2 error norm
subroutine pressure(err)
use var
real(kind=pr),dimension(ii,jj,ll) :: u1,v1,w1
- real(kind=pr) :: err,h2sum,de,dex,dey,dez,ra,dpa,zz,da,dpm, &
- erru,errv,errw, &
- Aw,Ap,Cp,rCp,dpa1,ua,va,wa
- integer :: i,j,l,n,nzz
+ real(kind=pr) :: err,h2sum,de,dex,dey,dez,ra,dpa,dpm, &
+ erru,errv,errw,syncSum
+ integer :: i,j,l
err=0._pr
@@ -21,65 +18,47 @@
u1=u; v1=v; w1=w
- do 1 i=imin,imax+1,1
- do 1 j=jmin,jmax,1
- do 1 l=lmin,lmax,1
-
- if (celltype(i,j,l).ne.active) cycle
-
- h2sum=0._pr
- if (celltype(i-1,j,l).eq.active) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2
- if (celltype(i+1,j,l).eq.active) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2
- if (celltype(i,j-1,l).eq.active) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2
- if (celltype(i,j+1,l).eq.active) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2
- if (celltype(i,j,l-1).eq.active) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2
- if (celltype(i,j,l+1).eq.active) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2
-
-
- call mass2(ra,i,j,l)
- de = (defect_c(i,j,l)+ra)/h2sum*urfp
- dex = de/(xf(i)-xf(i-1))
- dey = de/(yf(j)-yf(j-1))
- dez = de/(zf(l)-zf(l-1))
-
- if (i.ne.imax+1) u1(i,j,l) = u1(i,j,l)+dex
- if (i.ne.imin) u1(i-1,j,l) = u1(i-1,j,l)-dex
- if (j.ne.jmax) v1(i,j,l) = v1(i,j,l)+dey
- if (j.ne.jmin) v1(i,j-1,l) = v1(i,j-1,l)-dey
- if (l.ne.lmax) w1(i,j,l) = w1(i,j,l)+dez
- if (l.ne.lmin) w1(i,j,l-1) = w1(i,j,l-1)-dez
-
- dpa = de*(1._pr/dt+2._pr*nuf/((xf(i)-xf(i-1))**2+(yf(j)-yf(j-1))**2+(zf(l)-zf(l-1))**2))
-! ua = u(i,j,l)
-! va = 0.25_pr*(v(i,j,l)+v(i,j-1,l)+v(i+1,j,l)+v(i+1,j-1,l))
-! wa = 0.25_pr*(w(i,j,l)+w(i+1,j,l)+w(i,j,l-1)+w(i+1,j,l-1))
-
-! Aw=-nuf*rx2(i)
-! Ap=-1.5_pr/dt -abs(ua)*rx(i) -abs(va)*ry(j) -abs(wa)*rz(l) &
-! -2._pr*nuf*(rx2(i)+ry2(j)+rz2(l))
-! Cp=1._pr/rhof*rx(i)
-! rCp=rhof*hx(i)
-
-! dpa=(Aw-Ap)*rCp*dex
-
- p(i,j,l)= p(i,j,l)+dpa
-
- da=dpa
- if(dpm.lt.da) dpm=da
- erru=erru+dex*dex
- errv=errv+dey*dey
- errw=errw+dez*dez
- err=err+da*da
-
-
- 1 enddo
+ do i=imin,imax+1; do j=jmin,jmax+1; do l=lmin,lmax+1
+
+ if (celltype(i,j,l).eq.wall) cycle
+
+ h2sum=0._pr
+ if (celltype(i+1,j,l).ne.wall) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2
+ if (celltype(i-1,j,l).ne.wall) h2sum=h2sum+1._pr/(xf(i)-xf(i-1))**2
+ if (celltype(i,j+1,l).ne.wall) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2
+ if (celltype(i,j-1,l).ne.wall) h2sum=h2sum+1._pr/(yf(j)-yf(j-1))**2
+ if (celltype(i,j,l+1).ne.wall) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2
+ if (celltype(i,j,l-1).ne.wall) h2sum=h2sum+1._pr/(zf(l)-zf(l-1))**2
+
+ call mass2(ra,i,j,l)
+ de= (defect_c(i,j,l)+ra)/h2sum*urfp
+ dpa= de*rhof*(1._pr/dt+2._pr*nuf/((xf(i)-xf(i-1))**2+(yf(j)-yf(j-1))**2+(zf(l)-zf(l-1))**2))
+ dex= de/(xf(i)-xf(i-1))
+ dey= de/(yf(j)-yf(j-1))
+ dez= de/(zf(l)-zf(l-1))
+
+ if (celltype(i,j,l).eq.active) p(i,j,l)= p(i,j,l)+dpa
+ if (celltype(i,j,l).eq.active.and.celltype(i+1,j,l).ne.wall) u1(i,j,l)= u1(i,j,l)+dex
+ if (celltype(i-1,j,l).eq.active.and.celltype(i,j,l).ne.wall) u1(i-1,j,l)= u1(i-1,j,l)-dex
+ if (celltype(i,j,l).eq.active.and.celltype(i,j+1,l).ne.wall) v1(i,j,l)= v1(i,j,l)+dey
+ if (celltype(i,j-1,l).eq.active.and.celltype(i,j,l).ne.wall) v1(i,j-1,l)= v1(i,j-1,l)-dey
+ if (celltype(i,j,l).eq.active.and.celltype(i,j,l+1).ne.wall) w1(i,j,l)= w1(i,j,l)+dez
+ if (celltype(i,j,l-1).eq.active.and.celltype(i,j,l).ne.wall) w1(i,j,l-1)= w1(i,j,l-1)-dez
+
+ if(dpm.lt.dpa) dpm=dpa
+ erru=erru+dex*dex
+ errv=errv+dey*dey
+ errw=errw+dez*dez
+ err=err+dpa*dpa
+
+ enddo; enddo; enddo
u=u1; v=v1; w=w1
call sync(u); call sync(v); call sync(w); call sync(p)
- call bc
+ call bcUVWP
- call syncSum(err); err=(err/dimgptot)**(0.5_pr)
+ err=(syncSum(err)/dimgptot)**(0.5_pr)
! write(*,'(x,a,4(es9.2e2))') 'res. inner it. pressure corr. =',err
return
diff --git a/src/restart.f90 b/src/restart.f90
index 184d55ec2f785f8a027d46742f2dac4ffdf74b1e..e6a8cbfe6a56e638692bf07eedf7515dbd2cf1a3 100644
--- a/src/restart.f90
+++ b/src/restart.f90
@@ -22,7 +22,7 @@
subroutine saveField
use var
integer :: i,j,l,n
- character*70 filename
+ character(70) :: filename
write(filename,'(a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,i6.6)') &
'restart/fluidField_p',myid,'_',dimi,'_',dimj,'_',diml,'_',nt
@@ -47,7 +47,7 @@
subroutine readField
use var
integer :: i,j,l,n
- character*70 filename
+ character(70) :: filename
if (myid.eq.0) write(*,'(a)') 'read old data'
@@ -62,13 +62,19 @@
open(unit=12,file=filename,form='unformatted',status='old',err=9000)
rewind(12)
read(12) np
+ call allocateParticleArrays
read(12) (up(n),vp(n),wp(n),xp(n),yp(n),zp(n),uf(n),vf(n),wf(n), &
radp(n),partn(n),q_el(n),wcollnum(n),ppcollnum(n),n=1,np)
close(12)
ntstart=ntstart+1
ntend=ntend+1
-
+
+ do n=1,np
+ nseedLoc=nseedLoc+1
+ nGlob(n)=nseedLoc*1000+myid
+ enddo
+ call particlesToCells
return
diff --git a/src/var.f90 b/src/var.f90
index c913e6f45e843f302673118a40efc490abf080d1..20090b91f7654e62ea741b2c32822a873134c33f 100644
--- a/src/var.f90
+++ b/src/var.f90
@@ -4,169 +4,156 @@
module var
implicit none
-! parameter
+! hardcoded parameters which the user cannot change
integer, parameter :: &
precision=8, & !< number precision
- pr=selected_real_kind(precision)
+ pr=selected_real_kind(precision), &
+ elforceScheme=3, & !< 1=Gauss, 2=Coulomb, 3=Hybrid (Grosshans&Papalexandris, 2017)
+ itmax=100, & !< max number of iterations for 1 equation
+ gc= 3 !< number of ghost cells, for consistence always 3, also at walls
+
real(kind=pr), parameter :: &
eps_el= 8.85e-12_pr, & !< permittivity vacuum/air
-! pi=3.141592653589793_pr
- pi=4._pr*atan(1._pr)
- character*70 :: version='pafiX v1.0.0 (Copyright 2019 by H. Grosshans)'
+ pi= 4._pr*atan(1._pr), &
+ tol= 1.e-2_pr, & !< factor by which L2 shall decrease
+ urfu= .25_pr, & !< under relaxation factor variable u
+ urfv= urfu, &
+ urfw= urfv, &
+ urfp= .5_pr, &
+ restRatio= 1._pr !< particle material restitution ratio
+ real(kind=pr), parameter :: &
+ Ew=1.e11_pr, & !< duct Young's modulus (kg/s**2/m)
+ Ep=1.e8_pr, & !< particle Young's modulus (kg/s**2/m)
+ nyp=0.4_pr , & !< particle Poisson ratio
+ nyw=0.28_pr, & !< duct Poisson ratio
+ Qaccfactor= 0.1_pr !< artificially accelerate the charging rate
+
+ character(70) :: version='pafiX v1.1.0 (Copyright 2019 by H. Grosshans)'
! input
real(kind=pr) :: &
- cfl, & !< CFL number (-)
- ubulk, & !< streamwise bulk fluid velocity (m/s)
- prad, & !< particle radius (m)
- qelp, & !< particle charge (C)
- dimx, & !< dimension of local domain in x-direction (m)
- dimy,dimz, &
- rhof, & !< fluid density (kg/m**3)
- nuf, & !< fluid kinematic viscosity (m**2/s)
- pnd, & !< particle number density (-/m**3)
- rhop, & !< particle material density (kg/m**3)
- g(3), & !< gravity vector (m/s**2,m/s**2,m/s**2)
- dimxtot, & !< dimension of domain (m)
- delta, & !< duct helf width
- tau_w, & !< wall shear stress
- u_tau, & !< friction velocity
- Re, & !< bulk Reynolds number
- Re_tau, & !< friction Reynolds number
- delta_v, & !< viscous length scale (m)
- t_tau, & !< viscous time scale (s)
- ftt, & !< flow through time (s)
+ cfl, & !< CFL number (-)
+ ubulk, ubulk0, & !< streamwise bulk fluid velocity, initial/inlet (m/s)
+ prad, & !< particle radius (m)
+ qp0,qpmax, & !< initial and maximum particle charge (C,C)
+ dimx,dimy,dimz, & !< dimension of local domain in x-direction (m)
+ rhof, & !< fluid density (kg/m**3)
+ nuf, & !< fluid kinematic viscosity (m**2/s)
+ pnd, & !< particle number density (-/m**3)
+ rhop, & !< particle material density (kg/m**3)
+ g(3), & !< gravity vector (m/s**2,m/s**2,m/s**2)
+ dimxtot, & !< dimension of domain (m)
+ delta, & !< duct half width
+ tau_w, & !< wall shear stress
+ u_tau, & !< friction velocity
+ Re, & !< bulk Reynolds number
+ Re_tau, & !< friction Reynolds number
+ delta_v, & !< viscous length scale (m)
+ t_tau, & !< viscous time scale (s)
+ ftt, & !< flow through time (s)
rmsv,rmsw
integer :: &
- ntstart, & !< start time step
- ntime, & !< end time step
- dimi, & !< local number of cells in x-direction
- dimj,diml, &
- int_results, & ! intervall write out results
- int_restart, & ! intervall write out restart files
- dimitot, & !< global number of cells in x-direction
- ntinj !< time step particle injection
+ ntstart, & !< start time step
+ ntime, & !< end time step
+ dimi,dimj,diml, & !< local number of cells in x-direction
+ int_results, & !< interval write out results
+ int_restart, & !< interval write out restart files
+ dimitot, & !< global number of cells in x-direction
+ ntseed !< time step particle seeding
! run
real(kind=pr) :: &
- dt, & !< size of time step nt (s)
- tol, & !< factor by which L2 shall decrease
- urfu, & !< under relaxation factor variable u
- urfv,urfw,urfp, &
- dpdx, & !< streamwise pressure gradient (N/m**3)
- dt01, & !< time step size nt-1
- dt02, & !< time step size nt-2
+ dt, & !< size of time step nt (s)
+ dpdx, & !< streamwise pressure gradient (N/m**3)
+ dt01, & !< time step size nt-1
+ dt02, & !< time step size nt-2
tau01,tau02, &
- t !< physical time (s)
- integer :: nt, & !< time step
- ntend, & !< last time step of computation
- elforceScheme, & !< 1=Gauss, 2=Coulomb, 3=Hybrid (Grosshans, 2017)
- itmax !< max number of iterations for 1 equation
+ t !< physical time (s)
+ integer :: nt, & !< time step
+ ntend !< last time step of computation
! grid
+ character(1) :: &
+ bcx,bcy,bcz, & !< type of boundary condition [(w)all/(p)eriodic]
+ gridx,gridy,gridz !< grid [(u)niform/(s)tretch]
real(kind=pr), allocatable, dimension(:,:,:) :: &
- vol !< volume per cell (m**3)
+ vol !< volume per cell (m**3)
real(kind=pr), allocatable, dimension(:) :: &
- hx, & !< side length of a cell (m)
- hy,hz, &
- xc, & !< center point of a cell (m)
- yc,zc, &
- xf, & !< face of a cell in pos. direction from xc (m)
- yf,zf, &
- dxcdi, & !< derivative of grid mapping at cell center (m/m)
- dycdj,dzcdl, &
- dxfdi, & !< derivative of grid mapping at cell face (m/m)
- dyfdj,dzfdl
+ xc,yc,zc, & !< center point of a cell (m)
+ xf,yf,zf, & !< face of a cell in pos. direction from xc (m)
+ dxcdi,dycdj,dzcdl, & !< derivative of grid mapping at cell center (m/m)
+ dxfdi,dyfdj,dzfdl !< derivative of grid mapping at cell face (m/m)
real(kind=pr) :: &
- xmin, & !< lower local boundary in global coordinates (m)
- xmax, & !< upper local boundary in global coordinates (m)
- ymin,ymax,zmin,zmax, &
- deltax, & !< shift between procs (m)
- volTot !< total domain volume (m**3)
+ xmin,ymin,zmin, & !< lower local boundary in global coordinates (m)
+ xmax,ymax,zmax, & !< upper local boundary in global coordinates (m)
+ deltax, & !< shift between processors (m)
+ volTot !< total domain volume (m**3)
integer, allocatable, dimension(:,:,:) :: &
- celltype !< celltype
+ celltype !< celltype
+ integer, parameter :: active=1, passive=2, wall=3
- integer, parameter :: active=1, ghostcell=2, wall=3
integer :: &
- gp, & !< local number of grid cells
- dimgp, & !< local number of cells containing fluid
- dimgptot, & !< glocal number of cells containing fluid
- ii, & !< max index local grid
- jj,ll, &
- gc, & !< ghost cells between processors
- wc, & !< wall boundary cells
- imin, & !< local min index of cell containing fluid
- imax, & !< local max index of cell containing fluid
- jmin,jmax,lmin,lmax
+ gp, & !< local number of grid cells
+ dimgp, & !< local number of cells containing fluid
+ dimgptot, & !< global number of cells containing fluid
+ ii,jj,ll, & !< max index local grid
+ imin,jmin,lmin, & !< local min index of cell containing fluid
+ imax,jmax,lmax !< local max index of cell containing fluid
! fluid
- real(kind=pr) :: muf !< fluid dynamic viscosity (kg/s/m)
+ real(kind=pr) :: muf !< fluid dynamic viscosity (kg/s/m)
real(kind=pr), allocatable, dimension(:,:,:) :: &
- u, & !< fluid velocity (m/s)
- v,w, &
- p, & !< fluid pressure (N/m**2)
- dudt, & !< term for time integration (m/s)
- dvdt,dwdt, &
- u01, & !< fluid velocity previous time-step
- v01,w01, &
- fsx, & !< source term particles -> fluid (m/s**2)
- fsy,fsz, &
- defect_c, & !< defect correction term
+ u,v,w, & !< fluid velocity (m/s)
+ p, & !< fluid pressure (N/m**2)
+ dudt,dvdt,dwdt, & !< term for time integration (m/s)
+ u01,v01,w01, & !< fluid velocity previous time-step
+ fsx,fsy,fsz, & !< source term particles -> fluid (m/s**2)
+ defect_c, & !< defect correction terms
defect_u,defect_v,defect_w
! particles
integer :: &
- np, & !< local number of particles
- npp, & !< local number of particles before some operation
- maxnp !< maximum possible local number of particles
+ np, & !< local number of particles
+ npTot, & !< total number of particles
+ npp, & !< local number of particles before some operation
+ maxnp, & !< maximum possible local number of particles
+ nseedLoc !< local counter for seeded particles
real(kind=pr), allocatable, dimension(:) :: &
- up, & !< particle velocity (m/s)
- vp,wp, &
- uf, & !< average fluid velocity around a particle (m/s)
- vf,wf, &
- uf01, & !< fluid vel. around particle prev. time-step (m/s)
- vf01,wf01, &
- xp, & !< particle position (m)
- yp,zp, &
- radp, & !< particle radius (m)
- q_el, & !< particle charge (C)
- fx_el, & !< electrostatic force on particle (m/s**2)
- fy_el,fz_el, &
- partn !< number of particles represented by parcel
+ up,vp,wp, & !< particle velocity (m/s)
+ uf,vf,wf, & !< average fluid velocity around a particle (m/s)
+ uf01,vf01,wf01, & !< fluid vel. around particle prev. time-step (m/s)
+ xp,yp,zp, & !< particle position (m)
+ radp, & !< particle radius (m)
+ q_el, & !< particle charge (C)
+ fx_el,fy_el,fz_el, & !< electrostatic force on particle (m/s**2)
+ partn !< number of particles represented by parcel
integer, allocatable, dimension(:) :: &
- wcollnum, & !< number of particle-wall collisions
- ppcollnum !< number of particle-particle collisions
+ wcollnum, & !< number of particle-wall collisions
+ ppcollnum, & !< number of particle-particle collisions
+ nGlob !< global particle id
+ integer, allocatable, dimension(:,:,:) :: &
+ npic !< number of particles in a Eulerian cell
+ integer, allocatable, dimension(:,:,:,:) :: &
+ nic !< indices of particles in a Eulerian cell
+
! parallel
- real(kind=pr) :: timebeg,timecom,timenow,timeend
+ real(kind=pr) :: timebeg,timecom(10),timenow,timeend
integer :: &
mpi_pr, &
- myid, & !< id of current processor
- nrprocs, & !< total number of processors
+ myid, & !< id of current processor
+ nrprocs, & !< total number of processors
mpierr, &
- next, & !< id of next processor
- prev !< id of previous processor
+ next,prev !< id of next and previous processor
integer, allocatable, dimension(:) :: &
mpistatus
! electrostatics
real(kind=pr), allocatable, dimension(:,:,:) :: &
- rho_el, & !< electric charge density (C/m**3)
- phi_el, & !< electric potential (V)
- Ex_el, & !< electric field strength (V/m)
- Ey_el,Ez_el
- real(kind=pr) &
- Ew, & !< duct Young's modulus (kg/s**2/m)
- Ep, & !< particle Young's modulus (kg/s**2/m)
- h0, & !< particle-wall separation distance (m)
- r_ep, & !< particle elec. surface resistivity
- q0, & !< particle initial charge (C)
- Vc, & !< contact potential difference (V)
- nyp, & !< particle Poisson ratio
- nyw, & !< duct Poisson ratio
- restRatio, &!< particle material resitution ratio
- epsw !< relative permittivity system
+ rho_el, & !< electric charge density (C/m**3)
+ phi_el, & !< electric potential (V)
+ Ex_el,Ey_el,Ez_el !< electric field strength (V/m)
end module var
diff --git a/src/writedat_fluid_xy.f90 b/src/writedat_fluid_xy.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a3c0e89fbc509af679ccb1756eec5beac996be5c
--- /dev/null
+++ b/src/writedat_fluid_xy.f90
@@ -0,0 +1,93 @@
+!####################################################################
+!> @author Holger Grosshans
+ subroutine writedat_ufluid_xy
+ use var
+ character(70) :: filename,rowfmt
+ integer :: i,j,l
+
+100 format(es11.4e2)
+ write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+ l=int(ll/2)
+
+ write(filename,'(a,i3.3,a,i6.6,a)') &
+ 'results/fluid_u_xy_p',myid,'_',nt,'.dat'
+ open(10,file=filename)
+
+ write(10,'(a,es14.6e2,a,x,a)') '# u: t=',t,' s --',version
+
+ write(10,'(a11)',advance='no') '# y/x'
+ write(10,fmt=rowfmt) (xf(i),i=imin,imax+1)
+ do j=jmin-1,jmax+1
+ write(10,fmt=100,advance='no') yc(j)
+ write(10,fmt=rowfmt) (u(i,j,l),i=imin,imax+1)
+ enddo
+
+ close(10)
+
+ return
+
+ end
+
+!####################################################################
+!> @author Holger Grosshans
+ subroutine writedat_vfluid_xy
+ use var
+ character(70) :: filename,rowfmt
+ integer :: i,j,l
+
+100 format(es11.4e2)
+ write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+ l=int(ll/2)
+
+ write(filename,'(a,i3.3,a,i6.6,a)') &
+ 'results/fluid_v_xy_p',myid,'_',nt,'.dat'
+ open(10,file=filename)
+
+ write(10,'(a,es14.6e2,a,x,a)') '# v: t=',t,' s --',version
+
+ write(10,'(a11)',advance='no') '# y/x'
+ write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+ do j=jmin-1,jmax+1
+ write(10,fmt=100,advance='no') yf(j)
+ write(10,fmt=rowfmt) (v(i,j,l),i=imin,imax+1)
+ enddo
+
+ close(10)
+
+ return
+
+ end
+
+!####################################################################
+!> @author Holger Grosshans
+ subroutine writedat_wfluid_xy
+ use var
+ character(70) :: filename,rowfmt
+ integer :: i,j,l
+
+100 format(es11.4e2)
+ write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+
+ l=int(ll/2)
+
+ write(filename,'(a,i3.3,a,i6.6,a)') &
+ 'results/fluid_w_xy_p',myid,'_',nt,'.dat'
+ open(10,file=filename)
+
+ write(10,'(a,es14.6e2,a,x,a)') '# w: t=',t,' s --',version
+
+ write(10,'(a11)',advance='no') '# y/x'
+ write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+ do j=jmin-1,jmax+1
+ write(10,fmt=100,advance='no') yc(j)
+ write(10,fmt=rowfmt) (w(i,j,l),i=imin,imax+1)
+ enddo
+
+ close(10)
+
+ return
+
+ end
+
diff --git a/src/writedat_fluid.f90 b/src/writedat_fluid_xz.f90
similarity index 86%
rename from src/writedat_fluid.f90
rename to src/writedat_fluid_xz.f90
index 1f0f9d6b2a8141394d1c6f0837ceecc475d80bef..77989532278fb60fbf9f696ca40e8f5d23804038 100644
--- a/src/writedat_fluid.f90
+++ b/src/writedat_fluid_xz.f90
@@ -2,7 +2,7 @@
!> @author Holger Grosshans
subroutine writedat_ufluid_xz
use var
- character*70 filename,rowfmt
+ character(70) :: filename,rowfmt
integer :: i,j,l
100 format(es11.4e2)
@@ -14,7 +14,7 @@
'results/fluid_u_xz_p',myid,'_',nt,'.dat'
open(10,file=filename)
- write(10,'(a,es11.2e2,a,x,a)') '# u: t=',t,' s --',version
+ write(10,'(a,es14.6e2,a,x,a)') '# u: t=',t,' s --',version
write(10,'(a11)',advance='no') '# z/x'
write(10,fmt=rowfmt) (xf(i),i=imin,imax+1)
@@ -33,7 +33,7 @@
!> @author Holger Grosshans
subroutine writedat_vfluid_xz
use var
- character*70 filename,rowfmt
+ character(70) :: filename,rowfmt
integer :: i,j,l
100 format(es11.4e2)
@@ -45,7 +45,7 @@
'results/fluid_v_xz_p',myid,'_',nt,'.dat'
open(10,file=filename)
- write(10,'(a,es11.2e2,a,x,a)') '# v: t=',t,' s --',version
+ write(10,'(a,es14.6e2,a,x,a)') '# v: t=',t,' s --',version
write(10,'(a11)',advance='no') '# z/x'
write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
@@ -64,7 +64,7 @@
!> @author Holger Grosshans
subroutine writedat_wfluid_xz
use var
- character*70 filename,rowfmt
+ character(70) :: filename,rowfmt
integer :: i,j,l
100 format(es11.4e2)
@@ -76,7 +76,7 @@
'results/fluid_w_xz_p',myid,'_',nt,'.dat'
open(10,file=filename)
- write(10,'(a,es11.2e2,a,x,a)') '# w: t=',t,' s --',version
+ write(10,'(a,es14.6e2,a,x,a)') '# w: t=',t,' s --',version
write(10,'(a11)',advance='no') '# z/x'
write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
diff --git a/src/writedat_particles.f90 b/src/writedat_particles.f90
index 3cc7546adb6fe9c3d6257dc30f5ddbbbb4a0e247..187b225928fc3a0a85eb0e75c05445693723b5dd 100644
--- a/src/writedat_particles.f90
+++ b/src/writedat_particles.f90
@@ -5,7 +5,7 @@
use var
real(kind=pr) :: partmass
integer n
- character*80 filename
+ character(80) :: filename
100 format(10(es13.4e2,x),2(i6,x))
@@ -14,6 +14,8 @@
open(10,file=filename)
+ write(10,'(a,es14.6e2,a,x,a)') '# t=',t,' s --',version
+
do n=1,np
partmass=4._pr/3._pr*pi*rhop*partn(n)*radp(n)**3
write(10,100) xp(n),yp(n),zp(n),radp(n),up(n),vp(n),wp(n), &
diff --git a/src/writevtk_fluid_xy.f90 b/src/writevtk_fluid_xy.f90
new file mode 100644
index 0000000000000000000000000000000000000000..3ef111e191837c54e130d8b3eddff940733cf2f7
--- /dev/null
+++ b/src/writevtk_fluid_xy.f90
@@ -0,0 +1,144 @@
+!> @author Holger Grosshans
+ subroutine writevtk_fluid_xy
+ use var
+ integer :: i,j,l,m
+ character(70) :: filename,filename2,rowfmt,vecfmt
+ logical :: file_ex
+
+
+100 format(es11.4e2)
+ write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+ write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))'
+
+ l=int(ll/2.)+1
+
+ write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xy_p',myid,'_',nt,'.vtk'
+ open(10,file=filename)
+
+! write visit container file
+ filename2='results/fluid_xy.visit'
+ if(myid.eq.0) then
+ inquire(file=filename2,exist=file_ex)
+ if (file_ex.and.nt.gt.1) then
+ open(11,file=filename2,access='append')
+ else
+ open(11,file=filename2)
+ write(11,'(a8,x,i3)') '!NBLOCKS',nrprocs
+ endif
+ do m=1,nrprocs
+ write(11,'(a,i3.3,a,i6.6,a)') 'fluid_xy_p',(m-1),'_',nt,'.vtk'
+ enddo
+ close(11)
+ endif
+
+
+ write(10,'(a)') '# vtk DataFile Version 3.0'
+ write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
+ write(10,'(a)') 'ASCII'
+ write(10,'(a)') 'DATASET RECTILINEAR_GRID'
+ write(10,'(a10,i4,i4,a4)') &
+ 'DIMENSIONS',(imax-imin+2),(jmax-jmin+3),'1'
+
+ write(10,*) 'X_COORDINATES',(imax-imin+2),'float'
+ write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+
+ write(10,*) 'Y_COORDINATES',(jmax-jmin+3),'float'
+ write(10,fmt=rowfmt) (yc(j),j=jmin-1,jmax+1)
+
+ write(10,*) 'Z_COORDINATES 1 float'
+ write(10,100) zc(l)
+
+ write(10,*) 'POINT_DATA ',(imax-imin+2)*(jmax-jmin+3)
+
+ write(10,'(a)') 'SCALARS u float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)),i=imin,imax+1)
+ enddo
+
+ write(10,'(a)') 'SCALARS v float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)),i=imin,imax+1)
+ enddo
+
+ write(10,'(a)') 'SCALARS w float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
+ enddo
+
+! write(10,'(a)') 'SCALARS uraw float 1'
+! write(10,'(a)') 'LOOKUP_TABLE default '
+! do j=jmin-1,jmax+1
+! write(10,fmt=rowfmt) (u(i,j,l),i=imin,imax+1)
+! enddo
+!
+! write(10,'(a)') 'SCALARS vraw float 1'
+! write(10,'(a)') 'LOOKUP_TABLE default '
+! do j=jmin-1,jmax+1
+! write(10,fmt=rowfmt) (v(i,j,l),i=imin,imax+1)
+! enddo
+!
+! write(10,'(a)') 'SCALARS wraw float 1'
+! write(10,'(a)') 'LOOKUP_TABLE default '
+! do j=jmin-1,jmax+1
+! write(10,fmt=rowfmt) (w(i,j,l),i=imin,imax+1)
+! enddo
+
+ write(10,'(a)') 'SCALARS p float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (p(i,j,l),i=imin,imax+1)
+ enddo
+
+ write(10,'(a)') 'SCALARS rho_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (rho_el(i,j,l),i=imin,imax+1)
+ enddo
+
+ write(10,'(a)') 'SCALARS phi_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (phi_el(i,j,l),i=imin,imax+1)
+ enddo
+
+ write(10,'(a)') 'SCALARS Ex_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (Ex_el(i,j,l),i=imin,imax+1)
+ enddo
+
+ write(10,'(a)') 'SCALARS Ey_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (Ey_el(i,j,l),i=imin,imax+1)
+ enddo
+
+ write(10,'(a)') 'SCALARS Ez_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1
+ write(10,fmt=rowfmt) (Ez_el(i,j,l),i=imin,imax+1)
+ enddo
+
+! write(10,'(a)') 'SCALARS Emag_el float 1'
+! write(10,'(a)') 'LOOKUP_TABLE default '
+! do j=jmin-1,jmax+1
+! write(10,fmt=rowfmt) &
+! (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),i=imin,imax+1)
+! enddo
+!
+!
+! write(10,'(a)') 'VECTORS uvw float'
+! do j=jmin-1,jmax+1
+! write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), &
+! 0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), &
+! 0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
+! enddo
+
+
+ close(10)
+
+ return
+ end
diff --git a/src/writevtk_fluid_xyz.f90 b/src/writevtk_fluid_xyz.f90
new file mode 100644
index 0000000000000000000000000000000000000000..22a942959b9325705c337bed4589024a3706244a
--- /dev/null
+++ b/src/writevtk_fluid_xyz.f90
@@ -0,0 +1,138 @@
+!> @author Holger Grosshans
+ subroutine writevtk_fluid_xyz
+ use var
+ real(kind=pr) :: pp
+ integer :: i,j,l,m
+ character(70) :: filename,filename2,rowfmt,vecfmt
+ logical :: file_ex
+
+
+100 format(es11.4e2)
+ write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+ write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))'
+
+ pp=0._pr
+
+ write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xyz_p',myid,'_',nt,'.vtk'
+ open(10,file=filename)
+
+! write visit container file
+ filename2='results/fluid_xyz.visit'
+ if(myid.eq.0) then
+ inquire(file=filename2,exist=file_ex)
+ if (file_ex.and.nt.gt.1) then
+ open(11,file=filename2,access='append')
+ else
+ open(11,file=filename2)
+ write(11,'(a8,x,i3)') '!NBLOCKS',nrprocs
+ endif
+ do m=1,nrprocs
+ write(11,'(a,i3.3,a,i6.6,a)') 'fluid_xyz_p',(m-1),'_',nt,'.vtk'
+ enddo
+ close(11)
+ endif
+
+
+ write(10,'(a)') '# vtk DataFile Version 3.0'
+ write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
+ write(10,'(a)') 'ASCII'
+ write(10,'(a)') 'DATASET RECTILINEAR_GRID'
+ write(10,'(a10,3(i4))') &
+ 'DIMENSIONS',(imax-imin+2),(jmax-jmin+3),(lmax-lmin+3)
+
+ write(10,*) 'X_COORDINATES',(imax-imin+2),'float'
+ write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
+
+ write(10,*) 'Y_COORDINATES',(jmax-jmin+3),'float'
+ write(10,fmt=rowfmt) (yc(j),j=jmin-1,jmax+1)
+
+ write(10,*) 'Z_COORDINATES',(lmax-lmin+3),'float'
+ write(10,fmt=rowfmt) (zc(l),l=lmin-1,lmax+1)
+
+ write(10,*) 'POINT_DATA ',(imax-imin+2)*(jmax-jmin+3)*(lmax-lmin+3)
+
+ write(10,'(a)') 'SCALARS u float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS v float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS w float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS p float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (p(i,j,l),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS rho_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (rho_el(i,j,l),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS phi_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (phi_el(i,j,l),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS Ex_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (Ex_el(i,j,l),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS Ey_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (Ey_el(i,j,l),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS Ez_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) (Ez_el(i,j,l),i=imin,imax+1)
+ enddo; enddo
+
+ write(10,'(a)') 'SCALARS Emag_el float 1'
+ write(10,'(a)') 'LOOKUP_TABLE default '
+ do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+ write(10,fmt=rowfmt) &
+ (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),i=imin,imax+1)
+ enddo; enddo
+
+! write(10,'(a)') 'VECTORS uvw float'
+! do j=jmin-1,jmax+1; do l=lmin-1,lmax+1
+! write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), &
+! 0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), &
+! 0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
+! enddo; enddo
+
+
+! write(10,'(a)') 'VECTORS coordinates float 1'
+! write(10,'(a)') 'LOOKUP_TABLE default '
+! do l=lmin-1,lmax+1
+! do i=imin,imax+1
+! write(10,'(3(i8))') i,j,l
+! enddo
+! enddo
+
+
+ close(10)
+
+
+
+ return
+
+ end
diff --git a/src/writevtk_fluid_xz.f90 b/src/writevtk_fluid_xz.f90
index b626ab4a2926912bab0a506f5fd677b973eeda97..832508e46c8367f4b3fb8f6beeafa5f079ae6d4f 100644
--- a/src/writevtk_fluid_xz.f90
+++ b/src/writevtk_fluid_xz.f90
@@ -1,28 +1,27 @@
!> @author Holger Grosshans
subroutine writevtk_fluid_xz
use var
- real(kind=pr), dimension(ii) :: temp4
- real(kind=pr) :: temp,temp1,temp2,temp3,pp
+ real(kind=pr) :: pp
integer :: i,j,l,m
- character*70 :: filename,filename2,rowfmt
+ character(70) :: filename,filename2,rowfmt,vecfmt
logical :: file_ex
100 format(es11.4e2)
write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
+ write(vecfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,3(es11.4e2)))'
j=int(jj/2.)+1
pp=0._pr
- write(filename,'(a,i3.3,a,i6.6,a)') &
- 'results/fluid_xz_p',myid,'_',nt,'.vtk'
+ write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_xz_p',myid,'_',nt,'.vtk'
open(10,file=filename)
! write visit container file
filename2='results/fluid_xz.visit'
if(myid.eq.0) then
inquire(file=filename2,exist=file_ex)
- if (file_ex.and.nt.ne.0) then
+ if (file_ex.and.nt.gt.1) then
open(11,file=filename2,access='append')
else
open(11,file=filename2)
@@ -35,12 +34,12 @@
endif
- write(10,'(a)') '# vtk DataFile Version 2.0'
- write(10,'(a,es11.2e2,a,x,a)') 'fluid data: t=',t,' s --',version
+ write(10,'(a)') '# vtk DataFile Version 3.0'
+ write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
write(10,'(a)') 'ASCII'
write(10,'(a)') 'DATASET RECTILINEAR_GRID'
- write(10,'(a10,i8,a4)') &
- 'DIMENSIONS',(imax-imin+2)*(lmax-lmin+3),'1 1'
+ write(10,'(a10,i4,a4,i4)') &
+ 'DIMENSIONS',(imax-imin+2),'1',(lmax-lmin+3)
write(10,*) 'X_COORDINATES',(imax-imin+2),'float'
write(10,fmt=rowfmt) (xc(i),i=imin,imax+1)
@@ -125,95 +124,23 @@
write(10,fmt=rowfmt) (Ez_el(i,j,l),i=imin,imax+1)
enddo
- write(10,'(a)') 'SCALARS Emag_el float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- do l=lmin-1,lmax+1
- write(10,fmt=rowfmt) &
- (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),i=imin,imax+1)
- enddo
-
-
-!
-! write(10,'(a)') 'scalars emag_el float 1'
-! write(10,'(a)') 'lookup_table default '
-! do l=1,ll
-! do i=1,ii
-! io=lst(l)+ist(i)+j
-! temp=(ex_el(i,j,l)**2+ey_el(i,j,l)**2+ez_el(i,j,l)**2)**(0.5)
-! if(abs(temp).lt.1.0d-99) then
-! write(10,100) 0.0
-! else
-! write(10,100) temp
-! endif
-! enddo
-! enddo
-!
-!
-!
-! write(10,'(a)') 'vectors E_el float 1'
-! write(10,'(a)') 'lookup_table default '
-! do l=1,ll
-! do i=1,ii
-! io=lst(l)+ist(i)+j
-! temp1=ex_el(i,j,l)
-! if(abs(temp1).lt.1.0d-99) temp1=0.0
-! temp2=ey_el(i,j,l)
-! if(abs(temp1).lt.1.0d-99) temp2=0.0
-! temp3=ez_el(i,j,l)
-! if(abs(temp1).lt.1.0d-99) temp3=0.0
-! write(10,'(3(es10.2e1))') temp1,temp2,temp3
-! enddo
-! enddo
-
-! write(10,'(a)') 'VECTORS coordinates float 1'
+! write(10,'(a)') 'SCALARS Emag_el float 1'
! write(10,'(a)') 'LOOKUP_TABLE default '
! do l=lmin-1,lmax+1
-! do i=imin,imax+1
-! write(10,'(3(i8))') i,j,l
-! enddo
+! write(10,fmt=rowfmt) &
+! (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),i=imin,imax+1)
! enddo
-
-! write(10,'(a)') 'vectors uvw float 1'
-! write(10,'(a)') 'lookup_table default '
-! do l=1,ll
-! do i=1,ii
-! io=lst(l)+ist(i)+j
-!
-! if (i.eq.ii) then
-! ie=io
-! else
-! ie=lst(l)+ist(i+1)+j
-! endif
-! if (i.eq.1) then
-! iw=io
-! else
-! iw=lst(l)+ist(i-1)+j
-! endif
-!
-! is=lst(l)+ist(i)+j-1
-! no=lst(l)+ist(i)+j+1
!
-! if (l.eq.ll) then
-! nf=io
-! else
-! nf=lst(l+1)+ist(i)+j
-! endif
-! if (l.eq.1) then
-! nb=io
-! else
-! nb=lst(l-1)+ist(i)+j
-! endif
!
-! temp1=u(i,j,l)/(rho(i,j,l)+rho(i+1,j,l))+u(i-1,j,l)/(rho(i,j,l)+rho(i-1,j,l))
-! if(abs(temp1).lt.1.0e-9) temp1=0.0
-! temp2=v(i,j,l)/(rho(i,j,l)+rho(i,j+1,l))+v(i,j-1,l)/(rho(i,j,l)+rho(i,j-1,l))
-! if(abs(temp2).lt.1.0e-9) temp2=0.0
-! temp3=w(i,j,l)/(rho(i,j,l)+rho(i,j,l+1))+w(i,j,l-1)/(rho(i,j,l)+rho(i,j,l-1))
-! if(abs(temp3).lt.1.0e-9) temp3=0.0
-! write(10,'(3(es10.2e1))') temp1,temp2,temp3
-! enddo
+! write(10,'(a)') 'VECTORS uvw float'
+! do l=lmin-1,lmax+1
+! write(10,fmt='(3(es12.2e2))') (0.5_pr*(u(i,j,l)+u(max(i-1,1),j,l)), &
+! 0.5_pr*(v(i,j,l)+v(i,max(j-1,1),l)), &
+! 0.5_pr*(w(i,j,l)+w(i,j,max(l-1,1))),i=imin,imax+1)
! enddo
-
+
+
+
close(10)
diff --git a/src/writevtk_fluid_yz.f90 b/src/writevtk_fluid_yz.f90
index 9c537f38ff5f81d6e70f93f2114385ae767d0107..59f40d580f4a128435ca9cf0142d611211c78885 100644
--- a/src/writevtk_fluid_yz.f90
+++ b/src/writevtk_fluid_yz.f90
@@ -2,9 +2,9 @@
subroutine writevtk_fluid_yz
use var
real(kind=pr), dimension(jj) :: temp4
- real(kind=pr) :: temp,temp1,temp2,temp3,pp
+ real(kind=pr) :: temp,temp1,temp2,temp3
integer :: i,j,l,m
- character*70 :: filename,filename2,rowfmt
+ character(70) :: filename,filename2,rowfmt
logical :: file_ex
@@ -12,17 +12,15 @@
write(rowfmt,'(a,i4,a)') '(',(jmax-jmin+3),'(1x,es11.4e2))'
i=int(ii/2.)+1
- pp=0._pr
- write(filename,'(a,i3.3,a,i6.6,a)') &
- 'results/fluid_yz_p',myid,'_',nt,'.vtk'
+ write(filename,'(a,i3.3,a,i6.6,a)') 'results/fluid_yz_p',myid,'_',nt,'.vtk'
open(10,file=filename)
! write visit container file
filename2='results/fluid_yz.visit'
if(myid.eq.0) then
inquire(file=filename2,exist=file_ex)
- if (file_ex.and.nt.ne.1) then
+ if (file_ex.and.nt.gt.1) then
open(11,file=filename2,access='append')
else
open(11,file=filename2)
@@ -35,12 +33,12 @@
endif
- write(10,'(a)') '# vtk DataFile Version 2.0'
- write(10,'(a,es11.2e2,a,x,a)') 'fluid data: t=',t,' s --',version
+ write(10,'(a)') '# vtk DataFile Version 3.0'
+ write(10,'(a,es14.6e2,a,x,a)') 'fluid data: t=',t,' s --',version
write(10,'(a)') 'ASCII'
write(10,'(a)') 'DATASET RECTILINEAR_GRID'
- write(10,'(a10,i8,a4)') &
- 'DIMENSIONS',(jmax-jmin+3)*(lmax-lmin+3),'1 1'
+ write(10,'(a12,i8,i8)') &
+ 'DIMENSIONS 1',(jmax-jmin+3),(lmax-lmin+3)
write(10,*) 'X_COORDINATES 1 float'
write(10,100) xc(i)
@@ -125,99 +123,16 @@
write(10,fmt=rowfmt) (Ez_el(i,j,l),j=jmin-1,jmax+1)
enddo
- write(10,'(a)') 'SCALARS Emag_el float 1'
- write(10,'(a)') 'LOOKUP_TABLE default '
- do l=lmin-1,lmax+1
- write(10,fmt=rowfmt) &
- (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),j=jmin-1,jmax+1)
- enddo
-
-
-!
-! write(10,'(a)') 'scalars emag_el float 1'
-! write(10,'(a)') 'lookup_table default '
-! do l=1,ll
-! do i=1,ii
-! io=lst(l)+ist(i)+j
-! temp=(ex_el(i,j,l)**2+ey_el(i,j,l)**2+ez_el(i,j,l)**2)**(0.5)
-! if(abs(temp).lt.1.0d-99) then
-! write(10,100) 0.0
-! else
-! write(10,100) temp
-! endif
-! enddo
-! enddo
-!
-!
-!
-! write(10,'(a)') 'vectors E_el float 1'
-! write(10,'(a)') 'lookup_table default '
-! do l=1,ll
-! do i=1,ii
-! io=lst(l)+ist(i)+j
-! temp1=ex_el(i,j,l)
-! if(abs(temp1).lt.1.0d-99) temp1=0.0
-! temp2=ey_el(i,j,l)
-! if(abs(temp1).lt.1.0d-99) temp2=0.0
-! temp3=ez_el(i,j,l)
-! if(abs(temp1).lt.1.0d-99) temp3=0.0
-! write(10,'(3(es10.2e1))') temp1,temp2,temp3
-! enddo
-! enddo
-
-! write(10,'(a)') 'VECTORS coordinates float 1'
+! write(10,'(a)') 'SCALARS Emag_el float 1'
! write(10,'(a)') 'LOOKUP_TABLE default '
! do l=lmin-1,lmax+1
-! do i=imin,imax+1
-! write(10,'(3(i8))') i,j,l
-! enddo
+! write(10,fmt=rowfmt) &
+! (sqrt(Ex_el(i,j,l)**2+Ey_el(i,j,l)**2+Ez_el(i,j,l)**2)**(0.5),j=jmin-1,jmax+1)
! enddo
-! write(10,'(a)') 'vectors uvw float 1'
-! write(10,'(a)') 'lookup_table default '
-! do l=1,ll
-! do i=1,ii
-! io=lst(l)+ist(i)+j
-!
-! if (i.eq.ii) then
-! ie=io
-! else
-! ie=lst(l)+ist(i+1)+j
-! endif
-! if (i.eq.1) then
-! iw=io
-! else
-! iw=lst(l)+ist(i-1)+j
-! endif
-!
-! is=lst(l)+ist(i)+j-1
-! no=lst(l)+ist(i)+j+1
-!
-! if (l.eq.ll) then
-! nf=io
-! else
-! nf=lst(l+1)+ist(i)+j
-! endif
-! if (l.eq.1) then
-! nb=io
-! else
-! nb=lst(l-1)+ist(i)+j
-! endif
-!
-! temp1=u(i,j,l)/(rho(i,j,l)+rho(i+1,j,l))+u(i-1,j,l)/(rho(i,j,l)+rho(i-1,j,l))
-! if(abs(temp1).lt.1.0e-9) temp1=0.0
-! temp2=v(i,j,l)/(rho(i,j,l)+rho(i,j+1,l))+v(i,j-1,l)/(rho(i,j,l)+rho(i,j-1,l))
-! if(abs(temp2).lt.1.0e-9) temp2=0.0
-! temp3=w(i,j,l)/(rho(i,j,l)+rho(i,j,l+1))+w(i,j,l-1)/(rho(i,j,l)+rho(i,j,l-1))
-! if(abs(temp3).lt.1.0e-9) temp3=0.0
-! write(10,'(3(es10.2e1))') temp1,temp2,temp3
-! enddo
-! enddo
close(10)
-
return
-
end
diff --git a/src/writevtk_grid.f90 b/src/writevtk_grid.f90
index b822d337cb63e314850758859afe57160f5bf602..346d0e66ad9df4ee1b9ff654a0647eb175825946 100644
--- a/src/writevtk_grid.f90
+++ b/src/writevtk_grid.f90
@@ -2,7 +2,7 @@
subroutine writevtk_grid
use var
integer :: i,j,l
- character*70 filename,rowfmt
+ character(70) :: filename,rowfmt
100 format(es11.4e2)
write(rowfmt,'(a,i4,a)') '(',(imax-imin+2),'(1x,es11.4e2))'
diff --git a/src/writevtk_particles.f90 b/src/writevtk_particles.f90
index 1e74710fc9ce791737a11daa0b2e0544b8320b7c..b519f0691331ef0f9ec5b81bbf77a85b8febc4b9 100644
--- a/src/writevtk_particles.f90
+++ b/src/writevtk_particles.f90
@@ -4,7 +4,7 @@
subroutine writevtk_particles
use var
integer :: i,j,l,n,m
- character*80 :: filename,filename2,rowfmt,rowfm2
+ character(80) :: filename,filename2,rowfmt,rowfm2
logical :: file_ex
@@ -36,8 +36,8 @@
endif
- write(10,'(a)') '# vtk DataFile Version 2.0'
- write(10,'(a,x,a)') 'particle data --',version
+ write(10,'(a)') '# vtk DataFile Version 3.0'
+ write(10,'(a,es14.6e2,a,x,a)') 'particle data: t=',t,' s --',version
write(10,'(a)') 'ASCII'
write(10,'(a)') 'DATASET UNSTRUCTURED_GRID'
write(10,*)
@@ -109,7 +109,7 @@
write(10,'(a)') 'SCALARS np FLOAT 1'
write(10,'(a)') 'LOOKUP_TABLE default '
- write(10,fmt=rowfm2) (n,n=1,np)
+ write(10,fmt=rowfm2) (nGlob(n),n=1,np)
write(10,'(a)')
write(10,'(a)') 'VECTORS UPVPWP FLOAT'