diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 99e12e7..6b3b20a 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index 6477cff..479efb2 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -29,8 +29,7 @@ \makeglossaries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -\newacronym{fpakc}{fpakc}{Finite Element PArticle Code} +\newacronym{fpakc}{fpakc}{Finite element PArticle Code} \newacronym{mpi}{MPI}{Message Passing Interface} \newacronym{gpu}{GPU}{Graphics Processing Unit} \newacronym{cpu}{CPU}{Central Processing Unit} @@ -308,14 +307,13 @@ make \end{lstlisting} to compile the code. If everything is correct, an executable named \textit{fpakc} will be generated. - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Running the code} To run a case, simply execute: \begin{lstlisting} ./fpakc path/to/input-file.json \end{lstlisting} - in a command line from the root \acrshort{fpakc} folder. + in a command line from the root \acrshort{fpakc} folder. Substitute \lstinline|path/to/input-file.json| with the path to the input file of the case you want to run. The examples in the run directory are presented in Chapter \ref{ch:exampleRuns}. @@ -659,6 +657,10 @@ make The file needs to be located in the folder \textbf{output.folder}. If this value is not present, the mesh defined in \textbf{geometry.meshFile} is used for \acrshort{mcc}. The format of this mesh needs to be the same as the one defined in \textbf{geometry.meshType}. + \item \textbf{timeStep}: Real. + Units of $\unit{s}$. + Time step to calculate MCC. + If no time is provided, the minimum time step is used. \item \textbf{collisions}: Object. Array. Contains the different short range interactions (\acrshort{mcc}). diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 308a0b5..8bd909b 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -70,7 +70,7 @@ PROGRAM fpakc tColl = omp_get_wtime() !$OMP END SINGLE - IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions() + IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions(t) !$OMP SINGLE tColl = omp_get_wtime() - tColl diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index db52de4..13e25ff 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -42,7 +42,7 @@ MODULE moduleMeshOutput0D USE moduleOutput IMPLICIT NONE - CLASS(meshGeneric), INTENT(in):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER, INTENT(in):: t CHARACTER(:), ALLOCATABLE:: fileName @@ -56,7 +56,7 @@ MODULE moduleMeshOutput0D END IF OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE(20, "(ES20.6E3, I20)") REAL(t)*tauMin*ti_ref, mesh%vols(1)%obj%nColl + WRITE(20, "(ES20.6E3, I20)") REAL(t)*tauMin*ti_ref, self%vols(1)%obj%nColl CLOSE(20) END SUBROUTINE printColl0D diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 4159873..66fe345 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -95,7 +95,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleOutput IMPLICIT NONE - CLASS(meshGeneric), INTENT(in):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER, INTENT(in):: t INTEGER:: numEdges INTEGER:: n diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 86343f0..77170da 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -310,7 +310,7 @@ MODULE moduleMesh SUBROUTINE printColl_interface(self, t) IMPORT meshGeneric - CLASS(meshGeneric), INTENT(in):: self + CLASS(meshGeneric), INTENT(inout):: self INTEGER, INTENT(in):: t END SUBROUTINE printColl_interface @@ -637,7 +637,7 @@ MODULE moduleMesh END FUNCTION findCellBrute !Computes collisions in element - SUBROUTINE doCollisions(self) + SUBROUTINE doCollisions(self, t) USE moduleCollisions USE moduleSpecies USE moduleList @@ -646,6 +646,7 @@ MODULE moduleMesh IMPLICIT NONE CLASS(meshGeneric), INTENT(inout), TARGET:: self + INTEGER, INTENT(in):: t INTEGER:: e CLASS(meshVol), POINTER:: vol INTEGER:: nPart !Number of particles inside the cell @@ -657,49 +658,57 @@ MODULE moduleMesh REAL(8):: sigmaVrelMaxNew TYPE(pointerArray), ALLOCATABLE:: partTemp(:) - !$OMP DO SCHEDULE(DYNAMIC) - DO e=1, self%numVols - vol => self%vols(e)%obj - nPart = vol%listPart_in%amount - !Calculates number of collisions if there is more than one particle in the cell - IF (nPart > 1) THEN - !Probability of collision - pMax = vol%totalWeight*vol%sigmaVrelMax*tauMin/vol%volume + IF (MOD(t, everyColl) == 0) THEN + !Collisions need to be performed in this iteration + !$OMP DO SCHEDULE(DYNAMIC) + DO e=1, self%numVols + vol => self%vols(e)%obj + nPart = vol%listPart_in%amount - !Number of collisions in the cell - vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) + !Resets the number of collisions + vol%nColl = 0 - IF (vol%nColl > 0) THEN - !Converts the list of particles to an array for easy access - partTemp = vol%listPart_in%convert2Array() + !Calculates number of collisions if there is more than one particle in the cell + IF (nPart > 1) THEN + !Probability of collision + pMax = vol%totalWeight*vol%sigmaVrelMax*tauColl/vol%volume - END IF + !Number of collisions in the cell + vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) - DO n = 1, vol%nColl - !Select random numbers - rnd = random(1, nPart) - part_i => partTemp(rnd)%part - rnd = random(1, nPart) - part_j => partTemp(rnd)%part - ij = interactionIndex(part_i%species%n, part_j%species%n) - sigmaVrelMaxNew = 0.D0 - DO k = 1, interactionMatrix(ij)%amount - CALL interactionMatrix(ij)%collisions(k)%obj%collide(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) - - END DO - - !Update maximum cross section*v_rel per each collision - IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN - vol%sigmaVrelMax = sigmaVrelMaxNew + IF (vol%nColl > 0) THEN + !Converts the list of particles to an array for easy access + partTemp = vol%listPart_in%convert2Array() END IF - END DO + DO n = 1, vol%nColl + !Select random numbers + rnd = random(1, nPart) + part_i => partTemp(rnd)%part + rnd = random(1, nPart) + part_j => partTemp(rnd)%part + ij = interactionIndex(part_i%species%n, part_j%species%n) + sigmaVrelMaxNew = 0.D0 + DO k = 1, interactionMatrix(ij)%amount + CALL interactionMatrix(ij)%collisions(k)%obj%collide(vol%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) - END IF + END DO - END DO - !$OMP END DO + !Update maximum cross section*v_rel per each collision + IF (sigmaVrelMaxNew > vol%sigmaVrelMax) THEN + vol%sigmaVrelMax = sigmaVrelMaxNew + + END IF + + END DO + + END IF + + END DO + !$OMP END DO + + END IF END SUBROUTINE doCollisions diff --git a/src/modules/moduleCaseParam.f90 b/src/modules/moduleCaseParam.f90 index c8df0e4..8df3210 100644 --- a/src/modules/moduleCaseParam.f90 +++ b/src/modules/moduleCaseParam.f90 @@ -4,6 +4,7 @@ MODULE moduleCaseParam INTEGER:: tFinal, tInitial = 0 REAL(8), ALLOCATABLE:: tau(:) REAL(8):: tauMin + REAL(8):: tauColl END MODULE moduleCaseParam diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 6557fa7..2f22204 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -2,6 +2,9 @@ MODULE moduleCollisions USE moduleSpecies USE moduleTable + !Integer for when collisions are computed + INTEGER:: everyColl + !Abstract type for collision between two particles TYPE, ABSTRACT:: collisionBinary REAL(8):: rMass !Reduced mass diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 34e7507..19255dd 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -40,6 +40,11 @@ MODULE moduleInput CALL readSpecies(config) CALL checkStatus(config, "readSpecies") + !Reads case parameters + CALL verboseError('Reading Case parameters...') + CALL readCase(config) + CALL checkStatus(config, "readCase") + !Read interactions between species CALL verboseError('Reading interaction between species...') CALL readInteractions(config) @@ -55,10 +60,10 @@ MODULE moduleInput CALL readGeometry(config) CALL checkStatus(config, "readGeometry") - !Reads case parameters - CALL verboseError('Reading Case parameters...') - CALL readCase(config) - CALL checkStatus(config, "readCase") + !Read initial state for species + CALL verboseError('Reading Initial state...') + CALL readInitial(config) + CALL checkStatus(config, "readInitial") !Read injection of particles CALL verboseError('Reading injection of particles from boundaries...') @@ -233,11 +238,6 @@ MODULE moduleInput WRITE(tString, '(I1)') iterationDigits iterationFormat = "(I" // tString // "." // tString // ")" - !Read initial state for species - CALL verboseError('Reading Initial state...') - CALL readInitial(config) - CALL checkStatus(config, "readInitial") - END SUBROUTINE readCase !Reads the initial information for the species @@ -558,6 +558,8 @@ MODULE moduleInput USE moduleCollisions USE moduleErrors USE moduleMesh + USE moduleCaseParam + USE moduleRefParam USE OMP_LIB USE json_module IMPLICIT NONE @@ -595,6 +597,24 @@ MODULE moduleInput END IF + !Reads collision time step + CALL config%info('interactions.timeStep', found) + IF (found) THEN + CALL config%get('interactions.timeStep', tauColl, found) + tauColl = tauColl / ti_ref + + ELSE + tauColl = tauMin + + END IF + + IF (tauColl < tauMin) THEN + CALL criticalError('Collisional time step below minimum time step.', 'readInteractions') + + END IF + + everyColl = NINT(tauColl / tauMin) + !Inits the MCC matrix CALL initInteractionMatrix(interactionMatrix)