From ecb1364d6ad19b7b2dc6da797656cc40d26d01e6 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 6 Mar 2026 20:00:46 +0100 Subject: [PATCH] Boundaries (EM and particles) can now be updated at every iteration --- src/fpakc.f90 | 9 +++ src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 10 ++-- src/modules/mesh/moduleMesh.f90 | 39 +++++++++++-- src/modules/mesh/moduleMesh@boundaryEM.f90 | 35 ++++++++++-- .../mesh/moduleMesh@boundaryParticle.f90 | 56 ++++++++++++------- 5 files changed, 115 insertions(+), 34 deletions(-) diff --git a/src/fpakc.f90 b/src/fpakc.f90 index ae6e7bb..18e27f0 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -45,6 +45,9 @@ PROGRAM fpakc !$OMP SINGLE CALL verboseError("Calculating initial EM field...") + ! Update EM boundary models + call boundariesEM_update() + !$OMP END SINGLE CALL doEMField() !$OMP END PARALLEL @@ -62,6 +65,12 @@ PROGRAM fpakc ! Update global time step index timeStep = t + ! Update Particle boundary models + call boundariesParticle_update + + ! Update EM boundary models + call boundariesEM_update() + !Checks if a species needs to me moved in this iteration CALL solver%updatePushSpecies() diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index e7cc060..a68d1f7 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -25,11 +25,11 @@ MODULE moduleMesh3DCart CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() CONTAINS !meshEdge DEFERRED PROCEDURES - PROCEDURE, PASS:: init => initEdge3DCartTria - PROCEDURE, PASS:: getNodes => getNodes3DCartTria - PROCEDURE, PASS:: intersection => intersection3DCartTria - PROCEDURE, PASS:: randPos => randPosEdgeTria - procedure, pass:: center => centerEdgeTria + PROCEDURE, PASS:: init => initEdge3DCartTria + PROCEDURE, PASS:: getNodes => getNodes3DCartTria + PROCEDURE, PASS:: intersection => intersection3DCartTria + PROCEDURE, PASS:: randPos => randPosEdgeTria + procedure, pass:: center => centerEdgeTria !PARTICULAR PROCEDURES PROCEDURE, NOPASS:: fPsi => fPsiTria diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 179097f..9133cca 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -257,8 +257,8 @@ MODULE moduleMesh PROCEDURE(detJac_interface), DEFERRED, NOPASS:: detJac PROCEDURE(invJac_interface), DEFERRED, NOPASS:: invJac !Procedures to get specific values in the node - PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherElectricField - PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherMagneticField + PROCEDURE(gatherCellArray_interface), DEFERRED, PASS:: gatherElectricField + PROCEDURE(gatherCellArray_interface), DEFERRED, PASS:: gatherMagneticField !Compute K and F to solve PDE on the mesh PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF @@ -390,13 +390,13 @@ MODULE moduleMesh END FUNCTION invJac_interface - PURE FUNCTION gatherArray_interface(self, Xi) RESULT(array) + PURE FUNCTION gatherCellArray_interface(self, Xi) RESULT(array) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) REAL(8):: array(1:3) - END FUNCTION gatherArray_interface + END FUNCTION gatherCellArray_interface PURE FUNCTION elemK_interface(self, nNodes) RESULT(localK) IMPORT:: meshCell @@ -650,6 +650,7 @@ MODULE moduleMesh type, abstract, public:: boundaryParticleGeneric integer:: n character(:), allocatable:: name + procedure(updateParticle_interface), pointer, pass:: update => null() contains procedure(applyParticle_interface), deferred, pass:: apply @@ -698,6 +699,14 @@ MODULE moduleMesh end subroutine + ! Update the values of the particle boundary model + subroutine updateParticle_interface(self) + import boundaryParticleGeneric + + class(boundaryParticleGeneric), intent(inout):: self + + end subroutine updateParticle_interface + end interface !Reflecting boundary @@ -754,6 +763,7 @@ MODULE moduleMesh ! Ensures quasi-neutrality by changing the reflection coefficient type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality real(8):: alpha ! Reflection parameter + integer:: s_incident ! species index of the incident species type(meshEdgePointer), allocatable:: edges(:) !Array with edges contains procedure, pass:: apply => quasiNeutrality @@ -867,6 +877,14 @@ MODULE moduleMesh !Array for boundaries type(boundaryParticleCont), allocatable, target:: boundariesParticle(:) + ! Update the particle boundary models + interface + module subroutine boundariesParticle_update() + + end subroutine boundariesParticle_update + + end interface + ! BOUNDARY ELECTROMAGNETIC DEFINITIONS ! Generic type for electromagnetic boundary conditions type, public, abstract:: boundaryEMGeneric @@ -918,11 +936,11 @@ MODULE moduleMesh end subroutine applyEM_interface - ! Update the values of the boundary condition + ! Update the values of the EM boundary model subroutine updateEM_interface(self) import boundaryEMGeneric - class(boundaryEMGeneric), intent(in):: self + class(boundaryEMGeneric), intent(inout):: self end subroutine updateEM_interface @@ -939,6 +957,7 @@ MODULE moduleMesh TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime real(8):: potential + real(8):: timeFactor type(table1D):: temporalProfile contains ! boundaryEMGeneric DEFERRED PROCEDURES @@ -973,6 +992,14 @@ MODULE moduleMesh !Information of charge and reference parameters for rho vector REAL(8), ALLOCATABLE:: qSpecies(:) + ! Update the EM boundary models + interface + module subroutine boundariesEM_update() + + end subroutine boundariesEM_update + + end interface + ! PHYSICAL SURFACES LINKING TO MESH ELEMENTS ! Link physical surface to edges type:: physicalSurface diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index 95c8fb3..14fd0ee 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -89,6 +89,8 @@ submodule(moduleMesh) boundaryEM call self%temporalProfile%convert(1.D0/ti_ref, 1.D0) + self%update => updateDirichletTime + end select END SUBROUTINE initDirichletTime @@ -118,17 +120,42 @@ submodule(moduleMesh) boundaryEM CLASS(boundaryEMDirichletTime), INTENT(in):: self REAL(8), INTENT(inout):: vectorF(:) - REAL(8):: timeFactor integer:: n - timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin) - DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential * timeFactor + self%nodes(n)%obj%emData%phi = self%potential * self%timeFactor vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi END DO END SUBROUTINE applyDirichletTime + module subroutine updateDirichletTime(self) + implicit none + + class(boundaryEMGeneric), intent(inout):: self + + select type(self) + type is(boundaryEMDirichletTime) + self%timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin) + + end select + + end subroutine updateDirichletTime + + module subroutine boundariesEM_update() + implicit none + + integer:: b + + do b = 1, nBoundariesEM + if (associated(boundariesEM(b)%obj%update)) then + call boundariesEM(b)%obj%update + + end if + + end do + + end subroutine boundariesEM_update + end submodule boundaryEM diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index de4edee..040f633 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -107,6 +107,8 @@ submodule(moduleMesh) boundaryParticle boundary%alpha = 0.d0 allocate(boundary%edges(0)) + boundary%update => quasiNeutrality_update + end select end subroutine initQuasiNeutrality @@ -328,24 +330,10 @@ submodule(moduleMesh) boundaryParticle implicit none class(boundaryQuasiNeutrality), intent(in):: self - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - real(8), allocatable:: density(:) - class(meshCell), pointer:: cell - real(8):: EF_dir - real(8):: alpha + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part - alpha = 0.85d0 - - if (associated(edge%e1)) then - cell => edge%e1 - - else - cell => edge%e2 - - end if - - if (random() <= alpha) then + if (random() <= self%alpha) then call genericReflection(edge, part) else @@ -355,14 +343,29 @@ submodule(moduleMesh) boundaryParticle end subroutine quasiNeutrality + subroutine quasiNeutrality_update(self) + implicit none + + class(boundaryParticleGeneric), intent(inout):: self + + select type(self) + type is(boundaryQuasiNeutrality) + self%alpha = 0.1d0 + + print *, self%alpha + + end select + + end subroutine quasiNeutrality_update + ! Generic boundary conditions for internal use module subroutine genericReflection(edge, part) use moduleCaseParam use moduleSpecies implicit none - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part !rp = intersection between particle and edge !rpp = final position of particle !vpp = final velocity of particle @@ -393,4 +396,19 @@ submodule(moduleMesh) boundaryParticle end subroutine genericTransparent + module subroutine boundariesParticle_update() + implicit none + + integer:: b + + do b = 1, nBoundariesParticle + if (associated(boundariesParticle(b)%obj%update)) then + call boundariesParticle(b)%obj%update + + end if + + end do + + end subroutine boundariesParticle_update + end submodule boundaryParticle