From eb10ae92a621b4562b7e903256d4aa85387bb3f1 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 8 Apr 2026 15:17:30 +0200 Subject: [PATCH] It seems to work, but I don't think this is quite useful for the manuscript --- src/modules/init/moduleInput.f90 | 2 + src/modules/mesh/moduleMesh.f90 | 16 ++++++ src/modules/mesh/moduleMesh@boundaryEM.f90 | 57 +++++++++++++++++----- src/modules/mesh/moduleMesh@elements.f90 | 8 +++ 4 files changed, 71 insertions(+), 12 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index be92b31..575b372 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1367,6 +1367,8 @@ MODULE moduleInput ! If the boundary for the species is linked to the one analysing, add the edges if (associated(physicalSurfaces(ps)%EM, bound)) then bound%nodes = [bound%nodes, physicalSurfaces(ps)%nodes] + + bound%edges = [bound%edges, physicalSurfaces(ps)%edges] end if end do diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index d27c1db..aa20a2b 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -964,6 +964,7 @@ MODULE moduleMesh type(meshNodePointer), allocatable:: nodes(:) procedure(updateEM_interface), pointer, pass:: update => null() + procedure(printEM_interface), pointer, pass:: print => null() contains procedure(applyEM_interface), deferred, pass:: apply @@ -987,6 +988,15 @@ MODULE moduleMesh end subroutine updateEM_interface + ! Write the values of the particle boundary model + subroutine printEM_interface(self, fileID) + import boundaryEMGeneric + + class(boundaryEMGeneric), intent(inout):: self + integer, intent(in):: fileID + + end subroutine printEM_interface + end interface ! Extended types @@ -1049,6 +1059,8 @@ MODULE moduleMesh ! Floating: Floating surface, ie, zero net current type, extends(boundaryEMGeneric):: boundaryEMFloating real(8):: potential + real(8):: invC ! Inverse of the capacitance of the surface + type(meshEdgePointer), allocatable:: edges(:) !Array with edges contains ! boundaryEMGeneric DEFERRED PROCEDURES procedure, pass:: apply => applyFloating @@ -1097,6 +1109,10 @@ MODULE moduleMesh end subroutine boundariesEM_update + module subroutine boundariesEM_write() + + end subroutine boundariesEM_write + end interface ! PHYSICAL SURFACES LINKING TO MESH ELEMENTS diff --git a/src/modules/mesh/moduleMesh@boundaryEM.f90 b/src/modules/mesh/moduleMesh@boundaryEM.f90 index 3b52c78..e45f996 100644 --- a/src/modules/mesh/moduleMesh@boundaryEM.f90 +++ b/src/modules/mesh/moduleMesh@boundaryEM.f90 @@ -127,14 +127,15 @@ submodule(moduleMesh) boundaryEM ! Init module subroutine initFloating(self, config, object) use json_module - use moduleRefParam, only: Volt_ref + use moduleRefParam, only: Volt_ref, L_ref + use moduleConstParam, only: eps_0 use moduleErrors implicit none class(boundaryEMGeneric), allocatable, intent(inout):: self type(json_file), intent(inout):: config character(:), allocatable, intent(in):: object - real(8):: potential + real(8):: potential, capacitance logical:: found select type(self) @@ -148,8 +149,20 @@ submodule(moduleMesh) boundaryEM self%potential = potential / Volt_ref + call config%get(object // '.capacitance', capacitance, found) + if (.not. found) then + call warningError('Parameter "capacitance" for floating boundary condition not found. Using the default value of 1 nF.') + capacitance = 1.0d-9 + + end if + + ! Inverse of non-dimensional capacitance + self%invC = 1.0d0 / (capacitance / (eps_0 * L_ref)) + self%update => updateFloating + allocate(self%edges(0)) + end select end subroutine initFloating @@ -174,40 +187,60 @@ submodule(moduleMesh) boundaryEM ! Update subroutine updateFloating(self) use moduleMesh, only: qSpecies, meshNode + use moduleCaseParam, only: tauMin implicit none class(boundaryEMGeneric), intent(inout):: self - integer:: n, s + integer:: e, n, s integer, allocatable:: nodes(:) + class(meshEdge), pointer:: edge + real(8), allocatable:: mom_nodes(:) class(meshNode), pointer:: node - real(8):: momNode, momNodeTotal, totalCurrent + real(8):: mom_center, edgeDensityCurrent, totalCurrent totalCurrent = 0.0d0 select type(self) type is(boundaryEMFloating) - do n = 1, self%nNodes - momNodeTotal = 0.0d0 + do e = 1, size(self%edges) + edge => self%edges(e)%obj - node => self%nodes(n)%obj + edgeDensityCurrent = 0.0d0 + nodes = edge%getNodes(edge%nNodes) + allocate(mom_nodes(1:edge%nNodes)) do s = 1, nSpecies - print *, s, node%output(s)%mom - momNode = norm2(node%output(s)%mom) + mom_center = 0.0d0 + do n = 1, self%nNodes + node => mesh%nodes(nodes(n))%obj - momNodeTotal = momNodeTotal + qSpecies(s) * momNode + ! Minus sign to get the flux exiting the edge + mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal) + + end do + + mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) + + edgeDensityCurrent = edgeDensityCurrent + qSpecies(s) * mom_center end do - totalCurrent = totalCurrent + momNodeTotal + totalCurrent = totalCurrent + edgeDensityCurrent * edge%surface end do - print *, totalCurrent + self%potential = self%potential + self%invC * totalCurrent * tauMin + + ! print *, totalCurrent, self%potential end select end subroutine updateFloating + ! Write + subroutine writeFloating() + + end subroutine writeFloating + ! Get the index of the boundary based on the name module function boundaryEMName_to_Index(boundaryName) result(bp) use moduleErrors diff --git a/src/modules/mesh/moduleMesh@elements.f90 b/src/modules/mesh/moduleMesh@elements.f90 index 5eadf28..b947278 100644 --- a/src/modules/mesh/moduleMesh@elements.f90 +++ b/src/modules/mesh/moduleMesh@elements.f90 @@ -176,6 +176,14 @@ submodule(moduleMesh) elements END DO + TYPE IS(boundaryEMFloating) + DO n = 1, boundary%nNodes + ni = boundary%nodes(n)%obj%n + self%K(ni, :) = 0.D0 + self%K(ni, ni) = 1.D0 + + END DO + END SELECT END DO