JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Reorganization of solver

Browse source 
I started grouping similar modules in subfolders to ease the expansion
process.
This commit is contained in:
Jorge Gonzalez 2022-12-24 12:59:23 +01:00
commit d9a1869564
13 changed files with 341 additions and 297 deletions
Download patch file Download diff file Expand all files Collapse all files

21
src/modules/makefile
View file

@ -1,19 +1,24 @@
OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \
moduleOutput.o moduleInput.o moduleSolver.o \
output.o moduleInput.o solver.o \
moduleCollisions.o moduleTable.o moduleParallel.o \
moduleEM.o moduleRandom.o moduleMath.o \
moduleProbe.o
moduleRandom.o moduleMath.o moduleProbe.o
all: $(OBJS)
mesh.o: moduleCollisions.o moduleBoundary.o moduleAverage.o
output.o: moduleMath.o moduleRefParam.o
$(MAKE) -C output all
mesh.o: moduleCollisions.o moduleBoundary.o output.o
$(MAKE) -C mesh all
solver.o: moduleSpecies.o moduleProbe.o moduleRandom.o output.o mesh.o
$(MAKE) -C solver all
moduleCollisions.o: moduleList.o moduleMath.o moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleInput.o: moduleParallel.o moduleSolver.o moduleInject.o moduleInput.f90
moduleInput.o: moduleParallel.o solver.o moduleInject.o moduleInput.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleInject.o: moduleInject.f90
@ -22,12 +27,6 @@ moduleInject.o: moduleInject.f90
moduleList.o: moduleConstParam.o moduleErrors.o moduleCaseParam.o moduleSpecies.o moduleList.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleOutput.o: moduleMath.o moduleRefParam.o moduleOutput.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleSolver.o: moduleProbe.o moduleEM.o moduleAverage.o moduleSolver.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleProbe.o: mesh.o moduleProbe.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

8
src/modules/mesh/1DRad/moduleMesh1DRad.f90
View file

@ -265,11 +265,11 @@ MODULE moduleMesh1DRad
l = 2.D0*detJ
self%volume = r*l
!Computes volume per node
xi = (/-5.D-1, 0.D0, 0.D0/)
r = DOT_PRODUCT(self%fPsi(xi),self%r)
Xii = (/-5.D-1, 0.D0, 0.D0/)
r = DOT_PRODUCT(self%fPsi(Xii),self%r)
self%arNodes(1) = fPsi(1)*r*l
xi = (/ 5.D-1, 0.D0, 0.D0/)
r = DOT_PRODUCT(self%fPsi(xi),self%r)
Xii = (/ 5.D-1, 0.D0, 0.D0/)
r = DOT_PRODUCT(self%fPsi(Xii),self%r)
self%arNodes(2) = fPsi(2)*r*l
END SUBROUTINE areaRad

3
src/modules/mesh/makefile
View file

@ -26,3 +26,6 @@ moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90
inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
$(MAKE) -C inout all
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

7
src/modules/output/makefile Normal file
View file

@ -0,0 +1,7 @@
all: moduleAverage.o moduleOutput.o
moduleAverage.o: moduleOutput.o moduleAverage.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

0
src/modules/moduleAverage.f90 → src/modules/output/moduleAverage.f90
View file

0
src/modules/moduleOutput.f90 → src/modules/output/moduleOutput.f90
View file

8
src/modules/solver/electromagnetic/makefile Normal file
View file

@ -0,0 +1,8 @@
all: moduleEM.o
moduleEM.o: moduleEM.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

46
src/modules/moduleEM.f90 → src/modules/solver/electromagnetic/moduleEM.f90
View file

@ -19,6 +19,7 @@ MODULE moduleEM
REAL(8), ALLOCATABLE:: qSpecies(:)
CONTAINS
!Apply boundary conditions to the K matrix for Poisson's equation
SUBROUTINE apply(self, edge)
USE moduleMesh
IMPLICIT NONE
@ -79,6 +80,7 @@ MODULE moduleEM
END FUNCTION gatherMagnField
!Assemble the source vector based on the charge density to solve Poisson's equation
SUBROUTINE assembleSourceVector(vectorF)
USE moduleMesh
USE moduleRefParam
@ -142,4 +144,48 @@ MODULE moduleEM
END SUBROUTINE assembleSourceVector
!Solving the Poisson equation for electrostatic potential
SUBROUTINE solveElecField()
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
INTEGER, SAVE:: INFO
INTEGER:: n
REAL(8), ALLOCATABLE, SAVE:: tempF(:)
EXTERNAL:: dgetrs
!$OMP SINGLE
ALLOCATE(tempF(1:mesh%numNodes))
!$OMP END SINGLE
CALL assembleSourceVector(tempF)
!$OMP SINGLE
CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
mesh%IPIV, tempF, mesh%numNodes, info)
!$OMP END SINGLE
IF (info == 0) THEN
!Suscessful resolution of Poission equation
!$OMP DO
DO n = 1, mesh%numNodes
mesh%nodes(n)%obj%emData%phi = tempF(n)
END DO
!$OMP END DO
ELSE
!$OMP SINGLE
CALL criticalError('Poisson equation failed', 'solveElecField')
!$OMP END SINGLE
END IF
!$OMP SINGLE
DEALLOCATE(tempF)
!$OMP END SINGLE
END SUBROUTINE solveElecField
END MODULE moduleEM

13
src/modules/solver/makefile Normal file
View file

@ -0,0 +1,13 @@
all: moduleSolver.o electromagnetic.o pusher.o
electromagnetic.o:
$(MAKE) -C electromagnetic all
pusher.o:
$(MAKE) -C pusher all
moduleSolver.o: electromagnetic.o pusher.o moduleSolver.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

285
src/modules/moduleSolver.f90 → src/modules/solver/moduleSolver.f90
View file

@ -62,6 +62,7 @@ MODULE moduleSolver
!Init Pusher
SUBROUTINE initPusher(self, pusherType, tau, tauSp)
USE moduleErrors
USE modulePusher
USE moduleMesh, ONLY: mesh
IMPLICIT NONE
@ -128,6 +129,7 @@ MODULE moduleSolver
END SUBROUTINE initPusher
SUBROUTINE initEM(self, EMType)
USE moduleEM
IMPLICIT NONE
CLASS(solverGeneric), INTENT(inout):: self
@ -184,239 +186,6 @@ MODULE moduleSolver
END SUBROUTINE doPushes
!Push neutral particles in cartesian coordinates
PURE SUBROUTINE pushCartNeutral(part, tauIn)
USE moduleSPecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
part%r = part%r + part%v*tauIn
END SUBROUTINE pushCartNeutral
PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: qmEFt(1:3)
!Get the electric field at particle position
qmEFt = part%species%qm*gatherElecField(part)*tauIn
!Update velocity
part%v = part%v + qmEFt
!Update position
part%r = part%r + part%v*tauIn
END SUBROUTINE pushCartElectrostatic
PURE SUBROUTINE pushCartElectromagnetic(part, tauIn)
USE moduleSPecies
USE moduleEM
USE moduleMath
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: tauInHalf
REAL(8):: qmEFt(1:3)
REAL(8):: B(1:3), BNorm
REAL(8):: fn
REAL(8):: v_minus(1:3), v_prime(1:3), v_plus(1:3)
tauInHalf = tauIn *0.5D0
!Half of the force o f the electric field
qmEFt = part%species%qm*gatherElecField(part)*tauInHalf
!Half step for electrostatic
v_minus = part%v + qmEFt
!Full step rotation
B = gatherMagnField(part)
BNorm = NORM2(B)
IF (BNorm > 0.D0) THEN
fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
v_prime = v_minus + fn * crossProduct(v_minus, B)
v_plus = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)
END IF
!Half step for electrostatic
part%v = v_plus + qmEFt
!Update position
part%r = part%r + part%v*tauIn
END SUBROUTINE pushCartElectromagnetic
!Push one particle. Boris pusher for 2D Cyl Neutral particle
PURE SUBROUTINE push2DCylNeutral(part, tauIn)
USE moduleSpecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: v_p_oh_star(2:3)
part_temp = part
!z
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!r,theta
v_p_oh_star(2) = part%v(2)
x_new = part%r(2) + v_p_oh_star(2)*tauIn
v_p_oh_star(3) = part%v(3)
y_new = v_p_oh_star(3)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(2) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push2DCylNeutral
!Push one particle. Boris pusher for 2D Cyl Electrostatic particle
PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: v_p_oh_star(2:3)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass
part_temp = part
!Get electric field at particle position
qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn
!z
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!r,theta
v_p_oh_star(2) = part%v(2) + qmEFt(2)
x_new = part%r(2) + v_p_oh_star(2)*tauIn
v_p_oh_star(3) = part%v(3) + qmEFt(3)
y_new = v_p_oh_star(3)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(2) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push2DCylElectrostatic
!Push one particle. Boris pusher for 1D Radial Neutral particle
PURE SUBROUTINE push1DRadNeutral(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: v_p_oh_star(1:2)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
part_temp = part
!r,theta
v_p_oh_star(1) = part%v(1)
x_new = part%r(1) + v_p_oh_star(1)*tauIn
v_p_oh_star(2) = part%v(2)
y_new = v_p_oh_star(2)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(1) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(1) = cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push1DRadNeutral
!Push one particle. Boris pusher for 1D Radial Electrostatic particle
PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: v_p_oh_star(1:2)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass
part_temp = part
!Get electric field at particle position
qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin
!r,theta
v_p_oh_star(1) = part%v(1) + qmEFt(1)
x_new = part%r(1) + v_p_oh_star(1)*tauIn
v_p_oh_star(2) = part%v(2) + qmEFt(2)
y_new = v_p_oh_star(2)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(1) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(1) = cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push1DRadElectrostatic
!Dummy pusher for 0D geometry
PURE SUBROUTINE push0D(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
END SUBROUTINE push0D
!Takes the particles from a list and put them into an array
!nStart indicates the last fill index in the array
SUBROUTINE resetList(partList, partArray, nStart)
@ -606,55 +375,11 @@ MODULE moduleSolver
END SUBROUTINE doEMField
!Solving the Poisson equation for electrostatic potential
SUBROUTINE solveElecField()
USE moduleEM
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
INTEGER, SAVE:: INFO
INTEGER:: n
REAL(8), ALLOCATABLE, SAVE:: tempF(:)
EXTERNAL:: dgetrs
!$OMP SINGLE
ALLOCATE(tempF(1:mesh%numNodes))
!$OMP END SINGLE
CALL assembleSourceVector(tempF)
!$OMP SINGLE
CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
mesh%IPIV, tempF, mesh%numNodes, info)
!$OMP END SINGLE
IF (info == 0) THEN
!Suscessful resolution of Poission equation
!$OMP DO
DO n = 1, mesh%numNodes
mesh%nodes(n)%obj%emData%phi = tempF(n)
END DO
!$OMP END DO
ELSE
!$OMP SINGLE
CALL criticalError('Poisson equation failed', 'solveElecField')
!$OMP END SINGLE
END IF
!$OMP SINGLE
DEALLOCATE(tempF)
!$OMP END SINGLE
END SUBROUTINE solveElecField
!Modify particle weight as a function of cell volume and splits particle
!Split particles as a function of cell volume and splits particle
SUBROUTINE volumeWScheme(part, volOld, volNew)
USE moduleSpecies
USE moduleMesh
USE moduleRandom
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
@ -669,7 +394,7 @@ MODULE moduleSolver
!Calculate probability of splitting particle
pSplit = 1.D0 - DEXP(-fractionVolume)
IF (random() < pSplit THEN
IF (random() < pSplit) THEN
!Split particle in two
CALL splitParticle(part, 2, volNew)

4
src/modules/solver/pusher/makefile Normal file
View file

@ -0,0 +1,4 @@
all: modulePusher.o
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

237
src/modules/solver/pusher/modulePusher.f90 Normal file
View file

@ -0,0 +1,237 @@
MODULE modulePusher
CONTAINS
!Push neutral particles in cartesian coordinates
PURE SUBROUTINE pushCartNeutral(part, tauIn)
USE moduleSpecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
part%r = part%r + part%v*tauIn
END SUBROUTINE pushCartNeutral
PURE SUBROUTINE pushCartElectrostatic(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: qmEFt(1:3)
!Get the electric field at particle position
qmEFt = part%species%qm*gatherElecField(part)*tauIn
!Update velocity
part%v = part%v + qmEFt
!Update position
part%r = part%r + part%v*tauIn
END SUBROUTINE pushCartElectrostatic
PURE SUBROUTINE pushCartElectromagnetic(part, tauIn)
USE moduleSPecies
USE moduleEM
USE moduleMath
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: tauInHalf
REAL(8):: qmEFt(1:3)
REAL(8):: B(1:3), BNorm
REAL(8):: fn
REAL(8):: v_minus(1:3), v_prime(1:3), v_plus(1:3)
tauInHalf = tauIn *0.5D0
!Half of the force o f the electric field
qmEFt = part%species%qm*gatherElecField(part)*tauInHalf
!Half step for electrostatic
v_minus = part%v + qmEFt
!Full step rotation
B = gatherMagnField(part)
BNorm = NORM2(B)
IF (BNorm > 0.D0) THEN
fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
v_prime = v_minus + fn * crossProduct(v_minus, B)
v_plus = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)
END IF
!Half step for electrostatic
part%v = v_plus + qmEFt
!Update position
part%r = part%r + part%v*tauIn
END SUBROUTINE pushCartElectromagnetic
!Push one particle. Boris pusher for 2D Cyl Neutral particle
PURE SUBROUTINE push2DCylNeutral(part, tauIn)
USE moduleSpecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: v_p_oh_star(2:3)
part_temp = part
!z
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!r,theta
v_p_oh_star(2) = part%v(2)
x_new = part%r(2) + v_p_oh_star(2)*tauIn
v_p_oh_star(3) = part%v(3)
y_new = v_p_oh_star(3)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(2) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push2DCylNeutral
!Push one particle. Boris pusher for 2D Cyl Electrostatic particle
PURE SUBROUTINE push2DCylElectrostatic(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: v_p_oh_star(2:3)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass
part_temp = part
!Get electric field at particle position
qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn
!z
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!r,theta
v_p_oh_star(2) = part%v(2) + qmEFt(2)
x_new = part%r(2) + v_p_oh_star(2)*tauIn
v_p_oh_star(3) = part%v(3) + qmEFt(3)
y_new = v_p_oh_star(3)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(2) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push2DCylElectrostatic
!Push one particle. Boris pusher for 1D Radial Neutral particle
PURE SUBROUTINE push1DRadNeutral(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: v_p_oh_star(1:2)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
part_temp = part
!r,theta
v_p_oh_star(1) = part%v(1)
x_new = part%r(1) + v_p_oh_star(1)*tauIn
v_p_oh_star(2) = part%v(2)
y_new = v_p_oh_star(2)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(1) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(1) = cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push1DRadNeutral
!Push one particle. Boris pusher for 1D Radial Electrostatic particle
PURE SUBROUTINE push1DRadElectrostatic(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
REAL(8):: v_p_oh_star(1:2)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass
part_temp = part
!Get electric field at particle position
qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin
!r,theta
v_p_oh_star(1) = part%v(1) + qmEFt(1)
x_new = part%r(1) + v_p_oh_star(1)*tauIn
v_p_oh_star(2) = part%v(2) + qmEFt(2)
y_new = v_p_oh_star(2)*tauIn
r = DSQRT(x_new**2+y_new**2)
part_temp%r(1) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(1) = cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2)
part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2)
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push1DRadElectrostatic
!Dummy pusher for 0D geometry
PURE SUBROUTINE push0D(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
END SUBROUTINE push0D
END MODULE modulePusher