From 10dee05922e17d72258fbd432ed1c689f08571ba Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 12 Jul 2024 23:30:35 +0200 Subject: [PATCH] NOT WORKING: Compilation okay, but not Dirichlet BC The code compiles but the right BC is not being applied to the vectorF. I'll check this tomorrow. --- src/modules/init/moduleInput.f90 | 12 +- .../solver/electromagnetic/moduleEM.f90 | 126 ++++++++++++------ 2 files changed, 92 insertions(+), 46 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 6800ac1..eea1937 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1148,8 +1148,8 @@ MODULE moduleInput CHARACTER(:), ALLOCATABLE:: typeEM REAL(8):: potential INTEGER:: physicalSurface - CHARACTER(:), ALLOCATABLE:: timeProfile - INTEGER:: b, e, s, n, ni + CHARACTER(:), ALLOCATABLE:: temporalProfile + INTEGER:: b, s, n, ni CHARACTER(2):: bString INTEGER:: info EXTERNAL:: dgetrf @@ -1187,12 +1187,20 @@ MODULE moduleInput END IF + CALL config%get(object // '.temporalProfile', temporalProfile, 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 initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfile) + CASE DEFAULT CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary') diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index 30a23ad..b480330 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -66,20 +66,65 @@ MODULE moduleEM REAL(8), ALLOCATABLE:: qSpecies(:) CONTAINS - ! Initialize Dirichlet boundary condition - SUBROUTINE initDirichlet(self, physicalSurface, potential) + SUBROUTINE findNodes(self, physicalSurface) USE moduleMesh - USE moduleRefParam, ONLY: Volt_ref IMPLICIT NONE - CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self - INTEGER:: physicalSurface - REAL(8), INTENT(in):: potential + CLASS(boundaryEMGeneric), INTENT(inout):: self + INTEGER, INTENT(in):: physicalSurface CLASS(meshEdge), POINTER:: edge INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:) INTEGER:: nNodes, nodesNew INTEGER:: e, n + !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)) + DO n = 1, nNodes + self%nodes(n)%obj => mesh%nodes(nodes(n))%obj + + END DO + + END SUBROUTINE findNodes + + ! Initialize Dirichlet boundary condition + SUBROUTINE initDirichlet(self, physicalSurface, potential) + USE moduleRefParam, ONLY: Volt_ref + IMPLICIT NONE + + CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self + INTEGER, INTENT(in):: physicalSurface + REAL(8), INTENT(in):: potential + ! Allocate boundary edge ALLOCATE(boundaryEMDirichlet:: self) @@ -87,48 +132,39 @@ MODULE moduleEM TYPE IS(boundaryEMDirichlet) self%potential = potential / Volt_ref - !TODO: This is going into a function - !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)) - DO n = 1, nNodes - self%nodes(n)%obj => mesh%nodes(nodes(n))%obj - - END DO + CALL findNodes(self, physicalSurface) END SELECT END SUBROUTINE initDirichlet + ! Initialize Dirichlet boundary condition + SUBROUTINE initDirichletTime(self, physicalSurface, potential, temporalProfile) + USE moduleRefParam, ONLY: Volt_ref, ti_ref + IMPLICIT NONE + + CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self + INTEGER, INTENT(in):: physicalSurface + REAL(8), INTENT(in):: potential + CHARACTER(:), ALLOCATABLE, INTENT(in):: temporalProfile + + ! Allocate boundary edge + ALLOCATE(boundaryEMDirichletTime:: self) + + SELECT TYPE(self) + TYPE IS(boundaryEMDirichletTime) + self%potential = potential / Volt_ref + + CALL findNodes(self, physicalSurface) + + CALL self%temporalProfile%init(temporalProfile) + + CALL self%temporalProfile%convert(1.D0/ti_ref, 1.D0) + + END SELECT + + END SUBROUTINE initDirichletTime + !Apply Dirichlet boundary condition to the poisson equation SUBROUTINE applyDirichlet(self, vectorF) USE moduleMesh @@ -140,6 +176,7 @@ MODULE moduleEM DO n = 1, self%nNodes self%nodes(n)%obj%emData%phi = self%potential + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi END DO @@ -156,7 +193,8 @@ MODULE moduleEM INTEGER:: n, ni DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential + self%nodes(n)%obj%emData%phi = self%potential !TODO: Correct for time + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi END DO