From 4c72e68246771b2893d4b9f85d0036789b018692 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 18 Feb 2026 13:31:38 +0100 Subject: [PATCH] Boundaries for EM almost done --- src/makefile | 2 +- src/modules/init/moduleInput.f90 | 96 ++------ .../mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 | 2 +- .../mesh/inout/text/moduleMeshInputText.f90 | 2 +- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 2 +- src/modules/mesh/moduleMesh.f90 | 75 +++++-- src/modules/mesh/moduleMesh@boundaryEM.f90 | 212 ++++++++++++------ .../mesh/moduleMesh@boundaryParticle.f90 | 19 +- src/modules/mesh/moduleMesh@elements.f90 | 10 + 9 files changed, 249 insertions(+), 171 deletions(-) diff --git a/src/makefile b/src/makefile index de27c82..49c12cb 100644 --- a/src/makefile +++ b/src/makefile @@ -1,4 +1,4 @@ -OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/moduleCompTime.o \ +OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshCommon.o $(OBJDIR)/moduleCompTime.o \ $(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \ $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index f41ac07..ed34f3f 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -38,10 +38,15 @@ MODULE moduleInput CALL readSpecies(config) CALL checkStatus(config, "readSpecies") - !Read boundaries - CALL verboseError('Reading boundary conditions...') - CALL readBoundary(config) - CALL checkStatus(config, "readBoundary") + !Read particle boundaries + CALL verboseError('Reading particle boundary conditions...') + CALL readBoundaryParticle(config) + CALL checkStatus(config, "readBoundaryParticle") + + ! read EM boundaries + CALL verboseError('Reading EM boundary conditions...') + CALL readBoundaryEM(config) + CALL checkStatus(config, "readBoundaryEM") !Read Geometry CALL verboseError('Reading Geometry...') @@ -254,17 +259,7 @@ MODULE moduleInput IF (found) THEN CALL solver%initEM(EMType) SELECT CASE(EMType) - CASE("Electrostatic") - !Read BC - CALL readEMBoundary(config) - - CASE("ElectrostaticBoltzmann") - !Read BC - CALL readEMBoundary(config) - CASE("ConstantB") - !Read BC - CALL readEMBoundary(config) !Read constant magnetic field DO i = 1, 3 WRITE(iString, '(i2)') i @@ -791,7 +786,7 @@ MODULE moduleInput END SUBROUTINE readInteractions !Reads boundary conditions for the mesh - SUBROUTINE readBoundary(config) + SUBROUTINE readBoundaryParticle(config) use moduleMesh USE moduleErrors use moduleSpecies, only: nSpecies @@ -837,7 +832,7 @@ MODULE moduleInput !Init the list of particles from surfaces CALL OMP_INIT_LOCK(partSurfaces%lock) - END SUBROUTINE readBoundary + END SUBROUTINE readBoundaryParticle !Read the geometry (mesh) for the case SUBROUTINE readGeometry(config) @@ -1086,7 +1081,7 @@ MODULE moduleInput END SUBROUTINE readProbes - SUBROUTINE readEMBoundary(config) + SUBROUTINE readBoundaryEM(config) USE moduleMesh USE moduleOutput USE moduleErrors @@ -1107,68 +1102,22 @@ MODULE moduleInput INTEGER:: info EXTERNAL:: dgetrf - CALL config%info('boundaryEM', found, n_children = nBoundaryEM) + CALL config%info('boundaries.EM.models', found, n_children = nBoundariesEM) IF (found) THEN - ALLOCATE(boundaryEM(1:nBoundaryEM)) + ALLOCATE(boundariesEM(1:nBoundariesEM)) END IF DO b = 1, nBoundaryEM WRITE(bString, '(I2)') b - object = 'boundaryEM(' // TRIM(bString) // ')' + object = 'boundaries.EM.models(' // TRIM(bString) // ')' - CALL config%get(object // '.type', typeEM, found) + call boundariesEM(b)%init(config, object, b) - SELECT CASE(typeEM) - CASE ("dirichlet") - CALL config%get(object // '.potential', potential, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') - - END IF - - CALL config%get(object // '.physicalSurface', physicalSurface, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', & - 'readEMBoundary') - - END IF - - CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential) - - CASE ("dirichletTime") - CALL config%get(object // '.potential', potential, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', & - 'readEMBoundary') - - END IF - - CALL config%get(object // '.temporalProfile', temporalProfile, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', & - 'readEMBoundary') - - END IF - temporalProfilePath = path // temporalProfile - - CALL config%get(object // '.physicalSurface', physicalSurface, found) - IF (.NOT. found) THEN - CALL criticalError('Required parameter "physicalSurface" for Dirichlet Time boundary condition not found', & - 'readEMBoundary') - - END IF - - CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfilePath) - - CASE DEFAULT - CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary') - - END SELECT - - END DO + end do + ! TODO: Move this to the init of species ALLOCATE(qSpecies(1:nSpecies)) DO s = 1, nSpecies SELECT TYPE(sp => species(s)%obj) @@ -1182,9 +1131,10 @@ MODULE moduleInput END DO + ! TODO: Do this after the mesh has been read ! Modify K matrix due to boundary conditions - DO b = 1, nBoundaryEM - SELECT TYPE(boundary => boundaryEM(b)%obj) + DO b = 1, nBoundariesEM + SELECT TYPE(boundary => boundariesEM(b)%obj) TYPE IS(boundaryEMDirichlet) DO n = 1, boundary%nNodes ni = boundary%nodes(n)%obj%n @@ -1208,11 +1158,11 @@ MODULE moduleInput !Compute the PLU factorization of K once boundary conditions have been read CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info) IF (info /= 0) THEN - CALL criticalError('Factorization of K matrix failed', 'readEMBoundary') + CALL criticalError('Factorization of K matrix failed', 'readBoundaryEM') END IF - END SUBROUTINE readEMBoundary + END SUBROUTINE readBoundaryEM !Reads the injection of particles from the boundaries SUBROUTINE readInject(config) diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 3752bc1..a495ab9 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -199,7 +199,7 @@ MODULE moduleMeshInputGmsh2 END SELECT !Associate boundary condition procedure. - bt = physicalSurface_to_id(boundaryType) + bt = physicalSurface_to_id(boundariesParticleLinking,boundaryType) CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex) DEALLOCATE(p) diff --git a/src/modules/mesh/inout/text/moduleMeshInputText.f90 b/src/modules/mesh/inout/text/moduleMeshInputText.f90 index c00e615..3ada5b9 100644 --- a/src/modules/mesh/inout/text/moduleMeshInputText.f90 +++ b/src/modules/mesh/inout/text/moduleMeshInputText.f90 @@ -142,7 +142,7 @@ module moduleMeshInputText allocate(p(1)) p(1) = n !Associate boundary condition procedure. - bt = physicalSurface_to_id(physicalID) + bt = physicalSurface_to_id(boundariesParticleLinking, physicalID) call self%edges(physicalID)%obj%init(physicalID, p, boundariesParticleLinking(bt)%speciesIndex) deallocate(p) diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index 03699bb..c775db8 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -365,7 +365,7 @@ MODULE moduleMeshInputVTU END SELECT !Associate boundary condition procedure. - bt = physicalSurface_to_id(entitiesID(n)) + bt = physicalSurface_to_id(boundariesParticleLinking,entitiesID(n)) !Init edge CALL self%edges(e)%obj%init(n, p, boundariesParticleLinking(bt)%speciesIndex) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 190c50c..56738f5 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -72,9 +72,20 @@ MODULE moduleMesh TYPE:: meshNodePointer CLASS(meshNode), POINTER:: obj CONTAINS + procedure, pass:: meshNodePointer_equal + generic:: operator(==) => meshNodePointer_equal END TYPE meshNodePointer + interface + module function meshNodePointer_equal(self, other) result(isEqual) + class(meshNodePointer), intent(in):: self, other + logical:: isEqual + + end function meshNodePointer_equal + + end interface + !Parent of Edge element TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge !Nomber of nodes in the edge @@ -89,7 +100,6 @@ MODULE moduleMesh REAL(8):: surface = 0.D0 !Pointer to boundary per species TYPE(boundaryParticlePointer), allocatable:: boundariesParticle(:) - class(boundaryEMGeneric), pointer:: boundaryEM => null() !Physical surface for the edge CONTAINS !DEFERED PROCEDURES @@ -118,16 +128,17 @@ MODULE moduleMesh ABSTRACT INTERFACE !Inits the edge parameters - SUBROUTINE initEdge_interface(self, n, p, bt) + subroutine initEdge_interface(self, n, p, btPart, btEM) use moduleSpecies, only:nSpecies - IMPORT:: meshEdge + import:: meshEdge - CLASS(meshEdge), INTENT(out):: self - INTEGER, INTENT(in):: n - INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt(1:nSpecies) + class(meshEdge), intent(out):: self + integer, intent(in):: n + integer, intent(in):: p(:) + integer, intent(in):: btPart(1:nSpecies) + integer, intent(in):: btEM - END SUBROUTINE initEdge_interface + end subroutine initEdge_interface !Get nodes index from node PURE FUNCTION getNodesEdge_interface(self, nNodes) RESULT(n) @@ -799,18 +810,11 @@ MODULE moduleMesh !Array for linking boundaries type(boundaryParticleLinking), allocatable, dimension(:):: boundariesParticleLinking - interface - module function physicalSurface_to_id(physicalSurface) result(b) - integer:: physicalSurface - integer:: b - - end function physicalSurface_to_id - - end interface - ! BOUNDARY ELECTROMAGNETIC DEFINITIONS ! Generic type for electromagnetic boundary conditions type, public, abstract:: boundaryEMGeneric + integer:: n + character(:), allocatable:: name integer:: nNodes type(meshNodePointer), allocatable:: nodes(:) @@ -818,6 +822,7 @@ MODULE moduleMesh contains procedure, pass:: init => initBoundaryEM procedure(applyEM_interface), deferred, pass:: apply + procedure, pass:: addNodes end type boundaryEMGeneric @@ -832,6 +837,12 @@ MODULE moduleMesh end subroutine initBoundaryEM + module subroutine addNodes(self, nodes) + class(boundaryEMGeneric), intent(inout):: self + type(meshNodePointer), intent(in):: nodes(:) + + end subroutine addNodes + end interface abstract interface @@ -850,7 +861,7 @@ MODULE moduleMesh class(boundaryEMGeneric), intent(in):: self - end subroutine updateEM_interface + end subroutine updateEM_interface end interface @@ -894,10 +905,34 @@ MODULE moduleMesh END TYPE boundaryEMCont - INTEGER:: nBoundaryEM - TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:) + type:: boundaryEMLinking + integer:: physicalSurface + class(boundaryEMGeneric), pointer:: boundary => null() + + end type boundaryEMLinking + + INTEGER:: nBoundariesEM + TYPE(boundaryEMCont), ALLOCATABLE, target:: boundariesEM(:) !Information of charge and reference parameters for rho vector REAL(8), ALLOCATABLE:: qSpecies(:) + ! Get ID from array based on physical surface + interface physicalSurface_to_id + module function physicalSurface_to_boundaryParticleId(boundaryArray, physicalSurface) result(b) + type(boundaryParticleLinking):: boundaryArray(:) + integer:: physicalSurface + integer:: b + + end function physicalSurface_to_boundaryParticleId + + module function physicalSurface_to_boundaryEMId(boundaryArray, physicalSurface) result(b) + type(boundaryEMLinking):: boundaryArray(:) + integer:: physicalSurface + integer:: b + + end function physicalSurface_to_boundaryEMId + + end interface physicalSurface_to_id + END MODULE moduleMesh diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index 3878532..c4154e9 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -1,102 +1,152 @@ submodule(moduleMesh) boundaryEM CONTAINS - SUBROUTINE findNodes(self, physicalSurface) - USE moduleMesh - IMPLICIT NONE + subroutine addNodes(self, nodes) + implicit none - CLASS(boundaryEMGeneric), INTENT(inout):: self - INTEGER, INTENT(in):: physicalSurface - CLASS(meshEdge), POINTER:: edge - INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:) - INTEGER:: nNodes, nodesNew - INTEGER:: e, n + class(boundaryEMGeneric), intent(inout):: self + type(meshNodePointer), intent(in):: nodes(:) + INTEGER:: n, nn, nNodes + logical:: inArray - !Temporal array to hold nodes - ALLOCATE(nodes(0)) - ! Loop thorugh the edges and identify those that are part of the boundary - DO e = 1, mesh%numEdges - edge => mesh%edges(e)%obj - IF (edge%physicalSurface == physicalSurface) THEN - ! Edge is of the right boundary index - ! Get nodes in the edge - nNodes = edge%nNodes - nodesEdge = edge%getNodes(nNodes) - ! Collect all nodes that are not already in the temporal array - DO n = 1, nNodes - IF (ANY(nodes == nodesEdge(n))) THEN - ! Node already in array, skip - CYCLE - - ELSE - ! If not, add element to array of nodes - nodes = [nodes, nodesEdge(n)] - - END IF - - END DO - - END IF - - END DO - - ! Point boundary to nodes - nNodes = SIZE(nodes) - ALLOCATE(self%nodes(nNodes)) - self%nNodes = nNodes + nNodes = size(nodes) + ! Collect all nodes that are not already in the boundary DO n = 1, nNodes - self%nodes(n)%obj => mesh%nodes(nodes(n))%obj + inArray = .false. + ! I cannot use the procedure 'any' for this + do nn = 1 , size(self%nodes) + IF (self%nodes(nn) == nodes(n)) THEN + ! Node already in array + inArray = .true. + exit - END DO + end if - END SUBROUTINE findNodes + end do + + if (.not. inArray) then + ! If not, add element to array of nodes + self%nodes = [self%nodes, nodes(n)] + + end if + + end do + + end subroutine addNodes + + module subroutine initBoundaryEM(self, config, object, b) + use json_module + use moduleErrors + implicit none + + class(boundaryEMGeneric), allocatable, intent(out):: self + type(json_file), intent(inout):: config + character(:), allocatable, intent(in):: object + integer, intent(in):: b + character(:), allocatable:: bType + logical:: found + + self%n = b + call config%get(object // '.name', self%name, found) + if (.not. found) then + call criticalError('Required parameter "name" for EM boundary condition not found', & + 'initBoundaryEM') + + end if + + call config%get(object // '.type', bType, found) + if (.not. found) then + call criticalError('Required parameter "type" for EM boundary condition not found', & + 'initBoundaryEM') + + end if + + select case(bType) + case ("dirichlet") + ! Allocate boundary edge + ALLOCATE(boundaryEMDirichlet:: self) + + CALL initDirichlet(self, config, object) + + case ("dirichlettime") + ! Allocate boundary edge + ALLOCATE(boundaryEMDirichletTime:: self) + + call initDirichletTime(self, config, object) + + case default + call criticalError('Boundary type ' // bType // ' not supported', 'readBoundaryEM') + + end select + + end subroutine initBoundaryEM ! Initialize Dirichlet boundary condition - SUBROUTINE initDirichlet(self, physicalSurface, potential) + SUBROUTINE initDirichlet(self, config, object) + use json_module USE moduleRefParam, ONLY: Volt_ref + use moduleErrors IMPLICIT NONE CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self - INTEGER, INTENT(in):: physicalSurface - REAL(8), INTENT(in):: potential + type(json_file), intent(inout):: config + character(:), allocatable, intent(in):: object + REAL(8):: potential + logical:: found - ! Allocate boundary edge - ALLOCATE(boundaryEMDirichlet:: self) + select type(self) + type is(boundaryEMDirichlet) + call config%get(object // '.potential', potential, found) + if (.not. found) then + call criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readBoundaryEM') + + end if - SELECT TYPE(self) - TYPE IS(boundaryEMDirichlet) self%potential = potential / Volt_ref - CALL findNodes(self, physicalSurface) + end select - END SELECT - - END SUBROUTINE initDirichlet + end subroutine initDirichlet ! Initialize Dirichlet boundary condition - SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile) - USE moduleRefParam, ONLY: Volt_ref, ti_ref - IMPLICIT NONE + subroutine initDirichletTime(self, config, object) + use json_module + use moduleRefParam, ONLY: Volt_ref, ti_ref + use moduleErrors + implicit none - CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self - INTEGER, INTENT(in):: physicalSurface - REAL(8), INTENT(in):: potential - CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile + class(boundaryEMGeneric), allocatable, intent(out):: self + type(json_file), intent(inout):: config + character(:), allocatable, intent(in):: object + real(8):: potential + character(:), allocatable:: temporalProfile + logical:: found + character(:), allocatable:: temporalProfilePath - ! Allocate boundary edge - ALLOCATE(boundaryEMDirichletTime:: self) + select type(self) + type is(boundaryEMDirichletTime) + call config%get(object // '.potential', potential, found) + if (.not. found) then + call criticalError('Required parameter "potential" for Dirichlet Time boundary condition not found', & + 'readBoundaryEM') + + end if + + call config%get(object // '.temporalProfile', temporalProfile, found) + if (.not. found) then + call criticalError('Required parameter "temporalProfile" for Dirichlet Time boundary condition not found', & + 'readBoundaryEM') + + end if + temporalProfilePath = path // temporalProfile - SELECT TYPE(self) - TYPE IS(boundaryEMDirichletTime) self%potential = potential / Volt_ref - CALL findNodes(self, physicalSurface) + call self%temporalProfile%init(temporalProfile) - CALL self%temporalProfile%init(temporalProfile) + call self%temporalProfile%convert(1.D0/ti_ref, 1.D0) - CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0) - - END SELECT + end select END SUBROUTINE initDirichletTime @@ -138,4 +188,24 @@ submodule(moduleMesh) boundaryEM END SUBROUTINE applyDirichletTime + function physicalSurface_to_boundaryEMId(boundaryArray,physicalSurface) result(b) + implicit none + + type(boundaryEMLinking):: boundaryArray(:) + integer:: physicalSurface + integer:: b + integer:: bt + + do bt=1, nPhysicalSurfaces + if (boundaryArray(bt)%physicalSurface == physicalSurface) then + b = boundaryArray(bt)%boundary%n + + exit + + end if + + end do + + end function physicalSurface_to_boundaryEMId + end submodule boundaryEM diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90 index b3225b3..8ed7ac8 100644 --- a/src/modules/mesh/moduleMesh@boundaryParticle.f90 +++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90 @@ -48,7 +48,19 @@ submodule(moduleMesh) boundaryParticle self%n = b CALL config%get(object // '.name', self%name, found) + if (.not. found) then + call criticalError('Required parameter "name" for EM boundary condition not found', & + 'initBoundaryParticle') + + end if + CALL config%get(object // '.type', bType, found) + if (.not. found) then + call criticalError('Required parameter "type" for EM boundary condition not found', & + 'initBoundaryParticle') + + end if + SELECT CASE(bType) CASE('reflection') ALLOCATE(boundaryReflection:: self) @@ -483,15 +495,16 @@ submodule(moduleMesh) boundaryParticle end subroutine genericTransparent - function physicalSurface_to_id(physicalSurface) result(b) + function physicalSurface_to_boundaryParticleId(boundaryArray, physicalSurface) result(b) implicit none + type(boundaryParticleLinking):: boundaryArray(:) integer:: physicalSurface integer:: b integer:: bt do bt=1, nPhysicalSurfaces - if (boundariesParticleLinking(bt)%physicalSurface == physicalSurface) then + if (boundaryArray(bt)%physicalSurface == physicalSurface) then b = bt exit @@ -500,6 +513,6 @@ submodule(moduleMesh) boundaryParticle end do - end function physicalSurface_to_id + end function physicalSurface_to_boundaryParticleId end submodule boundaryParticle diff --git a/src/modules/mesh/moduleMesh@elements.f90 b/src/modules/mesh/moduleMesh@elements.f90 index a83a4e2..c05e7e7 100644 --- a/src/modules/mesh/moduleMesh@elements.f90 +++ b/src/modules/mesh/moduleMesh@elements.f90 @@ -18,6 +18,16 @@ submodule(moduleMesh) elements END SUBROUTINE resetOutput + module function meshNodePointer_equal(self, other) result(isEqual) + implicit none + + class(meshNodePointer), intent(in):: self, other + logical:: isEqual + + isEqual = self%obj%n == other%obj%n + + end function meshNodePointer_equal + !Constructs the global K matrix PURE module SUBROUTINE constructGlobalK(self) IMPLICIT NONE