From c1b6cf1b319d8aa60fdb33420e893fa0e4793506 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 10 Feb 2026 17:22:56 +0100 Subject: [PATCH] Now the deferred apply is pass and I made some generic boundaries to use internally --- src/modules/mesh/moduleMesh.f90 | 165 ++++++++------ src/modules/mesh/moduleMesh@boundary.f90 | 270 +++++++++++------------ 2 files changed, 216 insertions(+), 219 deletions(-) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 3295b57..a951eda 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -589,8 +589,8 @@ MODULE moduleMesh TYPE, abstract, PUBLIC:: boundaryGeneric character(:), allocatable:: name contains - procedure, pass:: init => initBoundary - procedure(boundary_interface), deferred, nopass:: apply + procedure, pass:: init => initBoundary + procedure(boundary_interface), deferred, pass:: apply END TYPE boundaryGeneric @@ -607,13 +607,13 @@ MODULE moduleMesh end interface abstract interface - subroutine boundary_interface(edge, part, self) + subroutine boundary_interface(self, edge, part) use moduleSpecies import boundaryGeneric, meshEdge - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self + class(boundaryGeneric), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part end subroutine @@ -622,28 +622,28 @@ MODULE moduleMesh !Reflecting boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryReflection CONTAINS - procedure, nopass:: apply => reflection + procedure, pass:: apply => reflection END TYPE boundaryReflection !Absorption boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAbsorption CONTAINS - procedure, nopass:: apply => absorption + procedure, pass:: apply => absorption END TYPE boundaryAbsorption !Transparent boundary TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryTransparent CONTAINS - procedure, nopass:: apply => transparent + procedure, pass:: apply => transparent END TYPE boundaryTransparent !Symmetry axis TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis CONTAINS - procedure, nopass:: apply => axis + procedure, pass:: apply => axis END TYPE boundaryAxis @@ -652,7 +652,7 @@ MODULE moduleMesh !Thermal velocity of the wall: square root(Wall temperature X specific heat) REAL(8):: vTh CONTAINS - procedure, nopass:: apply => wallTemperature + procedure, pass:: apply => wallTemperature END TYPE boundaryWallTemperature @@ -666,7 +666,7 @@ MODULE moduleMesh REAL(8):: eThreshold REAL(8):: deltaV CONTAINS - procedure, nopass:: apply => ionization + procedure, pass:: apply => ionization END TYPE boundaryIonization @@ -675,80 +675,99 @@ MODULE moduleMesh real(8):: alpha ! Reflection parameter integer, allocatable:: edges(:) !Array with edges contains - procedure, nopass:: apply => quasiNeutrality + procedure, pass:: apply => quasiNeutrality end type boundaryQuasiNeutrality !Wrapper for boundary types (one per species) interface - module subroutine pointBoundaryFunction(edge, s) + module subroutine reflection(self, edge, part) + use moduleSpecies + + class(boundaryReflection), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine reflection + + module subroutine absorption(self, edge, part) + use moduleSpecies + + class(boundaryAbsorption), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine absorption + + module subroutine transparent(self, edge, part) + use moduleSpecies + + class(boundaryTransparent), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine transparent + + module subroutine axis(self, edge, part) + use moduleSpecies + + class(boundaryAxis), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine axis + + module subroutine wallTemperature(self, edge, part) + use moduleSpecies + + class(boundaryWallTemperature), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine wallTemperature + + module subroutine ionization(self, edge, part) + use moduleSpecies + + class(boundaryIonization), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine ionization + + module subroutine quasiNeutrality(self, edge, part) + use moduleSpecies + + class(boundaryQuasiNeutrality), intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + end subroutine quasiNeutrality + + ! Generic basic boundary conditions to use internally in the code + module subroutine genericReflection(edge, part) + use moduleSpecies + class(meshEdge), intent(inout):: edge - integer, intent(in):: s !Species index + class(particle), intent(inout):: part - end subroutine pointBoundaryFunction + end subroutine genericReflection - module subroutine reflection(edge, part, self) - use moduleSpecies + module subroutine genericAbsorption(edge, part) + use moduleSpecies - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part - end subroutine reflection + end subroutine genericAbsorption - module subroutine absorption(edge, part, self) - use moduleSpecies + module subroutine genericTransparent(edge, part) + use moduleSpecies - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part - end subroutine absorption - - module subroutine transparent(edge, part, self) - use moduleSpecies - - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self - - end subroutine transparent - - module subroutine axis(edge, part, self) - use moduleSpecies - - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self - - end subroutine axis - - module subroutine wallTemperature(edge, part, self) - use moduleSpecies - - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self - - end subroutine wallTemperature - - module subroutine ionization(edge, part, self) - use moduleSpecies - - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self - - end subroutine ionization - - module subroutine quasiNeutrality(edge, part, self) - use moduleSpecies - - class(meshEdge), intent(inout):: edge - class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self - - end subroutine quasiNeutrality + end subroutine genericTransparent end interface diff --git a/src/modules/mesh/moduleMesh@boundary.f90 b/src/modules/mesh/moduleMesh@boundary.f90 index 5603d28..e3805a5 100644 --- a/src/modules/mesh/moduleMesh@boundary.f90 +++ b/src/modules/mesh/moduleMesh@boundary.f90 @@ -1,21 +1,6 @@ !moduleMeshBoundary: Boundary functions for the mesh edges submodule(moduleMesh) boundary - CONTAINS - ! FUNCTION getBoundaryId(physicalSurface) RESULT(id) - ! IMPLICIT NONE - - ! INTEGER:: physicalSurface - ! INTEGER:: id - ! INTEGER:: i - - ! id = 0 - ! DO i = 1, nBoundary - ! IF (physicalSurface == boundariesParticle(i)%obj%physicalSurface) id = boundariesParticle(i)%obj%n - - ! END DO - - ! END FUNCTION getBoundaryId - + contains module subroutine initBoundary(self, config, object) use json_module use moduleRefParam, only: m_ref @@ -185,14 +170,14 @@ submodule(moduleMesh) boundary end subroutine initQuasiNeutrality - module SUBROUTINE reflection(edge, part, self) + module SUBROUTINE reflection(self, edge, part) USE moduleCaseParam USE moduleSpecies IMPLICIT NONE + class(boundaryReflection), intent(in):: self CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - class(boundaryGeneric), optional, intent(in):: self !rp = intersection between particle and edge !rpp = final position of particle !vpp = final velocity of particle @@ -212,14 +197,14 @@ submodule(moduleMesh) boundary END SUBROUTINE reflection !Absoption in a surface - SUBROUTINE absorption(edge, part, self) + SUBROUTINE absorption(self, edge, part) USE moduleCaseParam USE moduleSpecies IMPLICIT NONE + class(boundaryAbsorption), intent(in):: self CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - class(boundaryGeneric), optional, intent(in):: self REAL(8):: rpp(1:3) !Position of particle projected to the edge REAL(8):: d !Distance from particle to edge @@ -249,13 +234,13 @@ submodule(moduleMesh) boundary END SUBROUTINE absorption !Transparent boundary condition - SUBROUTINE transparent(edge, part, self) + SUBROUTINE transparent(self, edge, part) USE moduleSpecies IMPLICIT NONE + class(boundaryTransparent), intent(in):: self CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - class(boundaryGeneric), optional, intent(in):: self !Removes particle from domain part%n_in = .FALSE. @@ -265,46 +250,42 @@ submodule(moduleMesh) boundary !Symmetry axis. Reflects particles. !Although this function should never be called, it is set as a reflective boundary !to properly deal with possible particles reaching a corner and selecting this boundary. - SUBROUTINE symmetryAxis(edge, part, self) + SUBROUTINE symmetryAxis(self, edge, part) USE moduleSpecies IMPLICIT NONE + class(boundaryAxis), intent(in):: self CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - class(boundaryGeneric), optional, intent(in):: self - CALL reflection(edge, part) + CALL genericReflection(edge, part) END SUBROUTINE symmetryAxis !Wall with temperature - SUBROUTINE wallTemperature(edge, part, self) + SUBROUTINE wallTemperature(self, edge, part) USE moduleSpecies USE moduleRandom IMPLICIT NONE + class(boundaryWallTemperature), intent(in):: self CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - class(boundaryGeneric), optional, intent(in):: self INTEGER:: i - select type(self) - type is(boundaryWallTemperature) - !Modifies particle velocity according to wall temperature - DO i = 1, 3 - part%v(i) = part%v(i) + self%vTh*randomMaxwellian() + !Modifies particle velocity according to wall temperature + DO i = 1, 3 + part%v(i) = part%v(i) + self%vTh*randomMaxwellian() - END DO - - CALL reflection(edge, part) - - end select + END DO + + CALL genericReflection(edge, part) END SUBROUTINE wallTemperature !Ionization surface: an electron will pass through the surface ! and create an ion-electron pair based on a neutral background - SUBROUTINE ionization(edge, part, self) + SUBROUTINE ionization(self, edge, part) USE moduleList USE moduleSpecies USE moduleMesh @@ -313,9 +294,9 @@ submodule(moduleMesh) boundary USE moduleMath IMPLICIT NONE + class(boundaryIonization), intent(in):: self CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - class(boundaryGeneric), optional, intent(in):: self REAL(8):: vRel, eRel, mRel !relative velocity, energy and mass INTEGER:: nIonizations !Number of ionizations based on eRel REAL(8):: pIonization !Probability of ionization of each event @@ -324,94 +305,90 @@ submodule(moduleMesh) boundary TYPE(particle), POINTER:: newElectron TYPE(particle), POINTER:: newIon - select type(self) - type is(boundaryIonization) - mRel = reducedMass(self%m0, part%species%m) - vRel = SUM(DABS(part%v-self%v0)) - eRel = mRel*vRel**2*5.D-1 + mRel = reducedMass(self%m0, part%species%m) + vRel = SUM(DABS(part%v-self%v0)) + eRel = mRel*vRel**2*5.D-1 - !Maximum number of possible ionizations based on relative energy - nIonizations = FLOOR(eRel/self%eThreshold) + !Maximum number of possible ionizations based on relative energy + nIonizations = FLOOR(eRel/self%eThreshold) - DO p = 1, nIonizations - !Get probability of ionization - pIonization = 1.D0 - DEXP(-self%n0*self%crossSection%get(eRel)*vRel*self%effectiveTime/REAL(nIonizations)) + DO p = 1, nIonizations + !Get probability of ionization + pIonization = 1.D0 - DEXP(-self%n0*self%crossSection%get(eRel)*vRel*self%effectiveTime/REAL(nIonizations)) - !If a random number is below the probability of ionization, create new pair of ion-electron - IF (random() < pIonization) THEN - !Assign random velocity to the neutral - v0(1) = self%v0(1) + self%vTh*randomMaxwellian() - v0(2) = self%v0(2) + self%vTh*randomMaxwellian() - v0(3) = self%v0(3) + self%vTh*randomMaxwellian() + !If a random number is below the probability of ionization, create new pair of ion-electron + IF (random() < pIonization) THEN + !Assign random velocity to the neutral + v0(1) = self%v0(1) + self%vTh*randomMaxwellian() + v0(2) = self%v0(2) + self%vTh*randomMaxwellian() + v0(3) = self%v0(3) + self%vTh*randomMaxwellian() - !Allocates the new particles - ALLOCATE(newElectron) - ALLOCATE(newIon) + !Allocates the new particles + ALLOCATE(newElectron) + ALLOCATE(newIon) - IF (ASSOCIATED(self%electronSecondary)) THEN - newElectron%species => self%electronSecondary + IF (ASSOCIATED(self%electronSecondary)) THEN + newElectron%species => self%electronSecondary - ELSE - newElectron%species => part%species - - END IF - newIon%species => self%species - - newElectron%v = v0 + (1.D0 + self%deltaV*v0/NORM2(v0)) - newIon%v = v0 - - newElectron%r = edge%randPos() - newIon%r = newElectron%r - - IF (ASSOCIATED(edge%e1)) THEN - newElectron%cell = edge%e1%n - - ELSEIF (ASSOCIATED(edge%e2)) THEN - newElectron%cell = edge%e2%n - - END IF - newIon%cell = newElectron%cell - - newElectron%Xi = mesh%cells(part%cell)%obj%phy2log(newElectron%r) - newIon%Xi = newElectron%Xi - - newElectron%weight = part%weight - newIon%weight = newElectron%weight - - newElectron%n_in = .TRUE. - newIon%n_in = .TRUE. - - !Add particles to list - CALL partSurfaces%setLock() - CALL partSurfaces%add(newElectron) - CALL partSurfaces%add(newIon) - CALL partSurfaces%unsetLock() - - !Electron loses energy due to ionization - eRel = eRel - self%eThreshold - vRel = 2.D0*DSQRT(eRel)/mRel - - !Reduce number of possible ionizations - nIonizations = nIonizations - 1 + ELSE + newElectron%species => part%species END IF + newIon%species => self%species - END DO + newElectron%v = v0 + (1.D0 + self%deltaV*v0/NORM2(v0)) + newIon%v = v0 - !Removes ionizing electron regardless the number of pair created - part%n_in = .FALSE. + newElectron%r = edge%randPos() + newIon%r = newElectron%r - end select + IF (ASSOCIATED(edge%e1)) THEN + newElectron%cell = edge%e1%n + + ELSEIF (ASSOCIATED(edge%e2)) THEN + newElectron%cell = edge%e2%n + + END IF + newIon%cell = newElectron%cell + + newElectron%Xi = mesh%cells(part%cell)%obj%phy2log(newElectron%r) + newIon%Xi = newElectron%Xi + + newElectron%weight = part%weight + newIon%weight = newElectron%weight + + newElectron%n_in = .TRUE. + newIon%n_in = .TRUE. + + !Add particles to list + CALL partSurfaces%setLock() + CALL partSurfaces%add(newElectron) + CALL partSurfaces%add(newIon) + CALL partSurfaces%unsetLock() + + !Electron loses energy due to ionization + eRel = eRel - self%eThreshold + vRel = 2.D0*DSQRT(eRel)/mRel + + !Reduce number of possible ionizations + nIonizations = nIonizations - 1 + + END IF + + END DO + + !Removes ionizing electron regardless the number of pair created + part%n_in = .FALSE. END SUBROUTINE ionization - subroutine quasiNeutrality(edge, part, self) + subroutine quasiNeutrality(self, edge, part) use moduleRandom implicit none + class(boundaryQuasiNeutrality), intent(in):: self class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part - class(boundaryGeneric), optional, intent(in):: self real(8), allocatable:: density(:) class(meshCell), pointer:: cell real(8):: EF_dir @@ -427,50 +404,51 @@ submodule(moduleMesh) boundary end if if (random() <= alpha) then - call reflection(edge, part) + call genericReflection(edge, part) else - call transparent(edge, part) + call genericTransparent(edge, part) end if end subroutine quasiNeutrality -! !Points the boundary function to specific type -! module SUBROUTINE pointBoundaryFunction(edge, s) -! USE moduleErrors -! IMPLICIT NONE -! -! CLASS(meshEdge), INTENT(inout):: edge -! INTEGER, INTENT(in):: s !Species index -! -! SELECT TYPE(obj => edge%boundary%bTypes(s)%obj) -! TYPE IS(boundaryAbsorption) -! edge%fBoundary(s)%apply => absorption -! -! TYPE IS(boundaryReflection) -! edge%fBoundary(s)%apply => reflection -! -! TYPE IS(boundaryTransparent) -! edge%fBoundary(s)%apply => transparent -! -! TYPE IS(boundaryAxis) -! edge%fBoundary(s)%apply => symmetryAxis -! -! TYPE IS(boundaryWallTemperature) -! edge%fBoundary(s)%apply => wallTemperature -! -! TYPE IS(boundaryIonization) -! edge%fBoundary(s)%apply => ionization -! -! type is(boundaryQuasiNeutrality) -! edge%fBoundary(s)%apply => quasiNeutrality -! -! CLASS DEFAULT -! CALL criticalError("Boundary type not defined", 'pointBoundaryFunction') -! -! END SELECT -! -! END SUBROUTINE pointBoundaryFunction + ! 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 + !rp = intersection between particle and edge + !rpp = final position of particle + !vpp = final velocity of particle + real(8), dimension(1:3):: rp, vpp + + !Reflect particle velocity + vpp = part%v - 2.D0*dot_product(part%v, edge%normal)*edge%normal + part%v = vpp + + rp = edge%intersection(part%r) + + part%r = 2.D0*(rp - part%r) + part%r + + !particle is assumed to be inside + part%n_in = .TRUE. + + end subroutine genericReflection + + subroutine genericTransparent(edge, part) + use moduleSpecies + implicit none + + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part + + !Removes particle from domain + part%n_in = .FALSE. + + end subroutine genericTransparent end submodule boundary