From b3ed954b28b1bcc3937ed7bb41128fc4bb4b3d08 Mon Sep 17 00:00:00 2001 From: Holger <holger.grosshans@ptb.de> Date: Mon, 20 Apr 2020 10:39:50 +0200 Subject: [PATCH] modified: src/electrostatics.f90 modified: src/makefile --- src/electrostatics.f90 | 93 ++++++++++++++++++++++++++---------------- src/makefile | 2 +- 2 files changed, 59 insertions(+), 36 deletions(-) diff --git a/src/electrostatics.f90 b/src/electrostatics.f90 index 4d4766f..746ebd5 100644 --- a/src/electrostatics.f90 +++ b/src/electrostatics.f90 @@ -215,51 +215,74 @@ 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,n1end - 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 3 i=1,ii + do 3 j=1,jj + do 3 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 +3 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 diff --git a/src/makefile b/src/makefile index 01e3e2a..148b351 100644 --- a/src/makefile +++ b/src/makefile @@ -15,7 +15,7 @@ INC = # gcc: CMP = mpifort #O3: optimization -FLAGS = -O3 -mcmodel=medium +FLAGS = -O3 -mcmodel=medium -pg #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 FFLAGS = -c $(FLAGS) -- GitLab