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.
This commit is contained in:
Jorge Gonzalez 2024-07-12 23:30:35 +02:00
commit 10dee05922
2 changed files with 92 additions and 46 deletions

View file

@ -1148,8 +1148,8 @@ MODULE moduleInput
CHARACTER(:), ALLOCATABLE:: typeEM CHARACTER(:), ALLOCATABLE:: typeEM
REAL(8):: potential REAL(8):: potential
INTEGER:: physicalSurface INTEGER:: physicalSurface
CHARACTER(:), ALLOCATABLE:: timeProfile CHARACTER(:), ALLOCATABLE:: temporalProfile
INTEGER:: b, e, s, n, ni INTEGER:: b, s, n, ni
CHARACTER(2):: bString CHARACTER(2):: bString
INTEGER:: info INTEGER:: info
EXTERNAL:: dgetrf EXTERNAL:: dgetrf
@ -1187,12 +1187,20 @@ MODULE moduleInput
END IF 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) CALL config%get(object // '.physicalSurface', physicalSurface, found)
IF (.NOT. found) THEN IF (.NOT. found) THEN
CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
END IF END IF
CALL initDirichletTime(boundaryEM(b)%obj, physicalSurface, potential, temporalProfile)
CASE DEFAULT CASE DEFAULT
CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary') CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary')

View file

@ -66,20 +66,65 @@ MODULE moduleEM
REAL(8), ALLOCATABLE:: qSpecies(:) REAL(8), ALLOCATABLE:: qSpecies(:)
CONTAINS CONTAINS
! Initialize Dirichlet boundary condition SUBROUTINE findNodes(self, physicalSurface)
SUBROUTINE initDirichlet(self, physicalSurface, potential)
USE moduleMesh USE moduleMesh
USE moduleRefParam, ONLY: Volt_ref
IMPLICIT NONE IMPLICIT NONE
CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self CLASS(boundaryEMGeneric), INTENT(inout):: self
INTEGER:: physicalSurface INTEGER, INTENT(in):: physicalSurface
REAL(8), INTENT(in):: potential
CLASS(meshEdge), POINTER:: edge CLASS(meshEdge), POINTER:: edge
INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:) INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:)
INTEGER:: nNodes, nodesNew INTEGER:: nNodes, nodesNew
INTEGER:: e, n 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 boundary edge
ALLOCATE(boundaryEMDirichlet:: self) ALLOCATE(boundaryEMDirichlet:: self)
@ -87,48 +132,39 @@ MODULE moduleEM
TYPE IS(boundaryEMDirichlet) TYPE IS(boundaryEMDirichlet)
self%potential = potential / Volt_ref self%potential = potential / Volt_ref
!TODO: This is going into a function CALL findNodes(self, physicalSurface)
!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 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
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 !Apply Dirichlet boundary condition to the poisson equation
SUBROUTINE applyDirichlet(self, vectorF) SUBROUTINE applyDirichlet(self, vectorF)
USE moduleMesh USE moduleMesh
@ -140,6 +176,7 @@ MODULE moduleEM
DO n = 1, self%nNodes DO n = 1, self%nNodes
self%nodes(n)%obj%emData%phi = self%potential self%nodes(n)%obj%emData%phi = self%potential
vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi
END DO END DO
@ -156,7 +193,8 @@ MODULE moduleEM
INTEGER:: n, ni INTEGER:: n, ni
DO n = 1, self%nNodes 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 END DO