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

Now the initial state has the same format as the mesh file, i.e., is

Browse source 
easy to use a file from a previous run without processing it into a
plain text file.

Although the previous method presented some updates for 1D small cases,
this is quite easy to do with any type of simulations and, in the
future, with different mesh formats.
This commit is contained in:
Jorge Gonzalez 2021-04-13 16:55:50 +02:00
commit e25b567d36
10 changed files with 180 additions and 94 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.

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

@ -618,25 +618,12 @@ make
Required values are: Required values are:
\begin{itemize} \begin{itemize}
\item \textbf{speciesName}: Character. \item \textbf{speciesName}: Character.
Name of species. Name of species as defined in the object \textbf{species}.
\item \textbf{initialState}: Character. \item \textbf{initialState}: Character.
Plain text file that contains the information about the species macroscopic properties in the grid. Output file from previous run used as an initial state for the species.
Initial particles are assumed to be Maxwellian. The file format must be the same as in \textbf{geometry.meshType}
Initial particles are assumed to have a Maxwellian distribution.
File must be located at \textbf{output.path}. File must be located at \textbf{output.path}.
The file must contain the following columns in this specific order:
\begin{itemize}
\item \textit{Element}: Integer.
Identifier of the volume in the mesh.
It is not required to specify empty volumes.
\item \textit{Density}: Real.
Species density in $\unit{m^-3}$.
\item \textit{Velocity}: Real.
Three colums representing the initial velocity in each direction.
Units are $\unit{m s^-1}$.
\item \textit{Temperature}: Real.
One column that represents the initial temperature in $\unit{K}$.
\end{itemize}
\end{itemize} \end{itemize}
\end{itemize} \end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7
src/modules/mesh/1DCart/moduleMesh1DCart.f90
View file

@ -36,19 +36,12 @@ MODULE moduleMesh1DCart
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ1DCart PROCEDURE, PASS:: detJac => detJ1DCart
PROCEDURE, PASS:: invJac => invJ1DCart PROCEDURE, PASS:: invJac => invJ1DCart
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol1DCart END TYPE meshVol1DCart
ABSTRACT INTERFACE ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8), ALLOCATABLE:: dPsi(:,:)

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

@ -36,19 +36,12 @@ MODULE moduleMesh1DRad
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ1DRad PROCEDURE, PASS:: detJac => detJ1DRad
PROCEDURE, PASS:: invJac => invJ1DRad PROCEDURE, PASS:: invJac => invJ1DRad
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol1DRad END TYPE meshVol1DRad
ABSTRACT INTERFACE ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8), ALLOCATABLE:: dPsi(:,:)

7
src/modules/mesh/2DCart/moduleMesh2DCart.f90
View file

@ -41,19 +41,12 @@ MODULE moduleMesh2DCart
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ2DCart PROCEDURE, PASS:: detJac => detJ2DCart
PROCEDURE, PASS:: invJac => invJ2DCart PROCEDURE, PASS:: invJac => invJ2DCart
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol2DCart END TYPE meshVol2DCart
ABSTRACT INTERFACE ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8), ALLOCATABLE:: dPsi(:,:)

7
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -41,19 +41,12 @@ MODULE moduleMesh2DCyl
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ2DCyl PROCEDURE, PASS:: detJac => detJ2DCyl
PROCEDURE, PASS:: invJac => invJ2DCyl PROCEDURE, PASS:: invJac => invJ2DCyl
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol2DCyl END TYPE meshVol2DCyl
ABSTRACT INTERFACE ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8), ALLOCATABLE:: dPsi(:,:)

19
src/modules/mesh/3DCart/moduleMesh3DCart.f90
View file

@ -35,19 +35,12 @@ MODULE moduleMesh3DCart
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ3DCart PROCEDURE, PASS:: detJac => detJ3DCart
PROCEDURE, PASS:: invJac => invJ3DCart PROCEDURE, PASS:: invJac => invJ3DCart
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol3DCart END TYPE meshVol3DCart
ABSTRACT INTERFACE ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xii) RESULT(fPsi)
REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION dPsi_interface(xii) RESULT(dPsi) PURE FUNCTION dPsi_interface(xii) RESULT(dPsi)
REAL(8), INTENT(in):: xii(1:3) REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8), ALLOCATABLE:: dPsi(:,:)
@ -327,18 +320,18 @@ MODULE moduleMesh3DCart
END SUBROUTINE volumeTetra END SUBROUTINE volumeTetra
!Computes element functions in point xii !Computes element functions in point xii
PURE FUNCTION fPsiTetra(xii) RESULT(fPsi) PURE FUNCTION fPsiTetra(xi) RESULT(fPsi)
IMPLICIT NONE IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3) REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:) REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:4)) ALLOCATE(fPsi(1:4))
fPsi(1) = 1.D0 - xii(1) - xii(2) - xii(3) fPsi(1) = 1.D0 - xi(1) - xi(2) - xi(3)
fPsi(2) = xii(1) fPsi(2) = xi(1)
fPsi(3) = xii(2) fPsi(3) = xi(2)
fPsi(4) = xii(3) fPsi(4) = xi(3)
END FUNCTION fPsiTetra END FUNCTION fPsiTetra

73
src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
View file

@ -14,6 +14,7 @@ MODULE moduleMeshInputGmsh2
TYPE IS(meshParticles) TYPE IS(meshParticles)
self%printOutput => printOutputGmsh2 self%printOutput => printOutputGmsh2
self%printEM => printEMGmsh2 self%printEM => printEMGmsh2
self%readInitial => readInitialGmsh2
END SELECT END SELECT
self%readMesh => readGmsh2 self%readMesh => readGmsh2
@ -288,4 +289,76 @@ MODULE moduleMeshInputGmsh2
END SUBROUTINE readGmsh2 END SUBROUTINE readGmsh2
!Reads the initial information from an output file for an species
SUBROUTINE readInitialGmsh2(sp, filename, density, velocity, temperature)
USE moduleRefParam
IMPLICIT NONE
INTEGER, INTENT(in):: sp
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
INTEGER:: i, e
INTEGER:: numNodes
OPEN(10, file = TRIM(filename))
!Skip first lines
DO i = 1, 11
READ(10, *)
END DO
!Reads number of nodes in file
READ(10, *) numNodes
ALLOCATE(density(1:numNodes))
ALLOCATE(velocity(1:numNodes, 1:3))
ALLOCATE(temperature(1:numNodes))
DO i = 1, numNodes
!Reads the density
READ(10, *), e, density(i)
END DO
DO i = 1, 10
READ(10, *)
END DO
DO i = 1, numNodes
!Reads the velocity
READ(10, *), e, velocity(i, 1:3)
END DO
!Skip uneccessary lines
DO i = 1, 10
READ(10, *)
END DO
!Assign density to nodes
DO i = 1, numNodes
!Skips pressure
READ(10, *)
END DO
!Skip uneccessary lines
DO i = 1, 10
READ(10, *)
END DO
!Assign density to nodes
DO i = 1, numNodes
!Skips pressure
READ(10, *) e, temperature(i)
END DO
END SUBROUTINE readInitialGmsh2
END MODULE moduleMeshInputGmsh2 END MODULE moduleMeshInputGmsh2

25
src/modules/mesh/moduleMesh.f90
View file

@ -161,6 +161,7 @@ MODULE moduleMesh
PROCEDURE(initVol_interface), DEFERRED, PASS:: init PROCEDURE(initVol_interface), DEFERRED, PASS:: init
PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes
PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter
PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF
PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK
@ -208,6 +209,12 @@ MODULE moduleMesh
END FUNCTION getNodesVol_interface END FUNCTION getNodesVol_interface
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION elemK_interface(self) RESULT(localK) PURE FUNCTION elemK_interface(self) RESULT(localK)
IMPORT:: meshVol IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self CLASS(meshVol), INTENT(in):: self
@ -271,9 +278,10 @@ MODULE moduleMesh
TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
!Array of volume elements !Array of volume elements
TYPE(meshVolCont), ALLOCATABLE:: vols(:) TYPE(meshVolCont), ALLOCATABLE:: vols(:)
PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL()
PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL()
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL()
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL()
CONTAINS CONTAINS
PROCEDURE, PASS:: doCollisions PROCEDURE, PASS:: doCollisions
@ -289,6 +297,15 @@ MODULE moduleMesh
END SUBROUTINE readMesh_interface END SUBROUTINE readMesh_interface
SUBROUTINE readInitial_interface(sp, filename, density, velocity, temperature)
INTEGER, INTENT(in):: sp
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
END SUBROUTINE readInitial_interface
!Connects volume and edges to the mesh !Connects volume and edges to the mesh
SUBROUTINE connectMesh_interface(self) SUBROUTINE connectMesh_interface(self)
IMPORT meshGeneric IMPORT meshGeneric
@ -655,8 +672,6 @@ MODULE moduleMesh
END SUBROUTINE doCollisions END SUBROUTINE doCollisions
SUBROUTINE doCoulomb(self) SUBROUTINE doCoulomb(self)
IMPORT meshParticles
CLASS(meshParticles), INTENT(inout):: self CLASS(meshParticles), INTENT(inout):: self
END SUBROUTINE doCoulomb END SUBROUTINE doCoulomb

108
src/modules/moduleInput.f90
View file

@ -206,7 +206,6 @@ MODULE moduleInput
CALL config%get(object // '.WeightingScheme', WSType, found) CALL config%get(object // '.WeightingScheme', WSType, found)
CALL solver%initWS(WSType) CALL solver%initWS(WSType)
!Makes tau(s) non-dimensional !Makes tau(s) non-dimensional
tau = tau / ti_ref tau = tau / ti_ref
tauMin = tauMin / ti_ref tauMin = tauMin / ti_ref
@ -219,7 +218,6 @@ MODULE moduleInput
!Read initial state for species !Read initial state for species
CALL verboseError('Reading Initial state...') CALL verboseError('Reading Initial state...')
CALL readInitial(config) CALL readInitial(config)
CALL checkStatus(config, "readInitial")
END SUBROUTINE readCase END SUBROUTINE readCase
@ -237,14 +235,21 @@ MODULE moduleInput
LOGICAL:: found LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: object CHARACTER(:), ALLOCATABLE:: object
INTEGER:: nInitial INTEGER:: nInitial
INTEGER:: i, p, e INTEGER:: i, j, p, e
CHARACTER(LEN=2):: iString CHARACTER(LEN=2):: iString
CHARACTER(:), ALLOCATABLE:: spName CHARACTER(:), ALLOCATABLE:: spName
INTEGER:: sp INTEGER:: sp
CHARACTER(:), ALLOCATABLE:: spFile CHARACTER(:), ALLOCATABLE:: spFile
INTEGER:: stat CHARACTER(:), ALLOCATABLE:: filename
CHARACTER(100):: dummy REAL(8), ALLOCATABLE, DIMENSION(:):: density, temperature
REAL(8):: density, velocity(1:3), temperature REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity
INTEGER, ALLOCATABLE, DIMENSION(:):: nodes
INTEGER:: nNodes
REAL(8), ALLOCATABLE, DIMENSION(:):: source, fPsi
!Density at the volume centroid
REAL(8):: densityCen
!Mean velocity and temperature at particle position
REAL(8):: velocityXi(1:3), temperatureXi
INTEGER:: nNewPart = 0.D0 INTEGER:: nNewPart = 0.D0
TYPE(particle), POINTER:: partNew TYPE(particle), POINTER:: partNew
REAL(8):: vTh REAL(8):: vTh
@ -260,38 +265,76 @@ MODULE moduleInput
CALL config%get(object // '.speciesName', spName, found) CALL config%get(object // '.speciesName', spName, found)
sp = speciesName2Index(spName) sp = speciesName2Index(spName)
CALL config%get(object // '.initialState', spFile, found) CALL config%get(object // '.initialState', spFile, found)
OPEN (10, FILE = path // spFile, ACTION = 'READ') !Reads node values at the nodes
DO filename = path // spFile
READ(10, '(A)', IOSTAT = stat) dummy CALL mesh%readInitial(sp, filename, density, velocity, temperature)
!If EoF, exit reading !For each volume in the node, create corresponding particles
IF (stat /= 0) EXIT DO e = 1, mesh%numVols
!If comment, skip
IF (INDEX(dummy,'#') /= 0) CYCLE
!Go up one line
BACKSPACE(10)
!Read information
READ(10, *) e, density, velocity, temperature
!Scale variables !Scale variables
!Particles in cell volume !Density at centroid of cell
nNewPart = INT(density * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) nodes = mesh%vols(e)%obj%getNodes()
!Non-dimensional velocity nNodes = SIZE(nodes)
velocity = velocity / v_ref !TODO: Procedure to obtain centroid from element (also for printing Electric Field)
!Non-dimensional temperature fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/))
temperature = temperature / T_ref ALLOCATE(source(1:nNodes))
!Non-dimensional thermal temperature DO j = 1, nNodes
vTh = DSQRT(temperature/species(sp)%obj%m) source(j) = density(nodes(j))
END DO
densityCen = DOT_PRODUCT(fPsi, source)
DEALLOCATE(fPsi)
!Calculate number of particles
nNewPart = INT(densityCen * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
!Allocate new particles !Allocate new particles
DO p = 1, nNewPart DO p = 1, nNewPart
ALLOCATE(partNew) ALLOCATE(partNew)
partNew%species => species(sp)%obj partNew%species => species(sp)%obj
partNew%v(1) = velocity(1) + vTh*randomMaxwellian()
partNew%v(2) = velocity(2) + vTh*randomMaxwellian()
partNew%v(3) = velocity(3) + vTh*randomMaxwellian()
partNew%vol = e
partNew%r = mesh%vols(e)%obj%randPos() partNew%r = mesh%vols(e)%obj%randPos()
partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r) partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r)
partNew%n_in = .TRUE. !Get mean velocity at particle position
partNew%weight = species(sp)%obj%weight fPsi = mesh%vols(e)%obj%fPsi(partNew%xi)
DO j = 1, nNodes
source(j) = velocity(nodes(j), 1)
END DO
velocityXi(1) = DOT_PRODUCT(fPsi, source)
DO j = 1, nNodes
source(j) = velocity(nodes(j), 2)
END DO
velocityXi(2) = DOT_PRODUCT(fPsi, source)
DO j = 1, nNodes
source(j) = velocity(nodes(j), 3)
END DO
velocityXi(3) = DOT_PRODUCT(fPsi, source)
velocityXi = velocityXi / v_ref
!Get temperature at particle position
DO j = 1, nNodes
source(j) = temperature(nodes(j))
END DO
temperatureXi = DOT_PRODUCT(fPsi, source)
temperatureXi = temperatureXi / T_ref
vTh = DSQRT(temperatureXi / species(sp)%obj%m)
partNew%v(1) = velocityXi(1) + vTh*randomMaxwellian()
partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian()
partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian()
partNew%vol = e
IF (ASSOCIATED(meshForMCC, mesh)) THEN
partNew%volColl = partNew%vol
ELSE
partNew%volColl = findCellBrute(meshColl, partNew%r)
END IF
partNew%n_in = .TRUE.
partNew%weight = species(sp)%obj%weight
!If charged species, add qm to particle !If charged species, add qm to particle
SELECT TYPE(sp => species(sp)%obj) SELECT TYPE(sp => species(sp)%obj)
TYPE IS (speciesCharged) TYPE IS (speciesCharged)
@ -305,8 +348,11 @@ MODULE moduleInput
!Assign particle to temporal list of particles !Assign particle to temporal list of particles
CALL partInitial%add(partNew) CALL partInitial%add(partNew)
END DO END DO
DEALLOCATE(source)
END DO END DO
END DO END DO