From cf4488976ae76b3d5a083c39e1f3114f2898f299 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 4 Apr 2026 20:03:20 +0200 Subject: [PATCH 1/3] Finally, the right values for velocity distributions --- src/modules/common/moduleRandom.f90 | 17 +-- src/modules/common/velocityDistribution.f90 | 3 +- src/modules/moduleInject.f90 | 114 ++++++-------------- 3 files changed, 44 insertions(+), 90 deletions(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index f4dcc46..a1bc3fc 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -67,19 +67,20 @@ MODULE moduleRandom end function randomMaxwellian !Returns a random number in a Maxwellian distribution of mean 0 and width 1 - FUNCTION randomHalfMaxwellian() RESULT(rnd) + function randomHalfMaxwellian() result(rnd) IMPLICIT NONE - REAL(8):: rnd - REAL(8):: x + real(8):: rnd + real(8):: v1 rnd = 0.D0 - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - END DO + v1 = 0.D0 + do while (v1 == 0.D0) + v1 = random() + + end do - rnd = DSQRT(-DLOG(x)) + rnd = sqrt(-2.d0*log(v1)) END FUNCTION randomHalfMaxwellian diff --git a/src/modules/common/velocityDistribution.f90 b/src/modules/common/velocityDistribution.f90 index 5b6c342..359eebd 100644 --- a/src/modules/common/velocityDistribution.f90 +++ b/src/modules/common/velocityDistribution.f90 @@ -53,12 +53,11 @@ module velocityDistribution function randomVelMaxwellian(self) result (v) use moduleRandom implicit none - class(velDistMaxwellian), intent(in):: self real(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian()/sqrt(2.d0) + v = self%vTh*randomMaxwellian() end function randomVelMaxwellian diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 39de1ad..7a1c5d6 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -12,7 +12,7 @@ MODULE moduleInject REAL(8):: temperature(1:3) !Temperature REAL(8):: n(1:3) !Direction of injection LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n - INTEGER:: nParticles !Number of particles to introduce each time step + INTEGER:: nParticles !Number of particles to inject each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges type(meshEdgePointer), allocatable:: edges(:) @@ -72,7 +72,6 @@ MODULE moduleInject USE moduleRefParam USE moduleConstParam USE moduleSpecies - USE moduleSolver USE moduleErrors use json_module IMPLICIT NONE @@ -211,8 +210,6 @@ MODULE moduleInject END DO - self%nParticles = SUM(self%particlesPerEdge) - ELSE ! No particles assigned per edge, use the species weight self%weightPerEdge = self%species%weight @@ -220,13 +217,14 @@ MODULE moduleInject self%particlesPerEdge(e) = max(1,FLOOR(fluxPerStep*self%edges(e)%obj%surface / self%species%weight)) END DO - self%nParticles = SUM(self%particlesPerEdge) - !Rescale weight to match flux self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles)) END IF + ! Total number of particles to inject + self%nParticles = SUM(self%particlesPerEdge) + !Scale particles for different species steps IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') @@ -314,126 +312,82 @@ MODULE moduleInject end subroutine updateQuasiNeutral - !Injection of particles - SUBROUTINE doInjects() - USE moduleSpecies - USE moduleSolver - IMPLICIT NONE - - INTEGER:: i - - !$OMP SINGLE - nPartInj = 0 - DO i = 1, nInject - IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN - nPartInj = nPartInj + injects(i)%obj%nParticles - - END IF - - END DO - - IF (ALLOCATED(partInj)) DEALLOCATE(partInj) - ALLOCATE(partInj(1:nPartInj)) - !$OMP END SINGLE - - DO i=1, nInject - IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN - CALL injects(i)%obj%addParticles() - - END IF - END DO - - END SUBROUTINE doInjects - !Add particles for the injection - SUBROUTINE addParticles(self) + SUBROUTINE addParticles(self, partArray) USE moduleSpecies - USE moduleSolver USE moduleMesh USE moduleRandom USE moduleErrors IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self - INTEGER, SAVE:: nMin + type(particle), allocatable, intent(inout):: partArray(:) INTEGER:: i, e INTEGER:: n, sp - CLASS(meshEdge), POINTER:: randomEdge + CLASS(meshEdge), POINTER:: edge REAL(8):: direction(1:3) + !$omp single + n = 0 + !$omp end single + !Insert particles - !$OMP SINGLE - nMin = 0 - DO i = 1, self%id -1 - IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN - nMin = nMin + injects(i)%obj%nParticles - - END IF - - END DO - nMin = nMin + 1 - !$OMP END SINGLE - !$OMP DO DO e = 1, self%nEdges ! Select edge for injection - randomEdge => self%edges(e)%obj + edge => self%edges(e)%obj ! Inject particles in edge DO i = 1, self%particlesPerEdge(e) - ! Index in the global partInj array - n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i - !Particle is considered to be outside the domain - partInj(n)%n_in = .FALSE. + ! Index in the partInj array + n = n + 1 + !Particle is considered to be outside the domain for now + partArray(n)%n_in = .FALSE. !Random position in edge - partInj(n)%r = randomEdge%randPos() + partArray(n)%r = edge%randPos() !Assign weight to particle. - partInj(n)%weight = self%weightPerEdge(e) + partArray(n)%weight = self%weightPerEdge(e) !Volume associated to the edge: - IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%cell = randomEdge%e1%n + IF (ASSOCIATED(edge%e1)) THEN + partArray(n)%cell = edge%e1%n - ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%cell = randomEdge%e2%n + ELSEIF (ASSOCIATED(edge%e2)) THEN + partArray(n)%cell = edge%e2%n ELSE CALL criticalError("No Volume associated to edge", 'addParticles') END IF - partInj(n)%cellColl = randomEdge%eColl%n + partArray(n)%cellColl = edge%eColl%n sp = self%species%n !Assign particle type - partInj(n)%species => self%species + partArray(n)%species => self%species if (all(self%n == 0.D0)) then - direction = randomEdge%normal + direction = edge%normal else direction = self%n end if - partInj(n)%v = 0.D0 + partArray(n)%v = 0.D0 - do while(dot_product(partInj(n)%v, direction) <= 0.d0) - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + ! do while(dot_product(partInj(n)%v, direction) <= 0.d0) + partArray(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) !If injecting a no-drift distribution and velocity is negative, reflect if ((self%vMod == 0.D0) .and. & - (dot_product(partInj(n)%v, direction) <= 0.D0)) then - partInj(n)%v = - partInj(n)%v + (dot_product(partArray(n)%v, direction) <= 0.D0)) then + partArray(n)%v = - partArray(n)%v end if - end do + ! end do !Obtain natural coordinates of particle in cell - partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) - !Push new particle with the minimum time step - CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) - !Assign cell to new particle - CALL solver%updateParticleCell(partInj(n)) + partArray(n)%Xi = mesh%cells(partArray(n)%cell)%obj%phy2log(partArray(n)%r) END DO @@ -500,7 +454,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistMaxwellian @@ -510,7 +464,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian From ed79eb018e695902ffadb84e3b0642b3dc684a85 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 4 Apr 2026 20:24:16 +0200 Subject: [PATCH 2/3] Import only what you need --- src/fpakc.f90 | 44 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 34 insertions(+), 10 deletions(-) diff --git a/src/fpakc.f90 b/src/fpakc.f90 index feee89a..cb68754 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -1,15 +1,39 @@ ! FPAKC main program PROGRAM fpakc - USE moduleCompTime - USE moduleCaseParam - USE moduleInput - USE moduleInject - USE moduleSolver - USE moduleMesh - USE moduleProbe - USE moduleErrors - USE OMP_LIB - IMPLICIT NONE + use moduleCompTime, only: tStep, & + tEMField, & + tCoul, & + tColl, & + tPush, & + tReset, & + tWeight + use omp_lib, only: omp_get_wtime + use moduleErrors, only: criticalError, & + verboseError + use moduleInput, only: readCOnfig, & + initOutput + use moduleCaseParam, only: timeStep, & + tInitial, & + tFinal + use moduleProbe, only: resetProbes + use moduleSolver, only: doScatter, & + doEMField, & + doOutput, & + solver, & + doPushes, & + doReset, & + doAverage, & + doInjects + use moduleMesh, only: boundariesEM_update, & + boundariesEM_update, & + boundariesParticle_update, & + mesh, & + meshForMCC, & + doMCCollisions, & + doCollisions, & + doCoulombScattering + use moduleInject, only: updateInjects + implicit none ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 From 6b72dbb108bdb8c01c9c250dbc433036d75bbbc0 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 4 Apr 2026 22:49:14 +0200 Subject: [PATCH 3/3] Organized injection --- src/modules/init/moduleInput.f90 | 7 ---- src/modules/makefile | 2 +- src/modules/moduleInject.f90 | 45 ++++++++++---------- src/modules/moduleList.f90 | 1 + src/modules/moduleSpecies.f90 | 3 -- src/modules/solver/moduleSolver.f90 | 65 ++++++++++++++++------------- 6 files changed, 63 insertions(+), 60 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index bcf4dcc..eb8bde5 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1431,7 +1431,6 @@ MODULE moduleInput CALL config%info('inject', found, n_children = nInject) ALLOCATE(injects(1:nInject)) - nPartInj = 0 DO i = 1, nInject WRITE(iString, '(i2)') i object = 'inject(' // trim(iString) // ')' @@ -1442,12 +1441,6 @@ MODULE moduleInput END DO - !Allocate array for injected particles - IF (nPartInj > 0) THEN - ALLOCATE(partInj(1:nPartInj)) - - END IF - END SUBROUTINE readInject SUBROUTINE readAverage(config) diff --git a/src/modules/makefile b/src/modules/makefile index 9f62c85..5789276 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -14,7 +14,7 @@ output.o: moduleSpecies.o common.o mesh.o: moduleList.o moduleSpecies.o moduleCollisions.o moduleCoulomb.o output.o common.o $(MAKE) -C mesh all -solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o +solver.o: moduleInject.o moduleSpecies.o moduleProbe.o common.o output.o mesh.o $(MAKE) -C solver all init.o: common.o solver.o moduleInject.o diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 7a1c5d6..942189c 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -313,23 +313,21 @@ MODULE moduleInject end subroutine updateQuasiNeutral !Add particles for the injection - SUBROUTINE addParticles(self, partArray) + SUBROUTINE addParticles(self) USE moduleSpecies USE moduleMesh USE moduleRandom USE moduleErrors + use moduleList, only: partInj IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self - type(particle), allocatable, intent(inout):: partArray(:) + type(particle), pointer:: part INTEGER:: i, e - INTEGER:: n, sp + INTEGER:: sp CLASS(meshEdge), POINTER:: edge - REAL(8):: direction(1:3) - !$omp single - n = 0 - !$omp end single + REAL(8):: direction(1:3) !Insert particles !$OMP DO @@ -338,30 +336,30 @@ MODULE moduleInject edge => self%edges(e)%obj ! Inject particles in edge DO i = 1, self%particlesPerEdge(e) + allocate(part) ! Index in the partInj array - n = n + 1 - !Particle is considered to be outside the domain for now - partArray(n)%n_in = .FALSE. + !Particle is considered to be inside the domain + part%n_in = .true. !Random position in edge - partArray(n)%r = edge%randPos() + part%r = edge%randPos() !Assign weight to particle. - partArray(n)%weight = self%weightPerEdge(e) + part%weight = self%weightPerEdge(e) !Volume associated to the edge: IF (ASSOCIATED(edge%e1)) THEN - partArray(n)%cell = edge%e1%n + part%cell = edge%e1%n ELSEIF (ASSOCIATED(edge%e2)) THEN - partArray(n)%cell = edge%e2%n + part%cell = edge%e2%n ELSE CALL criticalError("No Volume associated to edge", 'addParticles') END IF - partArray(n)%cellColl = edge%eColl%n + part%cellColl = edge%eColl%n sp = self%species%n !Assign particle type - partArray(n)%species => self%species + part%species => self%species if (all(self%n == 0.D0)) then direction = edge%normal @@ -371,23 +369,28 @@ MODULE moduleInject end if - partArray(n)%v = 0.D0 + part%v = 0.D0 ! do while(dot_product(partInj(n)%v, direction) <= 0.d0) - partArray(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + part%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) !If injecting a no-drift distribution and velocity is negative, reflect if ((self%vMod == 0.D0) .and. & - (dot_product(partArray(n)%v, direction) <= 0.D0)) then - partArray(n)%v = - partArray(n)%v + (dot_product(part%v, direction) <= 0.D0)) then + part%v = - part%v end if ! end do !Obtain natural coordinates of particle in cell - partArray(n)%Xi = mesh%cells(partArray(n)%cell)%obj%phy2log(partArray(n)%r) + part%Xi = mesh%cells(part%cell)%obj%phy2log(part%r) + + ! Add particle to global list + call partInj%setLock() + call partInj%add(part) + call partInj%unsetLock() END DO diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index b1dafdc..44524c6 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -23,6 +23,7 @@ MODULE moduleList END TYPE listNode + TYPE(listNode):: partInj !Particles comming from injections TYPE(listNode):: partWScheme !Particles comming from the nonAnalogue scheme TYPE(listNode):: partCollisions !Particles created in collisional process TYPE(listNode):: partSurfaces !Particles created in surface interactions diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index ab08f08..8c575da 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -48,11 +48,8 @@ MODULE moduleSpecies !Number of old particles INTEGER:: nPartOld - !Number of injected particles - INTEGER:: nPartInj !Arrays that contain the particles TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles CONTAINS FUNCTION speciesName2Index(speciesName) RESULT(sp) diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index f951839..e0ca06c 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -235,23 +235,21 @@ MODULE moduleSolver INTEGER:: nn, n, e INTEGER, SAVE:: nPartNew - INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces + INTEGER, SAVE:: nOldIn, nInj, nWScheme, nCollisions, nSurfaces TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:) INTEGER:: s !$OMP SECTIONS !$OMP SECTION - nInjIn = 0 - IF (ALLOCATED(partInj)) THEN - nInjIn = COUNT(partInj%n_in) - - END IF - !$OMP SECTION nOldIn = 0 IF (ALLOCATED(partOld)) THEN nOldIn = COUNT(partOld%n_in) END IF + + !$OMP SECTION + nInj = partInj%amount + !$OMP SECTION nWScheme = partWScheme%amount @@ -267,30 +265,14 @@ MODULE moduleSolver !$OMP SINGLE CALL MOVE_ALLOC(partOld, partTemp) - nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces + nPartNew = nInj + nOldIn + nWScheme + nCollisions + nSurfaces ALLOCATE(partOld(1:nPartNew)) !$OMP END SINGLE !$OMP SECTIONS - !$OMP SECTION - !Reset particles from injection - nn = 0 - DO n = 1, nPartInj - IF (partInj(n)%n_in) THEN - nn = nn + 1 - partOld(nn) = partInj(n) - IF (nProbes > 0) THEN - CALL doProbes(partOld(nn)) - - END IF - - END IF - - END DO - !$OMP SECTION !Reset particles from previous iteration - nn = nInjIn + nn = 0 DO n = 1, nPartOld IF (partTemp(n)%n_in) THEN nn = nn + 1 @@ -304,19 +286,24 @@ MODULE moduleSolver END DO + !$OMP SECTION + !Reset particles from injection + nn = nOldIn + call resetList(partInj, partOld, nn) + !$OMP SECTION !Reset particles from weighting scheme - nn = nInjIn + nOldIn + nn = nOldIn + nInj CALL resetList(partWScheme, partOld, nn) !$OMP SECTION !Reset particles from collisional process - nn = nInjIn + nOldIn + nWScheme + nn = nOldIn + nInj + nWScheme CALL resetList(partCollisions, partOld, nn) !$OMP SECTION !Reset particles from surface process - nn = nInjIn + nOldIn + nWScheme + nCollisions + nn = nOldIn + nInj + nWScheme + nCollisions CALL resetList(partSurfaces, partOld, nn) !$OMP SECTION @@ -473,6 +460,28 @@ MODULE moduleSolver END SUBROUTINE splitParticle + !Injection of particles + SUBROUTINE doInjects() + USE moduleSpecies + use moduleInject + IMPLICIT NONE + + INTEGER:: i, n, sp + + DO i=1, nInject + associate(inject => injects(i)%obj) + sp = inject%species%n + IF (solver%pusher(sp)%pushSpecies) THEN + + CALL inject%addParticles() + + END IF + end associate + + END DO + + END SUBROUTINE doInjects + SUBROUTINE updateParticleCell(self, part) USE moduleSpecies USE moduleMesh