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

Merge branch 'feature/2DCart' into 'development'

Browse source 
First implementation of 2D Cartesian space.

See merge request JorgeGonz/fpakc!5
This commit is contained in:
Jorge Gonzalez 2021-01-21 11:07:31 +00:00
commit 81e2202959
14 changed files with 2110 additions and 162 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
doc/user-manual/fpakc_UserManual.pdf
View file

Binary file not shown.

35
doc/user-manual/fpakc_UserManual.tex
View file

@ -128,15 +128,18 @@
When a 2D cylindrical geometry is used ($z$, $r$), a Boris solver\cite{boris1970relativistic} is used to move particles accounting for the effect of the symmetry axis.
This pusher removes the issue with particles going to infinite velocity when $r \rightarrow 0$ by pushing the particles in Cartesian space and then converting it to $r-z$ geometry.
Velocity in the $\theta$ direction is updated for collision processes, although the dynamic in the angular direction is assumed as symmetric.
Cross-sections are read from files.
These files contains two columns: one for the relative energy (in $\unit{eV}$) and another one for the cross-section (in $\unit{m^{-2}}$).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{1D Cartesian pusher}
\subsection{2D Cartesian pusher}
Moving particles in a simple 2D Cartesian space.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{1D Radial pusher}
Same implementation as 2D cylindrical pusher but direction $z$ is ignored.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{1D Cartesian pusher}
Moving particles in a simple 1D Cartesian space. Same implementation as in 2D Cartesian but $y$ direction is ignored.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Find new cell}
@ -338,12 +341,14 @@ make
Type of geometry.
Current accepted vaules are
\begin{itemize}
\item \textbf{2DCyl}: Two-dimensional grid (z-r) with symmetry axis at $r = 0$.
\item \textbf{2DCyl}: Two-dimensional grid ($z \hyphen r$) with symmetry axis at $r = 0$.
For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $z$ and $r$ respectively.
\item \textbf{1DCart}: One-dimensional grid (x) in Cartesian coordinates
For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $x$.
\item \textbf{1DRad}: One-dimensional grid (r) in radial coordinates
\item \textbf{2DCart}: Two-dimensional grid ($x \hyphen y$) in Cartesian coordinates..
For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $x$ and $y$ respectively.
\item \textbf{1DRad}: One-dimensional grid ($r$) in radial coordinates
For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$.
\item \textbf{1DCart}: One-dimensional grid ($x$) in Cartesian coordinates
For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $x$.
\end{itemize}
\item \textbf{meshType}: Character.
Format of mesh file.
@ -494,10 +499,14 @@ make
Array dimension 'number of species'.
Indicates the type of pusher used for each species:
\begin{itemize}
\item \textbf{2DCylNeutral}: Pushes particles in a 2D z-r space without any external force.
\item \textbf{2DCylCharged}: Pushes particles in a 2D z-r space including the effect of the electrostatic field.
\item \textbf{1DCartCharged}: Pushes particles in a 1D Cartesian space accounting the the electrostatic field.
\item \textbf{1DRadCharged}: Pushes particles in a 1D cylindrical space (r) accounting the the electrostatic field.
\item \textbf{2DCylNeutral}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) without any external force.
\item \textbf{2DCylCharged}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) including the effect of the electrostatic field.
\item \textbf{2DCartNeutral}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) without any external force.
\item \textbf{2DCartCharged}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) including the effect of the electrostatic field.
\item \textbf{1DRadNeutral}: Pushes particles in a 1D cylindrical space ($r$) without any external force.
\item \textbf{1DRadCharged}: Pushes particles in a 1D cylindrical space ($r$) accounting the the electrostatic field.
\item \textbf{1DCartNeutral}: Pushes particles in a 1D Cartesian space ($x$) without any external force.
\item \textbf{1DCartCharged}: Pushes particles in a 1D Cartesian space ($x$) accounting the the electrostatic field.
\end{itemize}
\item \textbf{WeightingScheme}: Character.
Indicates the variable weighting scheme to be used in the simulation.

3
src/makefile
View file

@ -4,7 +4,8 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolv
$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \
$(OBJDIR)/moduleMeshCyl.o $(OBJDIR)/moduleMeshCylRead.o $(OBJDIR)/moduleMeshCylBoundary.o \
$(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o $(OBJDIR)/moduleMesh2DCylBoundary.o \
$(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o $(OBJDIR)/moduleMesh2DCartBoundary.o \
$(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o $(OBJDIR)/moduleMesh1DCartBoundary.o \
$(OBJDIR)/moduleMesh1DRad.o $(OBJDIR)/moduleMesh1DRadRead.o $(OBJDIR)/moduleMesh1DRadBoundary.o

11
src/modules/mesh/2DCart/makefile Normal file
View file

@ -0,0 +1,11 @@
all : moduleMesh2DCart.o moduleMesh2DCartBoundary.o moduleMesh2DCartRead.o
moduleMesh2DCart.o: moduleMesh2DCart.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh2DCartBoundary.o: moduleMesh2DCart.o moduleMesh2DCartBoundary.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh2DCartRead.o: moduleMesh2DCart.o moduleMesh2DCartBoundary.o moduleMesh2DCartRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

1057
src/modules/mesh/2DCart/moduleMesh2DCart.f90 Normal file
View file

File diff suppressed because it is too large Load diff

154
src/modules/mesh/2DCart/moduleMesh2DCartBoundary.f90 Normal file
View file

@ -0,0 +1,154 @@
!moduleMesh2DCartBoundary: Boundary functions for cylindrical coordinates
SUBMODULE (moduleMesh2DCart) moduleMesh2DCartBoundary
USE moduleMesh2DCart
CONTAINS
SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
!TODO: Try to do this without select
SELECT TYPE(edge)
TYPE IS(meshEdge2DCart)
edgeNorm = DSQRT((edge%y(2)-edge%y(1))**2 + (edge%x(2)-edge%x(1))**2)
cosT = (edge%x(2)-edge%x(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - edge%x(1);
rp(2) = part%r(2) - edge%y(1);
rpp(1) = cosT*rp(1) - sinT*rp(2)
rpp(2) = sinT*rp(1) + cosT*rp(2)
rpp(2) = -rpp(2)
vpp(1) = cosT*part%v(1) - sinT*part%v(2)
vpp(2) = sinT*part%v(1) + cosT*part%v(2)
vpp(2) = -vpp(2)
part%r(1) = cosT*rpp(1) + sinT*rpp(2) + edge%x(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%y(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
END SELECT
part%n_in = .TRUE.
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: rEdge(1:2) !Position of particle projected to the edge
REAL(8):: a, b, c
REAL(8):: a2b2
REAL(8):: d !Distance from particle to edge
SELECT TYPE(edge)
TYPE IS(meshEdge2DCart)
a = (edge%x(1) - edge%x(2))
b = (edge%y(1) - edge%y(2))
c = edge%x(1)*edge%y(2) - edge%x(2)*edge%y(1)
a2b2 = a**2 + b**2
rEdge(1) = (b*( b*part%r(1) - a*part%r(2)) - a*c)/a2b2
rEdge(2) = (a*(-b*part%r(1) + a*part%r(2)) - b*c)/a2b2
d = NORM2(rEdge - part%r(1:2))
!Reduce weight of particle by the distance to the edge and move it to the edge
IF (d > 0.D0) THEN
part%weight = part%weight / d
part%r(1:2) = rEdge
END IF
!Scatter particle in associated volume
IF (ASSOCIATED(edge%e1)) THEN
CALL edge%e1%scatter(part)
ELSE
CALL edge%e2%scatter(part)
END IF
END SELECT
!Remove particle from the domain
part%n_in = .FALSE.
END SUBROUTINE absorption
!Transparent boundary condition
SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
!Removes particle from domain
part%n_in = .FALSE.
END SUBROUTINE transparent
!Wall with temperature
SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
USE moduleBoundary
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
INTEGER:: i
!Modifies particle velocity according to wall temperature
SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
TYPE IS(boundaryWallTemperature)
DO i = 1, 3
part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
END DO
END SELECT
!Reflects particle in the edge
SELECT TYPE(edge)
TYPE IS(meshEdge2DCart)
edgeNorm = DSQRT((edge%y(2)-edge%y(1))**2 + (edge%x(2)-edge%x(1))**2)
cosT = (edge%x(2)-edge%x(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - edge%x(1);
rp(2) = part%r(2) - edge%y(1);
rpp(1) = cosT*rp(1) - sinT*rp(2)
rpp(2) = sinT*rp(1) + cosT*rp(2)
rpp(2) = -rpp(2)
vpp(1) = cosT*part%v(1) - sinT*part%v(2)
vpp(2) = sinT*part%v(1) + cosT*part%v(2)
vpp(2) = -vpp(2)
part%r(1) = cosT*rpp(1) + sinT*rpp(2) + edge%x(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%y(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
END SELECT
part%n_in = .TRUE.
END SUBROUTINE wallTemperature
END SUBMODULE moduleMesh2DCartBoundary

599
src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 Normal file
View file

@ -0,0 +1,599 @@
MODULE moduleMesh2DCartRead
USE moduleMesh
USE moduleMesh2DCart
TYPE, EXTENDS(meshGeneric):: mesh2DCartGeneric
CONTAINS
PROCEDURE, PASS:: init => init2DCartMesh
PROCEDURE, PASS:: readMesh => readMesh2DCartGmsh
END TYPE
INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connected
CONTAINS
!Init mesh
SUBROUTINE init2DCartMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(mesh2DCartGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
SELECT CASE(meshFormat)
CASE ("gmsh")
self%printOutput => printOutputGmsh
self%printColl => printCollGmsh
self%printEM => printEMGmsh
CASE DEFAULT
CALL criticalError("Mesh type " // meshFormat // " not supported.", "init2DCartMesh")
END SELECT
END SUBROUTINE init2DCartMesh
!Read mesh from gmsh file
SUBROUTINE readMesh2DCartGmsh(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(mesh2DCartGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: x, y
INTEGER:: p(1:4)
INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Read msh
OPEN(10, FILE=TRIM(filename))
!Skip header
READ(10, *)
READ(10, *)
READ(10, *)
READ(10, *)
!Read number of nodes
READ(10, *) self%numNodes
!Allocate required matrices and vectors
ALLOCATE(self%nodes(1:self%numNodes))
ALLOCATE(self%K(1:self%numNodes,1:self%numNodes))
ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
self%K = 0.D0
self%IPIV = 0
!Read nodes cartesian coordinates (x=x, y=y, z=null)
DO e=1, self%numNodes
READ(10, *) n, x, y
ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/x, y, 0.D0 /))
END DO
!Skips comments
READ(10, *)
READ(10, *)
!Reads Totalnumber of elements
READ(10, *) TotalnumElem
!counts edges and volume elements
self%numEdges = 0
DO e=1, TotalnumElem
READ(10,*) eTemp, elemType
IF (elemType==1) THEN
self%numEdges=e
END IF
END DO
!Substract the number of edges to the total number of elements
!to obtain the number of volume elements
self%numVols = TotalnumElem - self%numEdges
!Allocates arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beggining to read elements
DO e=1, totalNumElem
BACKSPACE(10)
END DO
!Reads edges
DO e=1, self%numEdges
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType)
END DO
!Read and initialize volumes
DO e=1, self%numVols
READ(10,*) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE (2)
!Triangular element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVol2DCartTria:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3))
CASE (3)
!Quadrilateral element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVol2DCartQuad:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
END SELECT
END DO
CLOSE(10)
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
DO et = 1, self%numVols
IF (e /= et) THEN
CALL connected(self%vols(e)%obj, self%vols(et)%obj)
END IF
END DO
!Connectivity between vols and edges
DO et = 1, self%numEdges
CALL connected(self%vols(e)%obj, self%edges(et)%obj)
END DO
!Constructs the global K matrix
CALL constructGlobalK(self%K, self%vols(e)%obj)
END DO
END SUBROUTINE readMesh2DCartGmsh
!Selects type of elements to build connection
SUBROUTINE connectedVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshVol2DCartQuad)
!Element A is a quadrilateral
SELECT TYPE(elemB)
TYPE IS(meshVol2DCartQuad)
!Element B is a quadrilateral
CALL connectedQuadQuad(elemA, elemB)
TYPE IS(meshVol2DCartTria)
!Element B is a triangle
CALL connectedQuadTria(elemA, elemB)
END SELECT
TYPE IS(meshVol2DCartTria)
!Element A is a Triangle
SELECT TYPE(elemB)
TYPE IS(meshVol2DCartQuad)
!Element B is a quadrilateral
CALL connectedQuadTria(elemB, elemA)
TYPE IS(meshVol2DCartTria)
!Element B is a triangle
CALL connectedTriaTria(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolVol
SUBROUTINE connectedVolEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemB)
CLASS IS(meshEdge2DCart)
SELECT TYPE(elemA)
TYPE IS(meshVol2DCartQuad)
!Element A is a quadrilateral
CALL connectedQuadEdge(elemA, elemB)
TYPE IS(meshVol2DCartTria)
!Element A is a triangle
CALL connectedTriaEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
SUBROUTINE connectedQuadQuad(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n4%n) THEN
elemA%e2 => elemB
elemB%e4 => elemA
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
END IF
END SUBROUTINE connectedQuadQuad
SUBROUTINE connectedQuadTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e4 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedQuadTria
SUBROUTINE connectedTriaTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA
CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
END IF
END IF
END SUBROUTINE connectedTriaTria
SUBROUTINE connectedQuadEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA
CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedQuadEdge
SUBROUTINE connectedTriaEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol2DCartTria), INTENT(inout), TARGET:: elemA
CLASS(meshEdge2DCart), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedTriaEdge
SUBROUTINE constructGlobalK(K, elem)
IMPLICIT NONE
REAL(8), INTENT(inout):: K(1:,1:)
CLASS(meshVol), INTENT(in):: elem
REAL(8), ALLOCATABLE:: localK(:,:)
INTEGER:: nNodes, i, j
INTEGER, ALLOCATABLE:: n(:)
SELECT TYPE(elem)
TYPE IS(meshVol2DCartQuad)
nNodes = 4
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, &
elem%n3%n, elem%n4%n /)
TYPE IS(meshVol2DCartTria)
nNodes = 3
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, elem%n3%n /)
CLASS DEFAULT
nNodes = 0
ALLOCATE(localK(1:1, 1:1))
localK = 0.D0
ALLOCATE(n(1:1))
n = 0
END SELECT
DO i = 1, nNodes
DO j = 1, nNodes
K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
END DO
END DO
END SUBROUTINE constructGlobalK
END MODULE moduleMesh2DCartRead

8
src/modules/mesh/2DCyl/makefile
View file

@ -1,11 +1,11 @@
all : moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.o
all : moduleMesh2DCyl.o moduleMesh2DCylBoundary.o moduleMesh2DCylRead.o
moduleMeshCyl.o: moduleMeshCyl.f90
moduleMesh2DCyl.o: moduleMesh2DCyl.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f90
moduleMesh2DCylBoundary.o: moduleMesh2DCyl.o moduleMesh2DCylBoundary.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMeshCylRead.o: moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f90
moduleMesh2DCylRead.o: moduleMesh2DCyl.o moduleMesh2DCylBoundary.o moduleMesh2DCylRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

158
src/modules/mesh/2DCyl/moduleMeshCyl.f90 → src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -1,8 +1,8 @@
!moduleMeshCyl: 2D axial symmetric extension of generic mesh from GMSH format.
!moduleMesh2DCyl: 2D axial symmetric extension of generic mesh from GMSH format.
! x == z
! y == r
! z == theta (unused)
MODULE moduleMeshCyl
MODULE moduleMesh2DCyl
USE moduleMesh
IMPLICIT NONE
@ -14,26 +14,26 @@ MODULE moduleMeshCyl
REAL(8), PARAMETER:: xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /)
REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /)
TYPE, PUBLIC, EXTENDS(meshNode):: meshNodeCyl
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCyl
!Element coordinates
REAL(8):: r = 0.D0, z = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initNodeCyl
PROCEDURE, PASS:: getCoordinates => getCoordCyl
PROCEDURE, PASS:: init => initNode2DCyl
PROCEDURE, PASS:: getCoordinates => getCoord2DCyl
END TYPE meshNodeCyl
END TYPE meshNode2DCyl
TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdgeCyl
TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge2DCyl
!Element coordinates
REAL(8):: r(1:2) = 0.D0, z(1:2) = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
CONTAINS
PROCEDURE, PASS:: init => initEdgeCyl
PROCEDURE, PASS:: getNodes => getNodesCyl
PROCEDURE, PASS:: init => initEdge2DCyl
PROCEDURE, PASS:: getNodes => getNodes2DCyl
PROCEDURE, PASS:: randPos => randPosEdge
END TYPE meshEdgeCyl
END TYPE meshEdge2DCyl
!Boundary functions defined in the submodule Boundary
INTERFACE
@ -84,15 +84,15 @@ MODULE moduleMeshCyl
END INTERFACE
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCyl
CONTAINS
PROCEDURE, PASS:: detJac => detJCyl
PROCEDURE, PASS:: invJac => invJCyl
PROCEDURE, PASS:: detJac => detJ2DCyl
PROCEDURE, PASS:: invJac => invJ2DCyl
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVolCyl
END TYPE meshVol2DCyl
ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
@ -108,8 +108,8 @@ MODULE moduleMeshCyl
END FUNCTION dPsi_interface
PURE SUBROUTINE partialDer_interface(self, dPsi, dz, dr)
IMPORT meshVolCyl
CLASS(meshVolCyl), INTENT(in):: self
IMPORT meshVol2DCyl
CLASS(meshVol2DCyl), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
@ -118,7 +118,7 @@ MODULE moduleMeshCyl
END INTERFACE
!Quadrilateral volume element
TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylQuad
TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylQuad
!Element coordinates
REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0
!Connectivity to nodes
@ -128,7 +128,7 @@ MODULE moduleMeshCyl
REAL(8):: arNodes(1:4) = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initVolQuadCyl
PROCEDURE, PASS:: init => initVolQuad2DCyl
PROCEDURE, PASS:: randPos => randPosVolQuad
PROCEDURE, PASS:: area => areaQuad
PROCEDURE, NOPASS:: fPsi => fPsiQuad
@ -146,10 +146,10 @@ MODULE moduleMeshCyl
PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: nextElement => nextElementQuad
END TYPE meshVolCylQuad
END TYPE meshVol2DCylQuad
!Triangular volume element
TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylTria
TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylTria
!Element coordinates
REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0
!Connectivity to nodes
@ -157,11 +157,9 @@ MODULE moduleMeshCyl
!Connectivity to adjacent elements
CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
REAL(8):: arNodes(1:3) = 0.D0
!Derivatives in z,r real space
REAL(8), DIMENSION(1:3):: dPsiZ, dPsiR
CONTAINS
PROCEDURE, PASS:: init => initVolTriaCyl
PROCEDURE, PASS:: init => initVolTria2DCyl
PROCEDURE, PASS:: randPos => randPosVolTria
PROCEDURE, PASS:: area => areaTria
PROCEDURE, NOPASS:: fPsi => fPsiTria
@ -179,51 +177,51 @@ MODULE moduleMeshCyl
PROCEDURE, PASS:: phy2log => phy2logTria
PROCEDURE, PASS:: nextElement => nextElementTria
END TYPE meshVolCylTria
END TYPE meshVol2DCylTria
CONTAINS
!NODE FUNCTIONS
!Inits node element
SUBROUTINE initNodeCyl(self, n, r)
SUBROUTINE initNode2DCyl(self, n, r)
USE moduleSpecies
USE moduleRefParam
IMPLICIT NONE
CLASS(meshNodeCyl), INTENT(out):: self
CLASS(meshNode2DCyl), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
self%n = n
self%z = r(2)/L_ref
self%r = r(1)/L_ref
self%z = r(1)/L_ref
self%r = r(2)/L_ref
!Node volume, to be determined in mesh
self%v = 0.D0
!Allocates output:
ALLOCATE(self%output(1:nSpecies))
END SUBROUTINE initNodeCyl
END SUBROUTINE initNode2DCyl
!Get coordinates from node
PURE FUNCTION getCoordCyl(self) RESULT(r)
PURE FUNCTION getCoord2DCyl(self) RESULT(r)
IMPLICIT NONE
CLASS(meshNodeCyl), INTENT(in):: self
CLASS(meshNode2DCyl), INTENT(in):: self
REAL(8):: r(1:3)
r = (/self%z, self%r, 0.D0/)
END FUNCTION getCoordCyl
END FUNCTION getCoord2DCyl
!EDGE FUNCTIONS
!Inits edge element
SUBROUTINE initEdgeCyl(self, n, p, bt, physicalSurface)
SUBROUTINE initEdge2DCyl(self, n, p, bt, physicalSurface)
USE moduleSpecies
USE moduleBoundary
USE moduleErrors
IMPLICIT NONE
CLASS(meshEdgeCyl), INTENT(out):: self
CLASS(meshEdge2DCyl), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
@ -265,7 +263,7 @@ MODULE moduleMeshCyl
self%fBoundary(s)%apply => symmetryAxis
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'initEdgeCyl')
CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCyl')
END SELECT
@ -274,14 +272,14 @@ MODULE moduleMeshCyl
!Physical surface
self%physicalSurface = physicalSurface
END SUBROUTINE initEdgeCyl
END SUBROUTINE initEdge2DCyl
!Random position in quadrilateral volume
FUNCTION randPosVolQuad(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
@ -299,23 +297,23 @@ MODULE moduleMeshCyl
END FUNCTION randposVolQuad
!Get nodes from edge
PURE FUNCTION getNodesCyl(self) RESULT(n)
PURE FUNCTION getNodes2DCyl(self) RESULT(n)
IMPLICIT NONE
CLASS(meshEdgeCyl), INTENT(in):: self
CLASS(meshEdge2DCyl), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1:2))
n = (/self%n1%n, self%n2%n /)
END FUNCTION getNodesCyl
END FUNCTION getNodes2DCyl
!Calculates a random position in edge
FUNCTION randPosEdge(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdgeCyl), INTENT(in):: self
CLASS(meshEdge2DCyl), INTENT(in):: self
REAL(8):: rnd
REAL(8):: r(1:3)
REAL(8):: p1(1:2), p2(1:2)
@ -332,11 +330,11 @@ MODULE moduleMeshCyl
!VOLUME FUNCTIONS
!QUAD FUNCTIONS
!Inits quadrilateral element
SUBROUTINE initVolQuadCyl(self, n, p)
SUBROUTINE initVolQuad2DCyl(self, n, p)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(out):: self
CLASS(meshVol2DCylQuad), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
REAL(8), DIMENSION(1:3):: r1, r2, r3, r4
@ -365,14 +363,14 @@ MODULE moduleMeshCyl
CALL OMP_INIT_LOCK(self%lock)
END SUBROUTINE initVolQuadCyl
END SUBROUTINE initVolQuad2DCyl
!Computes element area
PURE SUBROUTINE areaQuad(self)
USE moduleConstParam, ONLY: PI8
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout):: self
CLASS(meshVol2DCylQuad), INTENT(inout):: self
REAL(8):: r, xi(1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:4)
@ -453,7 +451,7 @@ MODULE moduleMeshCyl
PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
@ -469,7 +467,7 @@ MODULE moduleMeshCyl
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8):: r, xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:2,1:4)
REAL(8):: ke(1:4,1:4)
@ -502,7 +500,7 @@ MODULE moduleMeshCyl
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: r, xi(1:3)
@ -558,7 +556,7 @@ MODULE moduleMeshCyl
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
CLASS(particle), INTENT(in):: part
TYPE(outputNode), POINTER:: vertex
REAL(8):: w_p(1:4)
@ -593,7 +591,7 @@ MODULE moduleMeshCyl
PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1:2,1:4)
REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates
@ -620,7 +618,7 @@ MODULE moduleMeshCyl
PURE FUNCTION getNodesQuad(self) RESULT(n)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1:4))
@ -632,7 +630,7 @@ MODULE moduleMeshCyl
PURE FUNCTION phy2logQuad(self,r) RESULT(xN)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2)
@ -662,7 +660,7 @@ MODULE moduleMeshCyl
SUBROUTINE nextElementQuad(self, xi, nextElement)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(in):: self
CLASS(meshVol2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), POINTER, INTENT(out):: nextElement
REAL(8):: xiArray(1:4)
@ -687,15 +685,14 @@ MODULE moduleMeshCyl
!TRIA ELEMENT
!Init tria element
SUBROUTINE initVolTriaCyl(self, n, p)
SUBROUTINE initVolTria2DCyl(self, n, p)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(out):: self
CLASS(meshVol2DCylTria), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
REAL(8), DIMENSION(1:3):: r1, r2, r3
REAL(8):: A
!Assign node index
self%n = n
@ -716,31 +713,18 @@ MODULE moduleMeshCyl
self%n2%v = self%n2%v + self%arNodes(2)
self%n3%v = self%n3%v + self%arNodes(3)
!Derivatives in z/r for shape functions (node independent)
A = self%z(2)*self%r(3) - self%z(3)*self%r(2) + &
self%z(3)*self%r(1) - self%z(1)*self%r(3) + &
self%z(1)*self%r(2) - self%z(2)*self%r(1)
self%dPsiZ = (/ self%r(2)-self%r(3), &
self%r(3)-self%r(1), &
self%r(1)-self%r(2) /)/A
self%dPsiR = (/ self%z(3)-self%z(2), &
self%z(1)-self%z(3), &
self%z(2)-self%z(1) /)/A
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
END SUBROUTINE initVolTriaCyl
END SUBROUTINE initVolTria2DCyl
!Random position in quadrilateral volume
FUNCTION randPosVolTria(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
@ -762,7 +746,7 @@ MODULE moduleMeshCyl
USE moduleConstParam, ONLY: PI
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(inout):: self
CLASS(meshVol2DCylTria), INTENT(inout):: self
REAL(8):: r, xi(1:3)
REAL(8):: detJ
REAL(8):: fPsi(1:3)
@ -837,7 +821,7 @@ MODULE moduleMeshCyl
PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
@ -853,7 +837,7 @@ MODULE moduleMeshCyl
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8):: r, xi(1:3)
REAL(8):: fPsi(1:3), dPsi(1:2,1:3)
REAL(8):: ke(1:3,1:3)
@ -883,7 +867,7 @@ MODULE moduleMeshCyl
USE moduleConstParam, ONLY: PI2
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:3)
@ -938,7 +922,7 @@ MODULE moduleMeshCyl
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
CLASS(particle), INTENT(in):: part
TYPE(outputNode), POINTER:: vertex
REAL(8):: w_p(1:3)
@ -968,7 +952,7 @@ MODULE moduleMeshCyl
PURE FUNCTION gatherEFTria(self,xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1:2,1:3)
REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates
@ -995,7 +979,7 @@ MODULE moduleMeshCyl
PURE FUNCTION getNodesTria(self) RESULT(n)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1:3))
@ -1007,7 +991,7 @@ MODULE moduleMeshCyl
PURE FUNCTION phy2logTria(self,r) RESULT(xi)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xi(1:3)
REAL(8):: invJ(1:2,1:2), detJ
@ -1027,7 +1011,7 @@ MODULE moduleMeshCyl
SUBROUTINE nextElementTria(self, xi, nextElement)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(in):: self
CLASS(meshVol2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), POINTER, INTENT(out):: nextElement
REAL(8):: xiArray(1:3)
@ -1049,10 +1033,10 @@ MODULE moduleMeshCyl
!COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS
!Computes element Jacobian determinant
PURE FUNCTION detJCyl(self, xi, dPsi_in) RESULT(dJ)
PURE FUNCTION detJ2DCyl(self, xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE
CLASS(meshVolCyl), INTENT(in):: self
CLASS(meshVol2DCyl), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
@ -1070,13 +1054,13 @@ MODULE moduleMeshCyl
CALL self%partialDer(dPsi, dz, dr)
dJ = dz(1)*dr(2)-dz(2)*dr(1)
END FUNCTION detJCyl
END FUNCTION detJ2DCyl
!Computes element Jacobian inverse matrix (without determinant)
PURE FUNCTION invJCyl(self,xi,dPsi_in) RESULT(invJ)
PURE FUNCTION invJ2DCyl(self,xi,dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshVolCyl), INTENT(in):: self
CLASS(meshVol2DCyl), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
@ -1095,6 +1079,6 @@ MODULE moduleMeshCyl
invJ(1,:) = (/ dr(2), -dz(2) /)
invJ(2,:) = (/ -dr(1), dz(1) /)
END FUNCTION invJCyl
END FUNCTION invJ2DCyl
END MODULE moduleMeshCyl
END MODULE moduleMesh2DCyl

14
src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 → src/modules/mesh/2DCyl/moduleMesh2DCylBoundary.f90
View file

@ -1,6 +1,6 @@
!moduleMeshCylBoundary: Boundary functions for cylindrical coordinates
SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary
USE moduleMeshCyl
!moduleMesh2DCylBoundary: Boundary functions for cylindrical coordinates
SUBMODULE (moduleMesh2DCyl) moduleMesh2DCylBoundary
USE moduleMesh2DCyl
CONTAINS
SUBROUTINE reflection(edge, part)
@ -13,7 +13,7 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary
!TODO: Try to do this without select
SELECT TYPE(edge)
TYPE IS(meshEdgeCyl)
TYPE IS(meshEdge2DCyl)
edgeNorm = DSQRT((edge%r(2)-edge%r(1))**2 + (edge%z(2)-edge%z(1))**2)
cosT = (edge%z(2)-edge%z(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
@ -53,7 +53,7 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary
REAL(8):: d !Distance from particle to edge
SELECT TYPE(edge)
TYPE IS(meshEdgeCyl)
TYPE IS(meshEdge2DCyl)
a = (edge%z(1) - edge%z(2))
b = (edge%r(1) - edge%r(2))
c = edge%z(1)*edge%r(2) - edge%z(2)*edge%r(1)
@ -124,7 +124,7 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary
!Reflects particle in the edge
SELECT TYPE(edge)
TYPE IS(meshEdgeCyl)
TYPE IS(meshEdge2DCyl)
edgeNorm = DSQRT((edge%r(2)-edge%r(1))**2 + (edge%z(2)-edge%z(1))**2)
cosT = (edge%z(2)-edge%z(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
@ -161,4 +161,4 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary
END SUBROUTINE symmetryAxis
END SUBMODULE moduleMeshCylBoundary
END SUBMODULE moduleMesh2DCylBoundary

78
src/modules/mesh/2DCyl/moduleMeshCylRead.f90 → src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90
View file

@ -1,11 +1,11 @@
MODULE moduleMeshCylRead
MODULE moduleMesh2DCylRead
USE moduleMesh
USE moduleMeshCyl
USE moduleMesh2DCyl
TYPE, EXTENDS(meshGeneric):: meshCylGeneric
TYPE, EXTENDS(meshGeneric):: mesh2DCylGeneric
CONTAINS
PROCEDURE, PASS:: init => initCylMesh
PROCEDURE, PASS:: readMesh => readMeshCylGmsh
PROCEDURE, PASS:: init => init2DCylMesh
PROCEDURE, PASS:: readMesh => readMesh2DCylGmsh
END TYPE
@ -16,12 +16,12 @@ MODULE moduleMeshCylRead
CONTAINS
!Init mesh
SUBROUTINE initCylMesh(self, meshFormat)
SUBROUTINE init2DCylMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(out):: self
CLASS(mesh2DCylGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
SELECT CASE(meshFormat)
@ -31,18 +31,18 @@ MODULE moduleMeshCylRead
self%printEM => printEMGmsh
CASE DEFAULT
CALL criticalError("Mesh type " // meshFormat // " not supported.", "initCylMesh")
CALL criticalError("Mesh type " // meshFormat // " not supported.", "init2DCylMesh")
END SELECT
END SUBROUTINE initCylMesh
END SUBROUTINE init2DCylMesh
!Read mesh from gmsh file
SUBROUTINE readMeshCylGmsh(self, filename)
SUBROUTINE readMesh2DCylGmsh(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(inout):: self
CLASS(mesh2DCylGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r, z
INTEGER:: p(1:4)
@ -68,8 +68,8 @@ MODULE moduleMeshCylRead
!Read nodes cartesian coordinates (x=z, y=r, z=null)
DO e=1, self%numNodes
READ(10, *) n, z, r
ALLOCATE(meshNodeCyl:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/r, z, 0.D0 /))
ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/z, r, 0.D0 /))
END DO
!Skips comments
@ -103,7 +103,7 @@ MODULE moduleMeshCylRead
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdgeCyl:: self%edges(e)%obj)
ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType)
@ -118,13 +118,13 @@ MODULE moduleMeshCylRead
CASE (2)
!Triangular element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVolCylTria:: self%vols(e)%obj)
ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3))
CASE (3)
!Quadrilateral element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVolCylQuad:: self%vols(e)%obj)
ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
END SELECT
@ -154,7 +154,7 @@ MODULE moduleMeshCylRead
END DO
END SUBROUTINE readMeshCylGmsh
END SUBROUTINE readMesh2DCylGmsh
!Selects type of elements to build connection
SUBROUTINE connectedVolVol(elemA, elemB)
@ -164,27 +164,27 @@ MODULE moduleMeshCylRead
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshVolCylQuad)
TYPE IS(meshVol2DCylQuad)
!Element A is a quadrilateral
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
TYPE IS(meshVol2DCylQuad)
!Element B is a quadrilateral
CALL connectedQuadQuad(elemA, elemB)
TYPE IS(meshVolCylTria)
TYPE IS(meshVol2DCylTria)
!Element B is a triangle
CALL connectedQuadTria(elemA, elemB)
END SELECT
TYPE IS(meshVolCylTria)
TYPE IS(meshVol2DCylTria)
!Element A is a Triangle
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
TYPE IS(meshVol2DCylQuad)
!Element B is a quadrilateral
CALL connectedQuadTria(elemB, elemA)
TYPE IS(meshVolCylTria)
TYPE IS(meshVol2DCylTria)
!Element B is a triangle
CALL connectedTriaTria(elemA, elemB)
@ -201,13 +201,13 @@ MODULE moduleMeshCylRead
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemB)
CLASS IS(meshEdgeCyl)
CLASS IS(meshEdge2DCyl)
SELECT TYPE(elemA)
TYPE IS(meshVolCylQuad)
TYPE IS(meshVol2DCylQuad)
!Element A is a quadrilateral
CALL connectedQuadEdge(elemA, elemB)
TYPE IS(meshVolCylTria)
TYPE IS(meshVol2DCylTria)
!Element A is a triangle
CALL connectedTriaEdge(elemA, elemB)
@ -220,8 +220,8 @@ MODULE moduleMeshCylRead
SUBROUTINE connectedQuadQuad(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB
CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
@ -264,8 +264,8 @@ MODULE moduleMeshCylRead
SUBROUTINE connectedQuadTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
@ -356,8 +356,8 @@ MODULE moduleMeshCylRead
SUBROUTINE connectedTriaTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
@ -429,8 +429,8 @@ MODULE moduleMeshCylRead
SUBROUTINE connectedQuadEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
@ -501,8 +501,8 @@ MODULE moduleMeshCylRead
SUBROUTINE connectedTriaEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
@ -564,7 +564,7 @@ MODULE moduleMeshCylRead
INTEGER, ALLOCATABLE:: n(:)
SELECT TYPE(elem)
TYPE IS(meshVolCylQuad)
TYPE IS(meshVol2DCylQuad)
nNodes = 4
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
@ -572,7 +572,7 @@ MODULE moduleMeshCylRead
n = (/ elem%n1%n, elem%n2%n, &
elem%n3%n, elem%n4%n /)
TYPE IS(meshVolCylTria)
TYPE IS(meshVol2DCylTria)
nNodes = 3
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
@ -596,4 +596,4 @@ MODULE moduleMeshCylRead
END SUBROUTINE constructGlobalK
END MODULE moduleMeshCylRead
END MODULE moduleMesh2DCylRead

5
src/modules/mesh/makefile
View file

@ -1,8 +1,11 @@
all: moduleMesh.o 2DCyl.o 1DRad.o 1DCart.o
all: moduleMesh.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
2DCyl.o:
$(MAKE) -C 2DCyl all
2DCart.o:
$(MAKE) -C 2DCart all
1DCart.o:
$(MAKE) -C 1DCart all

10
src/modules/moduleInput.f90
View file

@ -365,7 +365,6 @@ MODULE moduleInput
!Gets the basename of the folder
CALL config%get(object // '.folder', baseName, found)
PRINT *, baseName
IF (found) THEN
folder = baseName
@ -654,7 +653,8 @@ MODULE moduleInput
!Read the geometry (mesh) for the case
SUBROUTINE readGeometry(config)
USE moduleMesh
USE moduleMeshCylRead, ONLY: meshCylGeneric
USE moduleMesh2DCylRead, ONLY: mesh2DCylGeneric
USE moduleMesh2DCartRead, ONLY: mesh2DCartGeneric
USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric
USE moduleMesh1DRadRead, ONLY: mesh1DRadGeneric
USE moduleErrors
@ -672,7 +672,11 @@ MODULE moduleInput
SELECT CASE(geometryType)
CASE ("2DCyl")
!Creates a 2D cylindrical mesh
ALLOCATE(meshCylGeneric:: mesh)
ALLOCATE(mesh2DCylGeneric:: mesh)
CASE ("2DCart")
!Creates a 2D cylindrical mesh
ALLOCATE(mesh2DCartGeneric:: mesh)
CASE ("1DCart")
!Creates a 1D cartesian mesh

140
src/modules/moduleSolver.f90
View file

@ -67,15 +67,31 @@ MODULE moduleSolver
REAL(8):: tau, tauSp
SELECT CASE(pusherType)
!2D Cylindrical
CASE('2DCylNeutral')
self%pushParticle => pushCylNeutral
self%pushParticle => push2DCylNeutral
CASE('2DCylCharged')
self%pushParticle => pushCylCharged
self%pushParticle => push2DCylCharged
!2D Cartesian
CASE('2DCartNeutral')
self%pushParticle => push2DCartNeutral
CASE('2DCartCharged')
self%pushParticle => push2DCartCharged
!1D Cartesian
CASE('1DCartNeutral')
self%pushParticle => push1DCartNeutral
CASE('1DCartCharged')
self%pushParticle => push1DCartCharged
!1D Radial
CASE('1DRadNeutral')
self%pushParticle => push1DRadNeutral
CASE('1DRadCharged')
self%pushParticle => push1DRadCharged
@ -147,7 +163,7 @@ MODULE moduleSolver
END SUBROUTINE doPushes
!Push one particle. Boris pusher for 2D Cyl Neutral particle
PURE SUBROUTINE pushCylNeutral(part, tauIn)
PURE SUBROUTINE push2DCylNeutral(part, tauIn)
USE moduleSpecies
IMPLICIT NONE
@ -181,10 +197,10 @@ MODULE moduleSolver
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE pushCylNeutral
END SUBROUTINE push2DCylNeutral
!Push one particle. Boris pusher for 2D Cyl Charged particle
PURE SUBROUTINE pushCylCharged(part, tauIn)
PURE SUBROUTINE push2DCylCharged(part, tauIn)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
@ -222,7 +238,83 @@ MODULE moduleSolver
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE pushCylCharged
END SUBROUTINE push2DCylCharged
!Push neutral particles in 2D cartesian coordinates
PURE SUBROUTINE push2DCartNeutral(part, tauIn)
USE moduleSPecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
part_temp = part
!x
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!y
part_temp%v(2) = part%v(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push2DCartNeutral
!Push charged particles in 2D cartesian coordinates
PURE SUBROUTINE push2DCartCharged(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
REAL(8):: qmEFt(1:3)
part_temp = part
!Get the electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn
!x
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
!y
part_temp%v(2) = part%v(2) + qmEFt(2)
part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push2DCartCharged
!Push neutral particles in 1D cartesian coordinates
PURE SUBROUTINE push1DCartNeutral(part, tauIn)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8), INTENT(in):: tauIn
TYPE(particle):: part_temp
part_temp = part
!x
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push1DCartNeutral
!Push charged particles in 1D cartesian coordinates
PURE SUBROUTINE push1DCartCharged(part, tauIn)
@ -249,6 +341,41 @@ MODULE moduleSolver
END SUBROUTINE push1DCartCharged
!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)
part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE push1DRadNeutral
!Push one particle. Boris pusher for 1D Radial Charged particle
PURE SUBROUTINE push1DRadCharged(part, tauIn)
USE moduleSpecies
@ -263,7 +390,6 @@ MODULE moduleSolver
REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass
part_temp = part
!Time step for the species
!Get electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauMin
!r,theta