From 5dfc8d4ce9623b9929a145ce8a6fbecb41e89325 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 20 Feb 2026 11:24:58 +0100 Subject: [PATCH] Inject fixed. Now getting close to input --- src/modules/init/moduleInput.f90 | 58 +++++++++++++++++++++++++++++--- src/modules/moduleInject.f90 | 55 ++++++++---------------------- 2 files changed, 67 insertions(+), 46 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 4a7930b..30bb3aa 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1128,12 +1128,35 @@ MODULE moduleInput !Builds the K matrix for the Particles mesh CALL mesh%constructGlobalK() - !Assign the procedure to find a volume for meshColl + !Assign the procedure to find a cell for meshColl IF (doubleMesh) THEN findCellColl => findCellCollMesh + ! Link edges with cells in meshColl + DO e=1, mesh%numEdges + nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos()) + IF (nVolColl > 0) THEN + mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj + + ELSE + CALL criticalError("No connection between edge and meshColl", "initInject") + + END IF + + end do ELSE findCellColl => findCellSameMesh + ! Link edges with cells in meshColl + DO e=1, mesh%numEdges + IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN + mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1 + + ELSE + mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2 + + END IF + + end do END IF @@ -1142,10 +1165,35 @@ MODULE moduleInput do b = 1, nBoundariesParticle select type(bound => boundariesParticle(b)%obj) type is(boundaryQuasiNeutrality) - ! Loop over all physical surfaces - do ps = 1, nPhysicalSurfaces - do s = 1, nSpecies - if (associated(physicalSurfaces(ps)%particles(s), bound)) then + ! Loop over all physical surfaces + do ps = 1, nPhysicalSurfaces + ! Loop over all species + do s = 1, nSpecies + ! If the boundary for the species is linked to the one analysing, add the edges + if (associated(physicalSurfaces(ps)%particles(s), bound)) then + bound%edges = [bound%edges, physicalSurfaces(ps)%edges] + end if + + end do + + end do + + end select + + end do + + ! EM Boundaries + do b = 1, nBoundariesEM + select type(bound => boundariesEM(b)%obj) + type is(boundaryEMDirichlet) + ! Loop over all physical surfaces + do ps = 1, nPhysicalSurfaces + ! If the boundary for the species is linked to the one analysing, add the edges + if (associated(physicalSurfaces(ps)%EM, bound)) then + bound%nodes = [bound%nodes, physicalSurfaces(ps)%nodes] + end if + + end do end select diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 17013ae..44c0f80 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -2,6 +2,7 @@ MODULE moduleInject USE moduleSpecies use velocityDistribution + use moduleMesh, only: meshEdgePointer !Generic injection of particles TYPE:: injectGeneric @@ -14,7 +15,7 @@ MODULE moduleInject INTEGER:: nParticles !Number of particles to introduce each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges - INTEGER, ALLOCATABLE:: edges(:) !Array with edges + type(meshEdgePointer), allocatable:: edges(:) INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge REAL(8):: surface ! Total surface of injection @@ -43,6 +44,7 @@ MODULE moduleInject INTEGER, INTENT(in):: i REAL(8), INTENT(in):: v, n(1:3), temperature(1:3) INTEGER, INTENT(in):: sp, ps, particlesPerEdge + integer:: ps_index REAL(8):: tauInject REAL(8), INTENT(in):: flow CHARACTER(:), ALLOCATABLE, INTENT(in):: units @@ -55,50 +57,21 @@ MODULE moduleInject self%vMod = v / v_ref self%n = n / NORM2(n) self%temperature = temperature / T_ref - !Gets the edge elements from which particles are injected - DO e = 1, mesh%numEdges - phSurface(e) = mesh%edges(e)%obj%physicalSurface - END DO - self%nEdges = COUNT(phSurface == physicalSurface) - ALLOCATE(self%edges(1:self%nEdges)) - ALLOCATE(self%particlesPerEdge(1:self%nEdges)) - ALLOCATE(self%weightPerEdge(1:self%nEdges)) - et = 0 - DO e=1, mesh%numEdges - IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN - et = et + 1 - self%edges(et) = mesh%edges(e)%obj%n - !Assign connectivity between injection edge and meshColl volume - IF (doubleMesh) THEN - nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos()) - IF (nVolColl > 0) THEN - mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj + !Get array index of physical surface for inject + ps_index = physicalSurface_to_index(ps) - ELSE - CALL criticalError("No connection between edge and meshColl", "initInject") + ! Copy array of edges + self%edges = physicalSurfaces(ps_index)%edges + self%nEdges = size(self%edges) - END IF - - ELSE - IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN - mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1 - - ELSE - mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2 - - END IF - - END IF - - END IF - - END DO + allocate(self%particlesPerEdge(1:self%nEdges)) + allocate(self%weightPerEdge(1:self%nEdges)) !Calculates total area self%surface = 0.D0 DO et = 1, self%nEdges - self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface + self%surface = self%surface + self%edges(et)%obj%surface END DO @@ -149,7 +122,7 @@ MODULE moduleInject ! Particles per edge defined by the user self%particlesPerEdge = particlesPerEdge DO et = 1, self%nEdges - self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge) + self%weightPerEdge(et) = fluxPerStep * self%edges(et)%obj%surface / real(particlesPerEdge) END DO @@ -159,7 +132,7 @@ MODULE moduleInject ! No particles assigned per edge, use the species weight self%weightPerEdge = self%species%weight DO et = 1, self%nEdges - self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight)) + self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*self%edges(et)%obj%surface / self%species%weight)) END DO self%nParticles = SUM(self%particlesPerEdge) @@ -266,7 +239,7 @@ MODULE moduleInject !$OMP DO DO e = 1, self%nEdges ! Select edge for injection - randomEdge => mesh%edges(self%edges(e))%obj + randomEdge => self%edges(e)%obj ! Inject particles in edge DO i = 1, self%particlesPerEdge(e) ! Index in the global partInj array