From fa1de0accf08ee312e839ea5c771c6458d10ae7b Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Fri, 14 Apr 2023 12:25:32 +0200 Subject: [PATCH 01/11] bias voltage WIP --- .../features-and-models/BC-field-solver.md | 136 +++---- src/equations/poisson/equation.f90 | 169 +++++---- src/hdg/hdg.f90 | 344 ++++++++++++++++-- src/mesh/prepare_mesh.f90 | 2 +- 4 files changed, 492 insertions(+), 159 deletions(-) diff --git a/docs/documentation/userguide/features-and-models/BC-field-solver.md b/docs/documentation/userguide/features-and-models/BC-field-solver.md index f6b27093b..6385fcd3b 100644 --- a/docs/documentation/userguide/features-and-models/BC-field-solver.md +++ b/docs/documentation/userguide/features-and-models/BC-field-solver.md @@ -19,39 +19,39 @@ In this case the boundary type is changed from 4 (in the mesh file) to 5 in the The boundary conditions used for Maxwell's equations are defined by the first integer value in the *BoundaryType* vector (consisting of the *Type* and *State*) and include, periodic, Dirichlet, Silver-Mueller, perfectly conducting, symmetry and reference state boundaries as detailed in the following table. -| BoundaryType | Type | State | -| :----------: | :-------: | :------------------------------------------------------------------------------------------------------------------------ | -| (/1,1/) | periodic | 1: positive direction of the 1st periodicity vector | -| (/1,-1/) | periodic | -1: negative (opposite) direction of the 1st periodicity vector | -| | | | -| (/2,2/) | Dirichlet | 2: Coaxial waveguide | -| (/2,22/) | Dirichlet | 22: Coaxial waveguide BC (boundary condition or exact flux) | -| (/2,3/) | Dirichlet | 3: Resonator | -| (/2,4/) | Dirichlet | 4: Electromagnetic dipole (implemented via RHS source terms and shape function deposition) | -| (/2,40/) | Dirichlet | 40: Electromagnetic dipole without initial condition (implemented via RHS source terms and shape function deposition) | -| (/2,41/) | Dirichlet | 41: Pulsed Electromagnetic dipole (implemented via RHS source terms and shape function deposition) | -| (/2,5/) | Dirichlet | 5: Transversal Electric (TE) plane wave in a circular waveguide | -| (/2,7/) | Dirichlet | 7: Special manufactured Solution | -| (/2,10/) | Dirichlet | 10: Issautier 3D test case with source (Stock et al., div. correction paper), domain [0;1]^3 | -| (/2,12/) | Dirichlet | 12: Plane wave | -| (/2,121/) | Dirichlet | 121: Pulsed plane wave (infinite spot size) and temporal Gaussian | -| (/2,14/) | Dirichlet | 14: Gaussian pulse is initialized inside the domain (usually used as initial condition and not BC) | -| (/2,15/) | Dirichlet | 15: Gaussian pulse with optional delay time *tDelayTime* | -| (/2,16/) | Dirichlet | 16: Gaussian pulse which is initialized in the domain and used as a boundary condition for t>0 | -| (/2,50/) | Dirichlet | 50: Initialization and BC Gyrotron - including derivatives | -| (/2,51/) | Dirichlet | 51: Initialization and BC Gyrotron - including derivatives (nothing is set for z>eps) | -| | | | -| (/3,0/) | SM | 1st order absorbing BC (Silver-Mueller) - Munz et al. 2000 / Computer Physics Communication 130, 83-117 with fix | -| | | of div. correction field for low B-fields that only set the correction fields when ABS(B)>1e-10 | -| (/5,0/) | SM | 1st order absorbing BC (Silver-Mueller) - Munz et al. 2000 / Computer Physics Communication 130, 83-117 | -| (/6,0/) | SM | 1st order absorbing BC (Silver-Mueller) - Munz et al. 2000 / Computer Physics Communication 130, 83-117 with fix | -| | | of div. correction field for low B-fields that only set the correction fields when B is significantly large compared to E | -| | | | -| (/4,0/) | PEC | Perfectly electric conducting surface (Munz, Omnes, Schneider 2000, pp. 97-98) | -| | | | -| (/10,0/) | Symmetry | Symmetry BC (perfect MAGNETIC conductor, PMC) | -| | | | -| (/20,0/) | Ref | Use state that is read from .h5 file and interpolated to the BC | +| (/BCType,BCState/) | Type | State | +| :----------------: | :-------: | :------------------------------------------------------------------------------------------------------------------------ | +| (/1,1/) | periodic | 1: positive direction of the 1st periodicity vector | +| (/1,-1/) | periodic | -1: negative (opposite) direction of the 1st periodicity vector | +| | | | +| (/2,2/) | Dirichlet | 2: Coaxial waveguide | +| (/2,22/) | Dirichlet | 22: Coaxial waveguide BC (boundary condition or exact flux) | +| (/2,3/) | Dirichlet | 3: Resonator | +| (/2,4/) | Dirichlet | 4: Electromagnetic dipole (implemented via RHS source terms and shape function deposition) | +| (/2,40/) | Dirichlet | 40: Electromagnetic dipole without initial condition (implemented via RHS source terms and shape function deposition) | +| (/2,41/) | Dirichlet | 41: Pulsed Electromagnetic dipole (implemented via RHS source terms and shape function deposition) | +| (/2,5/) | Dirichlet | 5: Transversal Electric (TE) plane wave in a circular waveguide | +| (/2,7/) | Dirichlet | 7: Special manufactured Solution | +| (/2,10/) | Dirichlet | 10: Issautier 3D test case with source (Stock et al., div. correction paper), domain [0;1]^3 | +| (/2,12/) | Dirichlet | 12: Plane wave | +| (/2,121/) | Dirichlet | 121: Pulsed plane wave (infinite spot size) and temporal Gaussian | +| (/2,14/) | Dirichlet | 14: Gaussian pulse is initialized inside the domain (usually used as initial condition and not BC) | +| (/2,15/) | Dirichlet | 15: Gaussian pulse with optional delay time *tDelayTime* | +| (/2,16/) | Dirichlet | 16: Gaussian pulse which is initialized in the domain and used as a boundary condition for t>0 | +| (/2,50/) | Dirichlet | 50: Initialization and BC Gyrotron - including derivatives | +| (/2,51/) | Dirichlet | 51: Initialization and BC Gyrotron - including derivatives (nothing is set for z>eps) | +| | | | +| (/3,0/) | SM | 1st order absorbing BC (Silver-Mueller) - Munz et al. 2000 / Computer Physics Communication 130, 83-117 with fix | +| | | of div. correction field for low B-fields that only set the correction fields when ABS(B)>1e-10 | +| (/5,0/) | SM | 1st order absorbing BC (Silver-Mueller) - Munz et al. 2000 / Computer Physics Communication 130, 83-117 | +| (/6,0/) | SM | 1st order absorbing BC (Silver-Mueller) - Munz et al. 2000 / Computer Physics Communication 130, 83-117 with fix | +| | | of div. correction field for low B-fields that only set the correction fields when B is significantly large compared to E | +| | | | +| (/4,0/) | PEC | Perfectly electric conducting surface (Munz, Omnes, Schneider 2000, pp. 97-98) | +| | | | +| (/10,0/) | Symmetry | Symmetry BC (perfect MAGNETIC conductor, PMC) | +| | | | +| (/20,0/) | Ref | Use state that is read from .h5 file and interpolated to the BC | Dielectric -> type 100? @@ -61,41 +61,41 @@ The boundary conditions used for Maxwell's equations are defined by the first in include, periodic, Dirichlet (via pre-defined function, zero-potential or *RefState*), Neumann and reference state boundaries as detailed in the following table. -| BoundaryType | Type | State | -| :----------: | :-------: | :------------------------------------------------------------------------------------------------------------------------- | -| (/1,1/) | periodic | 1: positive direction of the 1st periodicity vector | -| (/1,-1/) | periodic | -1: negative (opposite) direction of the 1st periodicity vector | -| | | | -| (/2,0/) | Dirichlet | 0: Phi=0 | -| (/2,2/) | Dirichlet | 2: Automatic adjustment for Phi to meet const. input power, see {ref}`sec:fixed-coupled-power` | -| (/2,1001/) | Dirichlet | 1001: linear potential y-z via Phi = 2340y + 2340z | -| (/2,101/) | Dirichlet | 101: linear in z-direction: z=-1: 0, z=1, 1000 | -| (/2,103/) | Dirichlet | 103: dipole | -| (/2,104/) | Dirichlet | 104: solution to Laplace's equation: Phi_xx + Phi_yy + Phi_zz = 0 | -| | | $\Phi=(COS(x)+SIN(x))(COS(y)+SIN(y))(COSH(SQRT(2.0)z)+SINH(SQRT(2.0)z))$ | -| (/2,200/) | Dirichlet | 200: Dielectric Sphere of Radius R in constant electric field E_0 from book: John David Jackson, Classical Electrodynamics | -| (/2,300/) | Dirichlet | 300: Dielectric Slab in z-direction of half width R in constant electric field E_0: | -| | | adjusted from CASE(200) | -| (/2,301/) | Dirichlet | 301: like CASE=300, but only in positive z-direction the dielectric region is assumed | -| (/2,400/) | Dirichlet | 400: Point Source in Dielectric Region with | -| | | epsR_1 = 1 for x $<$ 0 (vacuum) | -| | | epsR_2 != 1 for x $>$ 0 (dielectric region) | -| | | | -| (/4,0/) | Dirichlet | zero-potential (Phi=0) | -| | | | -| (/5,1/) | Dirichlet | 1: use RefState Nbr 1, see details below | -| | | | -| (/6,1/) | Dirichlet | 1: use RefState Nbr 1, see details below | -| | | | -| (/7,1/) | Dirichlet | 1: use LinState Nbr 1, linear function for Phi, see {ref}`sec:linear-potential` | -| | | | -| (/8,1/) | Dirichlet | 8: Assign BC to EPC group nbr. 1 (different BCs can be assigned the same EPC), see {ref}`sec:electric-potential-condition` | -| | | | -| (/10,0/) | Neumann | zero-gradient (dPhi/dn=0) | -| | | | -| (/11,0/) | Neumann | q*n=1 | -| | | | -| (/20,1/) | FPC | 1: Assign BC to FPC group nbr. 1 (different BCs can be assigned the same FPC), see {ref}`sec:floating-boundary-condition` | +| (/BCType,BCState/) | Type | State | +| :----------------: | :-------: | :------------------------------------------------------------------------------------------------------------------------- | +| (/1,1/) | periodic | 1: positive direction of the 1st periodicity vector | +| (/1,-1/) | periodic | -1: negative (opposite) direction of the 1st periodicity vector | +| | | | +| (/2,0/) | Dirichlet | 0: Phi=0 | +| (/2,2/) | Dirichlet | 2: Automatic adjustment for Phi to meet const. input power, see {ref}`sec:fixed-coupled-power` | +| (/2,1001/) | Dirichlet | 1001: linear potential y-z via Phi = 2340y + 2340z | +| (/2,101/) | Dirichlet | 101: linear in z-direction: z=-1: 0, z=1, 1000 | +| (/2,103/) | Dirichlet | 103: dipole | +| (/2,104/) | Dirichlet | 104: solution to Laplace's equation: Phi_xx + Phi_yy + Phi_zz = 0 | +| | | $\Phi=(COS(x)+SIN(x))(COS(y)+SIN(y))(COSH(SQRT(2.0)z)+SINH(SQRT(2.0)z))$ | +| (/2,200/) | Dirichlet | 200: Dielectric Sphere of Radius R in constant electric field E_0 from book: John David Jackson, Classical Electrodynamics | +| (/2,300/) | Dirichlet | 300: Dielectric Slab in z-direction of half width R in constant electric field E_0: | +| | | adjusted from CASE(200) | +| (/2,301/) | Dirichlet | 301: like CASE=300, but only in positive z-direction the dielectric region is assumed | +| (/2,400/) | Dirichlet | 400: Point Source in Dielectric Region with | +| | | epsR_1 = 1 for x $<$ 0 (vacuum) | +| | | epsR_2 != 1 for x $>$ 0 (dielectric region) | +| | | | +| (/4,0/) | Dirichlet | zero-potential (Phi=0) | +| | | | +| (/5,1/) | Dirichlet | 1: use RefState Nbr 1 and $\cos(\omega t)$ function (see details below) | +| | | | +| (/6,1/) | Dirichlet | 1: use RefState Nbr 1 and $\cos(\omega t)+1$ function that does not switch sign (see details below) | +| | | | +| (/7,1/) | Dirichlet | 1: use LinState Nbr 1, linear function for Phi, see {ref}`sec:linear-potential` | +| | | | +| (/8,1/) | Dirichlet | 8: Assign BC to EPC group nbr. 1 (different BCs can be assigned the same EPC), see {ref}`sec:electric-potential-condition` | +| | | | +| (/10,0/) | Neumann | zero-gradient (dPhi/dn=0) | +| | | | +| (/11,0/) | Neumann | q*n=1 | +| | | | +| (/20,1/) | FPC | 1: Assign BC to FPC group nbr. 1 (different BCs can be assigned the same FPC), see {ref}`sec:floating-boundary-condition` | ### RefState boundaries {-} diff --git a/src/equations/poisson/equation.f90 b/src/equations/poisson/equation.f90 index 53ed9df98..c8bbd8c87 100644 --- a/src/equations/poisson/equation.f90 +++ b/src/equations/poisson/equation.f90 @@ -84,7 +84,7 @@ SUBROUTINE DefineParametersEquation() ! Special BC with floating potential that is defined by the absorbed power of the charged particles CALL prms%CreateRealArrayOption('CoupledPowerPotential' , 'Controlled power input: Supply vector of form (/min, start, max/) for the minimum, start (t=0) and maximum electric potential that is applied at BoundaryType = (/2,2/).', no=3 ) CALL prms%CreateRealOption( 'CoupledPowerTarget' , 'Controlled power input: Target input power to which the electric potential is adjusted for BoundaryType = (/2,2/)' ) -CALL prms%CreateRealOption( 'CoupledPowerRelaxFac' , 'Relaxation factor for calculation of new electric potential due to defined Target input power. Default = 0.05 (5%)', '0.05' ) +CALL prms%CreateRealOption( 'CoupledPowerRelaxFac' , 'Relaxation factor for calculation of new electric potential due to defined Target input power. Default = 0.05 (which is 5%)', '0.05' ) #endif /*defined(PARTICLES)*/ END SUBROUTINE DefineParametersEquation @@ -127,12 +127,6 @@ SUBROUTINE InitEquation() CHARACTER(LEN=255) :: BCName INTEGER :: nRefStateMax INTEGER :: nLinState,nLinStateMax -LOGICAL :: PCouplAlreadySet -#if defined(PARTICLES) -CHARACTER(255) :: ContainerName -LOGICAL :: CoupledPowerPotentialExists -REAL :: TmpArray2(3),CoupledPowerPotentialHDF5(3) -#endif /*defined(PARTICLES)*/ !=================================================================================================================================== IF((.NOT.InterpolationInitIsDone).OR.EquationInitIsDone)THEN LBWRITE(*,*) "InitPoisson not ready to be called or already called." @@ -140,6 +134,7 @@ SUBROUTINE InitEquation() END IF LBWRITE(UNIT_StdOut,'(132("-"))') LBWRITE(UNIT_stdOut,'(A)') ' INIT POISSON...' +IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) ! Read in boundary parameters IniExactFunc = GETINT('IniExactFunc') @@ -148,13 +143,6 @@ SUBROUTINE InitEquation() nRefStateMax = 0 nLinStateMax = 0 ! Defaults -#if defined(PARTICLES) -#if USE_LOADBALANCE -IF(.NOT.PerformLoadBalance)& -#endif /*USE_LOADBALANCE*/ - CalcPCouplElectricPotential=.FALSE. -#endif /*defined(PARTICLES)*/ -PCouplAlreadySet=.FALSE. DO i=1,nBCs BCType = BoundaryType(i,BC_TYPE) BCState = BoundaryType(i,BC_STATE) @@ -170,57 +158,13 @@ SUBROUTINE InitEquation() nRefStateMax = MAX(nRefStateMax,BCState) ELSEIF(BCType.EQ.7.AND.BCState.GT.0)THEN nLinStateMax = MAX(nLinStateMax,BCState) - ELSEIF(BCType.EQ.2)THEN + END IF +END DO #if defined(PARTICLES) - ! Special BC with that adjusts a floating potential in order to meet a user-specified power input - ! Only for BoundaryType = (/2,2/) - IF(BCState.EQ.2)THEN -#if USE_LOADBALANCE - ! Do not set during load balance in order to keep the old value - IF(PerformLoadBalance) CYCLE -#endif /*USE_LOADBALANCE*/ - ! Do not set twice - IF(PCouplAlreadySet) CYCLE - PCouplAlreadySet=.TRUE. - - ! Electric potential (/min, start, max/) Note that start is only required at t=0 and is used for the BC potential - CoupledPowerPotential = GETREALARRAY('CoupledPowerPotential',3) - - ! Restart: Only root reads state file to prevent access with a large number of processors - IF(DoRestart)THEN - ! Only root reads the values and distributes them via MPI Broadcast - IF(MPIRoot)THEN - CALL OpenDataFile(RestartFile,create=.FALSE.,single=.TRUE.,readOnly=.TRUE.) - ! Check old parameter name - ContainerName='CoupledPowerPotential' - CALL DatasetExists(File_ID,TRIM(ContainerName),CoupledPowerPotentialExists) - ! Check for new parameter name - IF(CoupledPowerPotentialExists)THEN - CALL ReadArray(TRIM(ContainerName) , 2 , (/1_IK , 3_IK/) , 0_IK , 1 , RealArray=TmpArray2) - CoupledPowerPotentialHDF5 = TmpArray2 - WRITE(UNIT_stdOut,'(3(A,ES10.2E3))') " Read CoupledPowerPotential from restart file ["//TRIM(RestartFile)//"] min[V]: ",& - CoupledPowerPotentialHDF5(1),", current[V]: ",CoupledPowerPotentialHDF5(2),", max[V]: ",CoupledPowerPotentialHDF5(3) - END IF ! CoupledPowerPotentialExists - CALL CloseDataFile() - END IF ! MPIRoot -#if USE_MPI - ! Broadcast from root to other processors - CALL MPI_BCAST(CoupledPowerPotentialHDF5,3, MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,iERROR) -#endif /*USE_MPI*/ - END IF ! DoRestart - - ! Over-write during restart with values from hdf5 state file - IF(DoRestart) CoupledPowerPotential = CoupledPowerPotentialHDF5 - - CoupledPowerTarget = GETREAL('CoupledPowerTarget') - IF(CoupledPowerTarget.LE.0.) CALL abort(__STAMP__,'CoupledPowerTarget must be > 0.') - CalcPCouplElectricPotential = .TRUE. - CoupledPowerRelaxFac = GETREAL('CoupledPowerRelaxFac') - END IF ! BCState.EQ.1 +! Coupled Power Potential: Adjust the electric potential on a BC to meet a specific power absorbed by the charged particles +CALL InitCoupledPowerPotential() #endif /*defined(PARTICLES)*/ - END IF -END DO ! Read linear potential boundaries nLinState = CountOption('LinPhi') @@ -312,6 +256,106 @@ SUBROUTINE InitEquation() END SUBROUTINE InitEquation +#if defined(PARTICLES) +!=================================================================================================================================== +!> Initialize containers for coupled power potential for adjusting the electric potential on a BC to meet a specific power absorbed +!> by the charged particles +!=================================================================================================================================== +SUBROUTINE InitCoupledPowerPotential() +! MODULES +! insert modules here +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! +! INPUT / OUTPUT VARIABLES +! Space-separated list of input and output types. Use: (int|real|logical|...)_(in|out|inout)_dim(n) +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER, PARAMETER :: BCTypeCPP(1:2) = (/2,52/) +#if defined(PARTICLES) +CHARACTER(255) :: ContainerName +LOGICAL :: CoupledPowerPotentialExists +REAL :: TmpArray2(3),CoupledPowerPotentialHDF5(3) +#endif /*defined(PARTICLES)*/ +INTEGER :: CPPBoundaries +!=================================================================================================================================== +! Get global number of coupled power potential boundaries in [1:nBCs] +CalcPCouplElectricPotential=.FALSE. +CPPBoundaries = 0 +DO iBC=1,nBCs + BCType = BoundaryType(iBC,BC_TYPE) + IF(.NOT.ANY(BCType.EQ.BCTypeCPP)) CYCLE ! Skip non-CPP boundaries + BCState = BoundaryType(iBC,BC_STATE) ! State is iFPC + IF((BCType.EQ.2).AND.(BCState.NE.2)) CYCLE ! BCType=2 must be combined with BCState=2 + CPPBoundaries=CPPBoundaries+1 + IF(BCState.LE.0) CALL CollectiveStop(__STAMP__,' BCState for FPC must be >0! BCState=',IntInfo=BCState) +END DO + +IF(CPPBoundaries.EQ.0) RETURN ! Already determined in HDG initialization + +CalcPCouplElectricPotential = .TRUE. + +#if USE_LOADBALANCE +! Do not set during load balance in order to keep the old value +IF(.NOT.PerformLoadBalance)THEN +#endif /*USE_LOADBALANCE*/ + ! Electric potential (/min, start, max/) Note that start is only required at t=0 and is used for the BC potential + CoupledPowerPotential = GETREALARRAY('CoupledPowerPotential',3) + ! Get target power value + CoupledPowerTarget = GETREAL('CoupledPowerTarget') + IF(CoupledPowerTarget.LE.0.) CALL abort(__STAMP__,'CoupledPowerTarget must be > 0.') + ! Get relaxation factor + CoupledPowerRelaxFac = GETREAL('CoupledPowerRelaxFac') + IF(CoupledPowerRelaxFac.LE.0.) CALL abort(__STAMP__,'CoupledPowerRelaxFac must be > 0.') +#if USE_LOADBALANCE +END IF ! .NOT.PerformLoadBalance +#endif /*USE_LOADBALANCE*/ + +! When restarting, load the last known values of CoupledPowerPotential(1:3) from the .h5 state file +CALL ReadCPPDataFromH5() + + + + + + ! Restart: Only root reads state file to prevent access with a large number of processors + IF(DoRestart)THEN + ! Only root reads the values and distributes them via MPI Broadcast + IF(MPIRoot)THEN + CALL OpenDataFile(RestartFile,create=.FALSE.,single=.TRUE.,readOnly=.TRUE.) + ! Check old parameter name + ContainerName='CoupledPowerPotential' + CALL DatasetExists(File_ID,TRIM(ContainerName),CoupledPowerPotentialExists) + ! Check for new parameter name + IF(CoupledPowerPotentialExists)THEN + CALL ReadArray(TRIM(ContainerName) , 2 , (/1_IK , 3_IK/) , 0_IK , 1 , RealArray=TmpArray2) + CoupledPowerPotentialHDF5 = TmpArray2 + WRITE(UNIT_stdOut,'(3(A,ES10.2E3))') " Read CoupledPowerPotential from restart file ["//TRIM(RestartFile)//"] min[V]: ",& + CoupledPowerPotentialHDF5(1),", current[V]: ",CoupledPowerPotentialHDF5(2),", max[V]: ",CoupledPowerPotentialHDF5(3) + END IF ! CoupledPowerPotentialExists + CALL CloseDataFile() + END IF ! MPIRoot +#if USE_MPI + ! Broadcast from root to other processors + CALL MPI_BCAST(CoupledPowerPotentialHDF5,3, MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,iERROR) +#endif /*USE_MPI*/ + END IF ! DoRestart + + ! Over-write during restart with values from hdf5 state file + IF(DoRestart) CoupledPowerPotential = CoupledPowerPotentialHDF5 + + END IF ! BCState.EQ.1 + + + + + + + + +END SUBROUTINE InitCoupledPowerPotential +#endif /*defined(PARTICLES)*/ + + SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState,iLinState,BCState) !=================================================================================================================================== ! Specifies all the initial conditions. The state in conservative variables is returned. @@ -806,7 +850,6 @@ PPURE SUBROUTINE CalcSourceHDG(i,j,k,iElem,resu, Phi, warning_linear, warning_li END SUBROUTINE CalcSourceHDG - FUNCTION shapefunc(r) !=================================================================================================================================== ! Implementation of (possibly several different) shapefunctions diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index 34da9f4f3..014fce67e 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -246,7 +246,7 @@ SUBROUTINE InitHDG() BCType =BoundaryType(BC(SideID),BC_TYPE) BCState=BoundaryType(BC(SideID),BC_STATE) SELECT CASE(BCType) - CASE(2,4,5,6,7,8) ! Dirichlet + CASE(2,4,5,6,7,8,50) ! Dirichlet nDirichletBCsides=nDirichletBCsides+1 CASE(10,11) ! Neumann nNeumannBCsides=nNeumannBCsides+1 @@ -263,6 +263,14 @@ SUBROUTINE InitHDG() ! Conductor: Initialize electric potential condition (resistive decharging of surface and electric potential calculation) CALL InitEPC() +! Bias Voltage: Initialize containers and sub-communicator +! BCType: 50,X for bias voltage +! BCType: 51,X for bias voltage + cos(wt) function +! BCType: 52,X for bias voltage + cos(wt) function + coupled power adjustment (for AC and not DC in this case) +!CALL InitBV() +write(*,*) "2.) CALL InitBV()" +IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) + ! Check if zero potential must be set on a boundary (or periodic side) CALL InitZeroPotential() @@ -282,7 +290,7 @@ SUBROUTINE InitHDG() BCType =BoundaryType(BC(SideID),BC_TYPE) BCState=BoundaryType(BC(SideID),BC_STATE) SELECT CASE(BCType) - CASE(2,4,5,6,7,8) ! Dirichlet + CASE(2,4,5,6,7,8,50) ! Dirichlet nDirichletBCsides=nDirichletBCsides+1 DirichletBC(nDirichletBCsides)=SideID MaskedSide(SideID)=1 @@ -558,7 +566,7 @@ SUBROUTINE InitZeroPotential() SELECT CASE(BCType) CASE(1) ! periodic ! do nothing - CASE(2,4,5,6,7,8) ! Dirichlet + CASE(2,4,5,6,7,8,50) ! Dirichlet ZeroPotentialSideID = 0 ! no zero potential required EXIT ! as soon as one Dirichlet BC is found, no zero potential must be used CASE(10,11) ! Neumann @@ -677,6 +685,7 @@ SUBROUTINE InitFPC() !----------------------------------------------------------------------------------------------------------------------------------! !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES +INTEGER, PARAMETER :: BCTypeFPC = 20 INTEGER :: BCType,BCState,iUniqueFPCBC INTEGER :: SideID,iBC #if USE_MPI @@ -699,7 +708,7 @@ SUBROUTINE InitFPC() FPC%nUniqueFPCBounds = 0 DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.20) CYCLE ! Skip non-FPC boundaries + IF(BCType.NE.BCTypeFPC) CYCLE ! Skip non-FPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! State is iFPC FPC%nFPCBounds=FPC%nFPCBounds+1 IF(BCState.LE.0) CALL CollectiveStop(__STAMP__,' BCState for FPC must be >0! BCState=',IntInfo=BCState) @@ -716,7 +725,7 @@ SUBROUTINE InitFPC() ! contribute to the total floating boundary condition. If yes, then this processor is part of the communicator DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.20) CYCLE ! Skip non-FPC boundaries + IF(BCType.NE.BCTypeFPC) CYCLE ! Skip non-FPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! State is iFPC WRITE(UNIT=hilf,FMT='(I0)') BCState WRITE(UNIT=hilf2,FMT='(I0)') FPC%nFPCBounds @@ -737,7 +746,7 @@ SUBROUTINE InitFPC() FPC%BCState = 0 DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.20) CYCLE + IF(BCType.NE.BCTypeFPC) CYCLE BCState = BoundaryType(iBC,BC_STATE) ! State is iFPC iUniqueFPCBC = FPC%Group(BCState,2) FPC%BCState(iUniqueFPCBC) = BCState @@ -773,7 +782,7 @@ SUBROUTINE InitFPC() DO SideID=1,nBCSides iBC = BC(SideID) BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.20) CYCLE ! Skip non-FPC boundaries + IF(BCType.NE.BCTypeFPC) CYCLE ! Skip non-FPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iFPC FPC%Group(BCState,3) = FPC%Group(BCState,3) + 1 END DO ! SideID=1,nBCSides @@ -790,7 +799,7 @@ SUBROUTINE InitFPC() DO SideID=1,nBCSides iBC = BC(SideID) BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.20) CYCLE ! Skip non-FPC boundaries + IF(BCType.NE.BCTypeFPC) CYCLE ! Skip non-FPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iFPC iUniqueFPCBC = FPC%Group(BCState,2) BConProc(iUniqueFPCBC) = .TRUE. @@ -829,7 +838,7 @@ SUBROUTINE InitFPC() ! Get boundary type BCType = BoundaryType(BCIndex,BC_TYPE) ! Check if FPC has been found - IF(BCType.EQ.20)THEN + IF(BCType.EQ.BCTypeFPC)THEN ! Check if the BC can be reached iGlobElemCenter(1:3) = (/ SUM(BoundsOfElem_Shared(1:2,1,iGlobElem)),& @@ -851,7 +860,7 @@ SUBROUTINE InitFPC() ! Go to next element CYCLE iCNElemLoop END IF ! VECNORM( ... - END IF ! BCType.EQ.20 + END IF ! BCType.EQ.BCTypeFPC END IF ! BCIndex.GT.0 END DO ! iSide = ElemInfo_Shared(ELEM_FIRSTSIDEIND,iGlobElem)+1,ElemInfo_Shared(ELEM_LASTSIDEIND,iGlobElem) END DO iCNElemLoop ! iCNElem = 1,nComputeNodeTotalElems @@ -927,15 +936,295 @@ END SUBROUTINE InitFPC !> subsequently an electric potential is created. !> !> 1.) Loop over all field BCs and check if the current processor is either the MPI root or has at least one of the BCs that -!> contribute to the total floating boundary condition. If yes, then this processor is part of the communicator -!> 2.) Create Mapping from floating boundary condition BC index to field BC index -!> 3.) Create Mapping from field BC index to floating boundary condition BC index +!> contribute to the total electric potential boundary condition. If yes, then this processor is part of the communicator +!> 2.) Create Mapping from electric potential boundary condition BC index to field BC index +!> 3.) Check if field BC is on current proc (or MPI root) +!> 3.1) Each processor loops over all of his elements +!> 3.2) Loop over all compute-node elements (every processor loops over all of these elements) +!> 4.) Create MPI sub-communicators +!=================================================================================================================================== +SUBROUTINE InitEPC() +! MODULES +USE MOD_Globals ! ,ONLY: MPIRoot,iError,myrank,UNIT_stdOut,MPI_COMM_WORLD +USE MOD_Preproc +USE MOD_Mesh_Vars ,ONLY: nBCs,BoundaryType +USE MOD_Analyze_Vars ,ONLY: DoFieldAnalyze +USE MOD_HDG_Vars ,ONLY: UseEPC,EPC +USE MOD_Mesh_Vars ,ONLY: nBCSides,BC +#if USE_LOADBALANCE +USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance +#endif /*USE_LOADBALANCE*/ +#if USE_MPI && defined(PARTICLES) +USE MOD_Mesh_Tools ,ONLY: GetGlobalElemID +USE MOD_Globals ,ONLY: ElementOnProc +USE MOD_Particle_Mesh_Vars ,ONLY: ElemInfo_Shared,BoundsOfElem_Shared,SideInfo_Shared +USE MOD_MPI_Shared_Vars ,ONLY: nComputeNodeTotalElems +USE MOD_Particle_MPI_Vars ,ONLY: halo_eps +USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem +#endif /*USE_MPI && defined(PARTICLES)*/ +USE MOD_ReadInTools ,ONLY: GETREALARRAY +IMPLICIT NONE +! INPUT / OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------! +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +INTEGER, PARAMETER :: BCTypeEPC = 8 +INTEGER :: BCType,BCState,iUniqueEPCBC +INTEGER :: SideID,iBC +#if USE_MPI +INTEGER :: color,WithSides +LOGICAL,ALLOCATABLE :: BConProc(:) +#if defined(PARTICLES) +INTEGER :: iElem,iCNElem +REAL :: iElemCenter(1:3),iGlobElemCenter(1:3) +REAL :: iElemRadius,iGlobElemRadius +INTEGER :: iGlobElem,BCIndex,iSide +#endif /*defined(PARTICLES)*/ +#endif /*USE_MPI*/ +CHARACTER(5) :: hilf,hilf2 +!=================================================================================================================================== + +! Get global number of EPC boundaries in [1:nBCs], they might belong to the same group (will be reduced to "nUniqueEPCBounds" below) +! EPC boundaries with the same BCState will be in the same group (electrically connected) +UseEPC = .FALSE. +EPC%nEPCBounds = 0 +EPC%nUniqueEPCBounds = 0 +DO iBC=1,nBCs + BCType = BoundaryType(iBC,BC_TYPE) + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries + BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC + EPC%nEPCBounds=EPC%nEPCBounds+1 + IF(BCState.LE.0) CALL CollectiveStop(__STAMP__,' BCState for EPC must be >0! BCState=',IntInfo=BCState) +END DO + +IF(EPC%nEPCBounds.EQ.0) RETURN ! Already determined in HDG initialization + +UseEPC = .TRUE. + +ALLOCATE(EPC%Group(1:EPC%nEPCBounds,3)) +EPC%Group = 0 ! Initialize + +! 1.) Loop over all field BCs and check if the current processor is either the MPI root or has at least one of the BCs that +! contribute to the total electric potential boundary condition. If yes, then this processor is part of the communicator +DO iBC=1,nBCs + BCType = BoundaryType(iBC,BC_TYPE) + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries + BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC + WRITE(UNIT=hilf,FMT='(I0)') BCState + WRITE(UNIT=hilf2,FMT='(I0)') EPC%nEPCBounds + IF(BCState.GT.EPC%nEPCBounds) CALL CollectiveStop(__STAMP__,& + 'BCState='//TRIM(hilf)//' must be smaller or equal than the total number of '//TRIM(hilf2)//' EPC boundaries!') + EPC%Group(BCState,1) = EPC%Group(BCState,1) + 1 + IF(EPC%Group(BCState,1).EQ.1) THEN + EPC%nUniqueEPCBounds = EPC%nUniqueEPCBounds +1 ! Only count once + EPC%Group(BCState,2) = EPC%nUniqueEPCBounds + END IF +END DO + +! Automatically activate surface model analyze flag +DoFieldAnalyze = .TRUE. + +! Read resistances for each unique EPC +EPC%Resistance = GETREALARRAY('EPC-Resistance',EPC%nUniqueEPCBounds) + +! 2.) Create Mapping from electric potential boundary condition BC index to BCState +ALLOCATE(EPC%BCState(EPC%nUniqueEPCBounds)) +EPC%BCState = 0 +DO iBC=1,nBCs + BCType = BoundaryType(iBC,BC_TYPE) + IF(BCType.NE.BCTypeEPC) CYCLE + BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC + iUniqueEPCBC = EPC%Group(BCState,2) + EPC%BCState(iUniqueEPCBC) = BCState +END DO + +! Allocate the containers +! This container is not deallocated for MPIRoot when performing load balance (only root needs this info to write it to .csv) +IF(.NOT.ALLOCATED(EPC%Voltage))THEN + ALLOCATE(EPC%Voltage(1:EPC%nUniqueEPCBounds)) + EPC%Voltage = 0. +END IF ! .NOT.ALLOCATED(EPC%Voltage) +ALLOCATE(EPC%VoltageProc(1:EPC%nUniqueEPCBounds)) +EPC%VoltageProc = 0. +! This container is not deallocated for MPIRoot when performing load balance as this process updates the information on the new +! sub-communicator processes during load balance +IF(.NOT.ALLOCATED(EPC%Charge))THEN + ALLOCATE(EPC%Charge(1:EPC%nUniqueEPCBounds)) + EPC%Charge = 0. +END IF ! .NOT.ALLOCATED(EPC%Charge) +ALLOCATE(EPC%ChargeProc(1:EPC%nUniqueEPCBounds)) +EPC%ChargeProc = 0. + +! Get processor-local number of EPC sides associated with each i-th EPC boundary +! Check local sides +DO SideID=1,nBCSides + iBC = BC(SideID) + BCType = BoundaryType(iBC,BC_TYPE) + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries + BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iEPC + EPC%Group(BCState,3) = EPC%Group(BCState,3) + 1 +END DO ! SideID=1,nBCSides + +#if USE_MPI +! 3) Check if field BC is on current proc (or MPI root) +ALLOCATE(BConProc(EPC%nUniqueEPCBounds)) +BConProc = .FALSE. +IF(MPIRoot)THEN + BConProc = .TRUE. +ELSE + + ! Check local sides + DO SideID=1,nBCSides + iBC = BC(SideID) + BCType = BoundaryType(iBC,BC_TYPE) + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries + BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iEPC + iUniqueEPCBC = EPC%Group(BCState,2) + BConProc(iUniqueEPCBC) = .TRUE. + END DO ! SideID=1,nBCSides + +#if defined(PARTICLES) + ! Check if all FPCs have already been found + IF(.NOT.(ALL(BConProc)))THEN + ! Particles might impact the EPC on another proc/node. Therefore check if a particle can travel from a local element to an + ! element that has at least one side, which is an EPC + ! 3.1) Each processor loops over all of his elements + iElemLoop: DO iElem = 1+offsetElem, nElems+offsetElem + + iElemCenter(1:3) = (/ SUM(BoundsOfElem_Shared(1:2,1,iElem)),& + SUM(BoundsOfElem_Shared(1:2,2,iElem)),& + SUM(BoundsOfElem_Shared(1:2,3,iElem)) /) / 2. + iElemRadius = VECNORM ((/ BoundsOfElem_Shared(2,1,iElem)-BoundsOfElem_Shared(1,1,iElem),& + BoundsOfElem_Shared(2,2,iElem)-BoundsOfElem_Shared(1,2,iElem),& + BoundsOfElem_Shared(2,3,iElem)-BoundsOfElem_Shared(1,3,iElem) /) / 2.) + + ! 3.2) Loop over all compute-node elements (every processor loops over all of these elements) + ! Loop ALL compute-node elements (use global element index) + iCNElemLoop: DO iCNElem = 1,nComputeNodeTotalElems + iGlobElem = GetGlobalElemID(iCNElem) + + ! Skip my own elements as they have already been tested when the local sides are checked + IF(ElementOnProc(iGlobElem)) CYCLE iCNElemLoop + + ! Check if one of the six sides of the compute-node element is a EPC + ! Note that iSide is in the range of 1:nNonUniqueGlobalSides + DO iSide = ElemInfo_Shared(ELEM_FIRSTSIDEIND,iGlobElem)+1,ElemInfo_Shared(ELEM_LASTSIDEIND,iGlobElem) + ! Get BC index of the global side index + BCIndex = SideInfo_Shared(SIDE_BCID,iSide) + ! Only check BC sides with BC index > 0 + IF(BCIndex.GT.0)THEN + ! Get boundary type + BCType = BoundaryType(BCIndex,BC_TYPE) + ! Check if EPC has been found + IF(BCType.EQ.BCTypeEPC)THEN + + ! Check if the BC can be reached + iGlobElemCenter(1:3) = (/ SUM(BoundsOfElem_Shared(1:2,1,iGlobElem)),& + SUM(BoundsOfElem_Shared(1:2,2,iGlobElem)),& + SUM(BoundsOfElem_Shared(1:2,3,iGlobElem)) /) / 2. + iGlobElemRadius = VECNORM ((/ BoundsOfElem_Shared(2,1,iGlobElem)-BoundsOfElem_Shared(1,1,iGlobElem),& + BoundsOfElem_Shared(2,2,iGlobElem)-BoundsOfElem_Shared(1,2,iGlobElem),& + BoundsOfElem_Shared(2,3,iGlobElem)-BoundsOfElem_Shared(1,3,iGlobElem) /) / 2.) + + ! check if compute-node element "iGlobElem" is within halo_eps of processor-local element "iElem" + IF (VECNORM( iElemCenter(1:3) - iGlobElemCenter(1:3) ) .LE. ( halo_eps + iElemRadius + iGlobElemRadius ) )THEN + BCState = BoundaryType(BCIndex,BC_STATE) ! BCState corresponds to iEPC + IF(BCState.LT.1) CALL abort(__STAMP__,'BCState cannot be <1',IntInfoOpt=BCState) + iUniqueEPCBC = EPC%Group(BCState,2) + ! Flag the i-th EPC + BConProc(iUniqueEPCBC) = .TRUE. + ! Check if all FPCs have been found -> exit complete loop + IF(ALL(BConProc)) EXIT iElemLoop + ! Go to next element + CYCLE iCNElemLoop + END IF ! VECNORM( ... + END IF ! BCType.EQ.BCTypeEPC + END IF ! BCIndex.GT.0 + END DO ! iSide = ElemInfo_Shared(ELEM_FIRSTSIDEIND,iGlobElem)+1,ElemInfo_Shared(ELEM_LASTSIDEIND,iGlobElem) + END DO iCNElemLoop ! iCNElem = 1,nComputeNodeTotalElems + END DO iElemLoop ! iElem = 1, nElems + END IF ! .NOT.(ALL(BConProc)) +#endif /*defined(PARTICLES)*/ + +END IF ! MPIRoot + +! 4.) Create MPI sub-communicators +ALLOCATE(EPC%COMM(EPC%nUniqueEPCBounds)) +DO iUniqueEPCBC = 1, EPC%nUniqueEPCBounds + ! Create new communicator + color = MERGE(iUniqueEPCBC, MPI_UNDEFINED, BConProc(iUniqueEPCBC)) + + ! Set communicator id + EPC%COMM(iUniqueEPCBC)%ID=iUniqueEPCBC + + ! Create new emission communicator for Electric potential boundary condition communication. + ! Pass MPI_INFO_NULL as rank to follow the original ordering + CALL MPI_COMM_SPLIT(MPI_COMM_WORLD, color, MPI_INFO_NULL, EPC%COMM(iUniqueEPCBC)%UNICATOR, iError) + + ! Find my rank on the shared communicator, comm size and proc name + IF(BConProc(iUniqueEPCBC))THEN + CALL MPI_COMM_RANK(EPC%COMM(iUniqueEPCBC)%UNICATOR, EPC%COMM(iUniqueEPCBC)%MyRank, iError) + CALL MPI_COMM_SIZE(EPC%COMM(iUniqueEPCBC)%UNICATOR, EPC%COMM(iUniqueEPCBC)%nProcs, iError) + + ! inform about size of emission communicator + IF (EPC%COMM(iUniqueEPCBC)%MyRank.EQ.0) THEN +#if USE_LOADBALANCE + IF(.NOT.PerformLoadBalance)& +#endif /*USE_LOADBALANCE*/ + WRITE(UNIT_StdOut,'(A,I0,A,I0,A,I0)') ' Electric potential boundary condition: Emission-Communicator ',iUniqueEPCBC,' on ',& + EPC%COMM(iUniqueEPCBC)%nProcs,' procs for BCState ',EPC%BCState(iUniqueEPCBC) + END IF + END IF ! BConProc(iUniqueEPCBC) +END DO ! iUniqueEPCBC = 1, EPC%nUniqueEPCBounds +DEALLOCATE(BConProc) + +! Get the number of procs that actually have a local BC side that is an EPC (required for voltage output to .csv) +! Procs might have zero EPC sides but are in the group because 1.) MPIRoot or 2.) the EPC is in the halo region +! Because only the MPI root process writes the .csv data, the information regarding the voltage on each EPC must be +! communicated with this process even though it might not be connected to each EPC boundary +DO iUniqueEPCBC = 1, EPC%nUniqueEPCBounds + ASSOCIATE( COMM => EPC%COMM(iUniqueEPCBC)%UNICATOR, nProcsWithSides => EPC%COMM(iUniqueEPCBC)%nProcsWithSides ) + IF(EPC%COMM(iUniqueEPCBC)%UNICATOR.NE.MPI_COMM_NULL)THEN + ! Check if the current processor is actually connected to the EPC via a BC side + IF(EPC%Group(EPC%BCState(iUniqueEPCBC),3).EQ.0)THEN + WithSides = 0 + ELSE + WithSides = 1 + END IF ! EPC%Group(EPC%BCState(iUniqueEPCBC),3).EQ.0 + ! Calculate the sum across the sub-communicator. Only the MPI root process needs this information + IF(MPIRoot)THEN + CALL MPI_REDUCE(WithSides, nProcsWithSides, 1 ,MPI_INTEGER, MPI_SUM, 0, COMM, iError) + ! Sanity check + IF(nProcsWithSides.EQ.0) CALL abort(__STAMP__,'Found EPC with no processors connected to it') + ELSE + CALL MPI_REDUCE(WithSides, 0 , 1 ,MPI_INTEGER, MPI_SUM, 0, COMM, IError) + END IF ! MPIRoot + END IF ! EPC%COMM(iUniqueEPCBC)%UNICATOR.NE.MPI_COMM_NULL + END ASSOCIATE +END DO ! iUniqueEPCBC = 1, EPC%nUniqueEPCBounds +#endif /*USE_MPI*/ + +! When restarting, load the deposited charge on each EPC from the .h5 state file +CALL ReadEPCDataFromH5() + +END SUBROUTINE InitEPC + + +!=================================================================================================================================== +!> Create containers and communicators for each bias-voltage boundary condition where impacting charges are removed and +!> subsequently an electric potential is created (the particle communication is part of the BPO analysis and required here). +!> +!> 1.) Loop over all field BCs and check if the current processor is either the MPI root or has at least one of the BCs that +!> contribute to the total bias-voltage boundary condition. If yes, then this processor is part of the communicator + +!> 2.) Create Mapping from electric potential boundary condition BC index to field BC index +!> 3.) Create Mapping from field BC index to electric potential boundary condition BC index !> 4.0) Check if field BC is on current proc (or MPI root) !> 4.1.) Each processor loops over all of his elements !> 4.2.) Loop over all compute-node elements (every processor loops over all of these elements) !> 5.) Create MPI sub-communicators !=================================================================================================================================== -SUBROUTINE InitEPC() +SUBROUTINE InitBV() ! MODULES USE MOD_Globals ! ,ONLY: MPIRoot,iError,myrank,UNIT_stdOut,MPI_COMM_WORLD USE MOD_Preproc @@ -960,6 +1249,7 @@ SUBROUTINE InitEPC() !----------------------------------------------------------------------------------------------------------------------------------! !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES +INTEGER, PARAMETER :: BCTypeEPC = 8 INTEGER :: BCType,BCState,iUniqueEPCBC INTEGER :: SideID,iBC #if USE_MPI @@ -982,7 +1272,7 @@ SUBROUTINE InitEPC() EPC%nUniqueEPCBounds = 0 DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.8) CYCLE ! Skip non-EPC boundaries + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC EPC%nEPCBounds=EPC%nEPCBounds+1 IF(BCState.LE.0) CALL CollectiveStop(__STAMP__,' BCState for EPC must be >0! BCState=',IntInfo=BCState) @@ -996,10 +1286,10 @@ SUBROUTINE InitEPC() EPC%Group = 0 ! Initialize ! 1.) Loop over all field BCs and check if the current processor is either the MPI root or has at least one of the BCs that -! contribute to the total floating boundary condition. If yes, then this processor is part of the communicator +! contribute to the total electric potential boundary condition. If yes, then this processor is part of the communicator DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.8) CYCLE ! Skip non-EPC boundaries + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC WRITE(UNIT=hilf,FMT='(I0)') BCState WRITE(UNIT=hilf2,FMT='(I0)') EPC%nEPCBounds @@ -1018,12 +1308,12 @@ SUBROUTINE InitEPC() ! Read resistances for each unique EPC EPC%Resistance = GETREALARRAY('EPC-Resistance',EPC%nUniqueEPCBounds) -! 2.) Create Mapping from floating boundary condition BC index to BCState +! 2.) Create Mapping from electric potential boundary condition BC index to BCState ALLOCATE(EPC%BCState(EPC%nUniqueEPCBounds)) EPC%BCState = 0 DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.8) CYCLE + IF(BCType.NE.BCTypeEPC) CYCLE BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC iUniqueEPCBC = EPC%Group(BCState,2) EPC%BCState(iUniqueEPCBC) = BCState @@ -1046,7 +1336,7 @@ SUBROUTINE InitEPC() ALLOCATE(EPC%ChargeProc(1:EPC%nUniqueEPCBounds)) EPC%ChargeProc = 0. -!! 3.) Create Mapping from field BC index to floating boundary condition BC index +!! 3.) Create Mapping from field BC index to electric potential boundary condition BC index !ALLOCATE(EPC%BCIDToEPCBCID(nBCs)) !EPC%BCIDToEPCBCID = -1 !DO iEPCBC = 1, EPC%NBoundaries @@ -1059,7 +1349,7 @@ SUBROUTINE InitEPC() DO SideID=1,nBCSides iBC = BC(SideID) BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.8) CYCLE ! Skip non-EPC boundaries + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iEPC EPC%Group(BCState,3) = EPC%Group(BCState,3) + 1 END DO ! SideID=1,nBCSides @@ -1076,7 +1366,7 @@ SUBROUTINE InitEPC() DO SideID=1,nBCSides iBC = BC(SideID) BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.8) CYCLE ! Skip non-EPC boundaries + IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iEPC iUniqueEPCBC = EPC%Group(BCState,2) BConProc(iUniqueEPCBC) = .TRUE. @@ -1115,7 +1405,7 @@ SUBROUTINE InitEPC() ! Get boundary type BCType = BoundaryType(BCIndex,BC_TYPE) ! Check if EPC has been found - IF(BCType.EQ.8)THEN + IF(BCType.EQ.BCTypeEPC)THEN ! Check if the BC can be reached iGlobElemCenter(1:3) = (/ SUM(BoundsOfElem_Shared(1:2,1,iGlobElem)),& @@ -1137,7 +1427,7 @@ SUBROUTINE InitEPC() ! Go to next element CYCLE iCNElemLoop END IF ! VECNORM( ... - END IF ! BCType.EQ.8 + END IF ! BCType.EQ.BCTypeEPC END IF ! BCIndex.GT.0 END DO ! iSide = ElemInfo_Shared(ELEM_FIRSTSIDEIND,iGlobElem)+1,ElemInfo_Shared(ELEM_LASTSIDEIND,iGlobElem) END DO iCNElemLoop ! iCNElem = 1,nComputeNodeTotalElems @@ -1156,7 +1446,7 @@ SUBROUTINE InitEPC() ! set communicator id EPC%COMM(iUniqueEPCBC)%ID=iUniqueEPCBC - ! create new emission communicator for floating boundary condition communication. Pass MPI_INFO_NULL as rank to follow the original ordering + ! create new emission communicator for electric potential boundary condition communication. Pass MPI_INFO_NULL as rank to follow the original ordering CALL MPI_COMM_SPLIT(MPI_COMM_WORLD, color, MPI_INFO_NULL, EPC%COMM(iUniqueEPCBC)%UNICATOR, iError) ! Find my rank on the shared communicator, comm size and proc name @@ -1169,7 +1459,7 @@ SUBROUTINE InitEPC() #if USE_LOADBALANCE IF(.NOT.PerformLoadBalance)& #endif /*USE_LOADBALANCE*/ - WRITE(UNIT_StdOut,'(A,I0,A,I0,A,I0)') ' Floating boundary condition: Emission-Communicator ',iUniqueEPCBC,' on ',& + WRITE(UNIT_StdOut,'(A,I0,A,I0,A,I0)') ' electric potential boundary condition: Emission-Communicator ',iUniqueEPCBC,' on ',& EPC%COMM(iUniqueEPCBC)%nProcs,' procs for BCState ',EPC%BCState(iUniqueEPCBC) END IF END IF ! BConProc(iUniqueEPCBC) @@ -1205,7 +1495,7 @@ SUBROUTINE InitEPC() ! When restarting, load the deposited charge on each EPC from the .h5 state file CALL ReadEPCDataFromH5() -END SUBROUTINE InitEPC +END SUBROUTINE InitBV !=================================================================================================================================== diff --git a/src/mesh/prepare_mesh.f90 b/src/mesh/prepare_mesh.f90 index 9b1cd8e23..2715f4300 100644 --- a/src/mesh/prepare_mesh.f90 +++ b/src/mesh/prepare_mesh.f90 @@ -1041,7 +1041,7 @@ SUBROUTINE fillMeshInfo() SELECT CASE(BCType) CASE(1) !periodic CALL abort(__STAMP__,'SideID.LE.nBCSides and SideID is periodic should not happen') - CASE(2,4,5,6,7,8) ! Dirichlet + CASE(2,4,5,6,7,8,50) ! Dirichlet ! do not consider this side CASE(10,11) ! Neumann HDGSides = HDGSides + 1 From 99e53065c292eb199cf7ff522b1d8e5fbf451055 Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Tue, 18 Apr 2023 14:31:40 +0200 Subject: [PATCH 02/11] Re-structured coupled power potential (CPP) to use sub-cumminicator for broadcasts from MPIRoot to other processes that are connected CPP boundaries. --- src/equations/poisson/equation.f90 | 206 +++++++++++++----- src/equations/poisson/equation_vars.f90 | 7 - src/hdg/hdg.f90 | 4 +- src/hdg/hdg_vars.f90 | 18 +- src/io_hdf5/hdf5_output_state.f90 | 2 +- src/particles/analyze/particle_analyze.f90 | 4 +- .../analyze/particle_analyze_tools.f90 | 2 +- src/particles/particle_init.f90 | 2 +- 8 files changed, 176 insertions(+), 69 deletions(-) diff --git a/src/equations/poisson/equation.f90 b/src/equations/poisson/equation.f90 index c8bbd8c87..6ef2e007f 100644 --- a/src/equations/poisson/equation.f90 +++ b/src/equations/poisson/equation.f90 @@ -107,11 +107,6 @@ SUBROUTINE InitEquation() #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ -#if defined(PARTICLES) -USE MOD_IO_HDF5 ,ONLY: OpenDataFile,CloseDataFile,File_ID -USE MOD_Restart_Vars ,ONLY: DoRestart,RestartFile -USE MOD_HDF5_Input ,ONLY: DatasetExists,ReadArray -#endif /*defined(PARTICLES)*/ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -134,7 +129,7 @@ SUBROUTINE InitEquation() END IF LBWRITE(UNIT_StdOut,'(132("-"))') LBWRITE(UNIT_stdOut,'(A)') ' INIT POISSON...' -IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) +!IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) ! Read in boundary parameters IniExactFunc = GETINT('IniExactFunc') @@ -258,33 +253,49 @@ END SUBROUTINE InitEquation #if defined(PARTICLES) !=================================================================================================================================== -!> Initialize containers for coupled power potential for adjusting the electric potential on a BC to meet a specific power absorbed -!> by the charged particles +!> Initialize containers for coupled power potential (CPP) for adjusting the electric potential on a BC to meet a specific power +!> absorbed by the charged particles +!> 1.) Get global number of coupled power potential boundaries in [1:nBCs] +!> 2.) Get parameters +!> 3.) Establish sub-communicator (all BCs directly connected to a coupled power potential boundary); MPIRoot is always in the group) +!> 4.) When restarting, load the last known values of CoupledPowerPotential(1:3) from the .h5 state file !=================================================================================================================================== SUBROUTINE InitCoupledPowerPotential() ! MODULES -! insert modules here +#if USE_LOADBALANCE +USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance +#endif /*USE_LOADBALANCE*/ +USE MOD_Globals ,ONLY: CollectiveStop +USE MOD_HDG_Vars ,ONLY: CalcPCouplElectricPotential,CoupledPowerPotential,CoupledPowerTarget,CoupledPowerRelaxFac +USE MOD_ReadInTools ,ONLY: GETREALARRAY,GETREAL,GETINT,CountOption +USE MOD_Mesh_Vars ,ONLY: BoundaryType,nBCs +#if USE_MPI +USE MOD_Globals ,ONLY: IERROR,MPI_COMM_NULL,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,MPI_INFO_NULL,MPI_UNDEFINED,MPIRoot +USE MOD_Globals ,ONLY: UNIT_StdOut +USE MOD_HDG_Vars ,ONLY: CPPCOMM +USE MOD_Mesh_Vars ,ONLY: nBCSides,BC +#endif /*USE_MPI*/ IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------! ! INPUT / OUTPUT VARIABLES -! Space-separated list of input and output types. Use: (int|real|logical|...)_(in|out|inout)_dim(n) !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER, PARAMETER :: BCTypeCPP(1:2) = (/2,52/) -#if defined(PARTICLES) -CHARACTER(255) :: ContainerName -LOGICAL :: CoupledPowerPotentialExists -REAL :: TmpArray2(3),CoupledPowerPotentialHDF5(3) -#endif /*defined(PARTICLES)*/ -INTEGER :: CPPBoundaries +INTEGER, PARAMETER :: BCTypeCPP(1:2) = (/2,52/) ! BCType which allows coupled power potential control +INTEGER :: iBC,CPPBoundaries,BCType,BCState +#if USE_MPI +LOGICAL :: BConProc +INTEGER :: color,SideID +#endif /*USE_MPI*/ !=================================================================================================================================== -! Get global number of coupled power potential boundaries in [1:nBCs] -CalcPCouplElectricPotential=.FALSE. + +! 1.) Get global number of coupled power potential boundaries in [1:nBCs] +CalcPCouplElectricPotential = .FALSE. CPPBoundaries = 0 +! Loop over global BC list DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) IF(.NOT.ANY(BCType.EQ.BCTypeCPP)) CYCLE ! Skip non-CPP boundaries - BCState = BoundaryType(iBC,BC_STATE) ! State is iFPC + BCState = BoundaryType(iBC,BC_STATE) ! BCState of the corresponding BCType IF((BCType.EQ.2).AND.(BCState.NE.2)) CYCLE ! BCType=2 must be combined with BCState=2 CPPBoundaries=CPPBoundaries+1 IF(BCState.LE.0) CALL CollectiveStop(__STAMP__,' BCState for FPC must be >0! BCState=',IntInfo=BCState) @@ -294,6 +305,7 @@ SUBROUTINE InitCoupledPowerPotential() CalcPCouplElectricPotential = .TRUE. +! 2.) Get parameters #if USE_LOADBALANCE ! Do not set during load balance in order to keep the old value IF(.NOT.PerformLoadBalance)THEN @@ -302,57 +314,143 @@ SUBROUTINE InitCoupledPowerPotential() CoupledPowerPotential = GETREALARRAY('CoupledPowerPotential',3) ! Get target power value CoupledPowerTarget = GETREAL('CoupledPowerTarget') - IF(CoupledPowerTarget.LE.0.) CALL abort(__STAMP__,'CoupledPowerTarget must be > 0.') + IF(CoupledPowerTarget.LE.0.) CALL CollectiveStop(__STAMP__,'CoupledPowerTarget must be > 0.') ! Get relaxation factor CoupledPowerRelaxFac = GETREAL('CoupledPowerRelaxFac') - IF(CoupledPowerRelaxFac.LE.0.) CALL abort(__STAMP__,'CoupledPowerRelaxFac must be > 0.') + IF(CoupledPowerRelaxFac.LE.0.) CALL CollectiveStop(__STAMP__,'CoupledPowerRelaxFac must be > 0.') #if USE_LOADBALANCE +ELSE + CoupledPowerPotential(1:3) = 0. END IF ! .NOT.PerformLoadBalance #endif /*USE_LOADBALANCE*/ -! When restarting, load the last known values of CoupledPowerPotential(1:3) from the .h5 state file -CALL ReadCPPDataFromH5() - - - - - - ! Restart: Only root reads state file to prevent access with a large number of processors - IF(DoRestart)THEN - ! Only root reads the values and distributes them via MPI Broadcast - IF(MPIRoot)THEN - CALL OpenDataFile(RestartFile,create=.FALSE.,single=.TRUE.,readOnly=.TRUE.) - ! Check old parameter name - ContainerName='CoupledPowerPotential' - CALL DatasetExists(File_ID,TRIM(ContainerName),CoupledPowerPotentialExists) - ! Check for new parameter name - IF(CoupledPowerPotentialExists)THEN - CALL ReadArray(TRIM(ContainerName) , 2 , (/1_IK , 3_IK/) , 0_IK , 1 , RealArray=TmpArray2) - CoupledPowerPotentialHDF5 = TmpArray2 - WRITE(UNIT_stdOut,'(3(A,ES10.2E3))') " Read CoupledPowerPotential from restart file ["//TRIM(RestartFile)//"] min[V]: ",& - CoupledPowerPotentialHDF5(1),", current[V]: ",CoupledPowerPotentialHDF5(2),", max[V]: ",CoupledPowerPotentialHDF5(3) - END IF ! CoupledPowerPotentialExists - CALL CloseDataFile() - END IF ! MPIRoot #if USE_MPI - ! Broadcast from root to other processors - CALL MPI_BCAST(CoupledPowerPotentialHDF5,3, MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,iERROR) +! 3.) Establish sub-communicator (all BCs directly connected to a coupled power potential boundary); MPIRoot is always in the group +BConProc = .FALSE. +IF(MPIRoot)THEN + BConProc = .TRUE. +ELSE + ! Check local sides + DO SideID=1,nBCSides + iBC = BC(SideID) + BCType = BoundaryType(iBC,BC_TYPE) + IF(.NOT.ANY(BCType.EQ.BCTypeCPP)) CYCLE ! Skip other boundaries + BCState = BoundaryType(iBC,BC_STATE) ! BCState of the corresponding BCType + IF((BCType.EQ.2).AND.(BCState.NE.2)) CYCLE ! BCType=2 must be combined with BCState=2 + BConProc = .TRUE. + END DO ! SideID=1,nBCSides +END IF ! MPIRoot + +! create new communicator +color = MERGE(CPPBoundaries, MPI_UNDEFINED, BConProc) + +! set communicator id +CPPCOMM%ID = CPPBoundaries + +! create new emission communicator for electric potential boundary condition communication. Pass MPI_INFO_NULL as rank to follow the original ordering +CALL MPI_COMM_SPLIT(MPI_COMM_WORLD, color, MPI_INFO_NULL, CPPCOMM%UNICATOR, iError) + +! Find my rank on the shared communicator, comm size and proc name +IF(BConProc)THEN + CALL MPI_COMM_RANK(CPPCOMM%UNICATOR, CPPCOMM%MyRank, iError) + CALL MPI_COMM_SIZE(CPPCOMM%UNICATOR, CPPCOMM%nProcs, iError) + + ! inform about size of emission communicator + IF (CPPCOMM%MyRank.EQ.0) THEN +#if USE_LOADBALANCE + IF(.NOT.PerformLoadBalance)& +#endif /*USE_LOADBALANCE*/ + WRITE(UNIT_StdOut,'(A,I0,A)') ' Coupled power potential boundary condition: Emission-Communicator on ',& + CPPCOMM%nProcs,' procs' + END IF +END IF ! BConProc #endif /*USE_MPI*/ - END IF ! DoRestart - ! Over-write during restart with values from hdf5 state file - IF(DoRestart) CoupledPowerPotential = CoupledPowerPotentialHDF5 +! 4.) When restarting, load the last known values of CoupledPowerPotential(1:3) from the .h5 state file +CALL ReadCPPDataFromH5() - END IF ! BCState.EQ.1 +END SUBROUTINE InitCoupledPowerPotential +!=================================================================================================================================== +!> Read the coupled power potential (CPP) data from a .h5 state file. +!> 1. The MPI root process reads the info and checks data consistency +!> 2. The MPI root process distributes the information among the sub-communicator processes connected to the CPP boundary. +!=================================================================================================================================== +SUBROUTINE ReadCPPDataFromH5() +! MODULES +USE MOD_io_hdf5 +USE MOD_Globals ,ONLY: UNIT_stdOut,MPIRoot,IK,abort +#if USE_LOADBALANCE +USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance,UseH5IOLoadBalance +#endif /*USE_LOADBALANCE*/ +USE MOD_IO_HDF5 ,ONLY: OpenDataFile,CloseDataFile,File_ID +USE MOD_Restart_Vars ,ONLY: DoRestart,RestartFile +USE MOD_HDF5_Input ,ONLY: DatasetExists,ReadArray,GetDataSize +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! +! INPUT / OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +CHARACTER(255) :: ContainerName +LOGICAL :: CPPExists +REAL :: CPPDataHDF5(1:3) +!=================================================================================================================================== +! Only required during restart +IF(.NOT.DoRestart) RETURN +#if USE_LOADBALANCE +! Do not try to read the data from .h5 if load balance is performed without creating a .h5 restart file +IF(PerformLoadBalance.AND..NOT.(UseH5IOLoadBalance)) RETURN +#endif /*USE_LOADBALANCE*/ +! 1. The MPI root process reads the info and checks data consistency +! Only root reads the values and distributes them via MPI Broadcast +IF(MPIRoot)THEN + CALL OpenDataFile(RestartFile,create=.FALSE.,single=.TRUE.,readOnly=.TRUE.) + ! Check old parameter name + ContainerName='CoupledPowerPotential' + CALL DatasetExists(File_ID,TRIM(ContainerName),CPPExists) + ! Check for new parameter name + IF(CPPExists)THEN + CALL ReadArray(TRIM(ContainerName) , 2 , (/1_IK , 3_IK/) , 0_IK , 1 , RealArray=CPPDataHDF5) + WRITE(UNIT_stdOut,'(3(A,ES10.2E3))') " Read CoupledPowerPotential from restart file ["//TRIM(RestartFile)//"] min[V]: ",& + CPPDataHDF5(1),", current[V]: ",CPPDataHDF5(2),", max[V]: ",CPPDataHDF5(3) + CoupledPowerPotential = CPPDataHDF5 + END IF ! CPPExists + CALL CloseDataFile() +END IF ! MPIRoot +#if USE_MPI +! 2. The MPI root process distributes the information among the sub-communicator processes for each EPC +CALL SynchronizeCPP() +#endif /*USE_MPI*/ +END SUBROUTINE ReadCPPDataFromH5 -END SUBROUTINE InitCoupledPowerPotential +#if USE_MPI +!=================================================================================================================================== +!> Communicate the CPP values from MPIRoot to sub-communicator processes +!=================================================================================================================================== +SUBROUTINE SynchronizeCPP() +! MODULES +USE MOD_Globals ,ONLY: IERROR,MPI_COMM_NULL,MPI_DOUBLE_PRECISION +USE MOD_HDG_Vars ,ONLY: CPPCOMM +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! +! INPUT / OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +!=================================================================================================================================== +IF(CPPCOMM%UNICATOR.NE.MPI_COMM_NULL)THEN + ! Broadcast from root to other processors on the sub-communicator + CALL MPI_BCAST(CoupledPowerPotential, 3, MPI_DOUBLE_PRECISION, 0, CPPCOMM%UNICATOR, IERROR) +END IF +END SUBROUTINE SynchronizeCPP +#endif /*USE_MPI*/ #endif /*defined(PARTICLES)*/ @@ -365,7 +463,7 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState,iLinState,BCState) USE MOD_Globals_Vars ,ONLY: PI,ElementaryCharge,eps0 USE MOD_Equation_Vars ,ONLY: IniCenter,IniHalfwidth,IniAmplitude,RefState,LinPhi,LinPhiHeight,LinPhiNormal,LinPhiBasePoint #if defined(PARTICLES) -USE MOD_Equation_Vars ,ONLY: CoupledPowerPotential +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential USE MOD_Particle_Vars ,ONLY: Species,nSpecies #endif /*defined(PARTICLES)*/ USE MOD_Dielectric_Vars ,ONLY: DielectricRatio,Dielectric_E_0,DielectricRadiusValue,DielectricEpsR diff --git a/src/equations/poisson/equation_vars.f90 b/src/equations/poisson/equation_vars.f90 index 83c395547..0147a1a06 100644 --- a/src/equations/poisson/equation_vars.f90 +++ b/src/equations/poisson/equation_vars.f90 @@ -67,12 +67,5 @@ MODULE MOD_Equation_Vars REAL,ALLOCATABLE :: LinPhiHeight(:) REAL,ALLOCATABLE :: LinPhi(:) -#if defined(PARTICLES) -! Special BC with floating potential that is defined by the absorbed power of the charged particles -REAL :: CoupledPowerPotential(3) ! (/min, start, max/) electric potential at all BoundaryType = (/2,2/) -REAL :: CoupledPowerTarget ! Target input power at all BoundaryType = (/2,2/) -REAL :: CoupledPowerRelaxFac ! Relaxation factor for calculation of new electric potential -LOGICAL :: CalcPCouplElectricPotential! Switch calculation on/off -#endif /*defined(PARTICLES)*/ !=================================================================================================================================== END MODULE MOD_Equation_Vars diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index 014fce67e..113c8af3a 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -268,8 +268,8 @@ SUBROUTINE InitHDG() ! BCType: 51,X for bias voltage + cos(wt) function ! BCType: 52,X for bias voltage + cos(wt) function + coupled power adjustment (for AC and not DC in this case) !CALL InitBV() -write(*,*) "2.) CALL InitBV()" -IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) +!write(*,*) "2.) CALL InitBV()" +!IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) ! Check if zero potential must be set on a boundary (or periodic side) CALL InitZeroPotential() diff --git a/src/hdg/hdg_vars.f90 b/src/hdg/hdg_vars.f90 index c0fa03481..5a11d1793 100644 --- a/src/hdg/hdg_vars.f90 +++ b/src/hdg/hdg_vars.f90 @@ -203,7 +203,7 @@ MODULE MOD_HDG_Vars TYPE(tFPC) :: FPC !=================================================================================================================================== -!-- Electric Potential Condition (for decharging ) +!-- Electric Potential Condition (for decharging) !=================================================================================================================================== LOGICAL :: UseEPC !< Automatic flag when EPCs are active @@ -232,6 +232,22 @@ MODULE MOD_HDG_Vars END TYPE TYPE(tEPC) :: EPC +#if defined(PARTICLES) +!=================================================================================================================================== +!-- Coupled Power Potential (CPP) +!-- Special BC with floating potential that is defined by the absorbed power of the charged particles +!=================================================================================================================================== + +LOGICAL :: CalcPCouplElectricPotential! Switch calculation on/off + +#if USE_MPI +TYPE(tMPIGROUP) :: CPPCOMM !< communicator and ID for parallel execution +#endif /*USE_MPI*/ + +REAL :: CoupledPowerPotential(3) !< (/min, start, max/) electric potential at all BoundaryType = (/2,2/) +REAL :: CoupledPowerTarget !< Target input power at all BoundaryType = (/2,2/) +REAL :: CoupledPowerRelaxFac !< Relaxation factor for calculation of new electric potential +#endif /*defined(PARTICLES)*/ !=================================================================================================================================== #if USE_MPI diff --git a/src/io_hdf5/hdf5_output_state.f90 b/src/io_hdf5/hdf5_output_state.f90 index 9194dbc24..ec948f4a3 100644 --- a/src/io_hdf5/hdf5_output_state.f90 +++ b/src/io_hdf5/hdf5_output_state.f90 @@ -101,7 +101,7 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) USE MOD_Particle_Analyze_Tools ,ONLY: AllocateElectronIonDensityCell,AllocateElectronTemperatureCell USE MOD_Particle_Analyze_Tools ,ONLY: CalculateElectronIonDensityCell,CalculateElectronTemperatureCell USE MOD_HDF5_Output_Particles ,ONLY: AddBRElectronFluidToPartSource -USE MOD_Equation_Vars ,ONLY: CoupledPowerPotential,CalcPCouplElectricPotential +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,CalcPCouplElectricPotential #endif /*PARTICLES*/ #endif /*USE_HDG*/ USE MOD_Analyze_Vars ,ONLY: OutputTimeFixed diff --git a/src/particles/analyze/particle_analyze.f90 b/src/particles/analyze/particle_analyze.f90 index f4b110536..97a309536 100644 --- a/src/particles/analyze/particle_analyze.f90 +++ b/src/particles/analyze/particle_analyze.f90 @@ -138,7 +138,7 @@ SUBROUTINE InitParticleAnalyze() USE MOD_Particle_Analyze_Tools,ONLY: AllocateElectronIonDensityCell,AllocateElectronTemperatureCell,AllocateCalcElectronEnergy #if USE_HDG USE MOD_HDG_Vars ,ONLY: CalcBRVariableElectronTemp,BRAutomaticElectronRef -USE MOD_Equation_Vars ,ONLY: CalcPCouplElectricPotential +USE MOD_HDG_Vars ,ONLY: CalcPCouplElectricPotential #endif /*USE_HDG*/ USE MOD_Particle_Analyze_Tools,ONLY: CalcNumberDensityBGGasDistri #if USE_LOADBALANCE @@ -852,7 +852,7 @@ SUBROUTINE AnalyzeParticles(Time) #if USE_HDG USE MOD_HDG_Vars ,ONLY: BRNbrOfRegions,CalcBRVariableElectronTemp,BRAutomaticElectronRef,RegionElectronRef USE MOD_Globals_Vars ,ONLY: BoltzmannConst,ElementaryCharge -USE MOD_Equation_Vars ,ONLY: CalcPCouplElectricPotential,CoupledPowerPotential +USE MOD_HDG_Vars ,ONLY: CalcPCouplElectricPotential,CoupledPowerPotential USE MOD_Particle_Analyze_Tools ,ONLY: CalculatePCouplElectricPotential #endif /*USE_HDG*/ USE MOD_Globals_Vars ,ONLY: eV2Kelvin diff --git a/src/particles/analyze/particle_analyze_tools.f90 b/src/particles/analyze/particle_analyze_tools.f90 index 879800ca8..62bc3b295 100644 --- a/src/particles/analyze/particle_analyze_tools.f90 +++ b/src/particles/analyze/particle_analyze_tools.f90 @@ -3094,7 +3094,7 @@ END SUBROUTINE CalcCoupledPowerPart SUBROUTINE CalculatePCouplElectricPotential() ! MODULES USE MOD_Globals -USE MOD_Equation_Vars ,ONLY: CoupledPowerPotential,CoupledPowerTarget,CoupledPowerRelaxFac +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,CoupledPowerTarget,CoupledPowerRelaxFac USE MOD_Particle_Analyze_Vars ,ONLY: PCoupl IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------! diff --git a/src/particles/particle_init.f90 b/src/particles/particle_init.f90 index a6ae504b8..ed1a162be 100644 --- a/src/particles/particle_init.f90 +++ b/src/particles/particle_init.f90 @@ -412,7 +412,7 @@ SUBROUTINE InitializeVariables() USE MOD_TimeDisc_Vars ,ONLY: ManualTimeStep,useManualTimeStep #if defined(PARTICLES) && USE_HDG USE MOD_Part_BR_Elecron_Fluid ,ONLY: InitializeVariablesElectronFluidRegions -USE MOD_Equation_Vars ,ONLY: CalcPCouplElectricPotential +USE MOD_HDG_Vars ,ONLY: CalcPCouplElectricPotential #endif /*defined(PARTICLES) && USE_HDG*/ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE From 541ae9f912d6c6bf96506c24c36ced8bf594398f Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 19 Apr 2023 00:35:28 +0200 Subject: [PATCH 03/11] Added structure for bias voltage with load balance synchronization, restart functionality and output to SurfaceAnalyze.csv. --- src/analyze/analyzefield.f90 | 9 +- src/equations/poisson/equation.f90 | 23 +- src/hdg/hdg.f90 | 454 ++++++++---------- src/hdg/hdg_vars.f90 | 32 +- src/io_hdf5/hdf5_output_state.f90 | 20 +- src/mesh/prepare_mesh.f90 | 2 +- .../surfacemodel/surfacemodel_analyze.f90 | 96 +++- src/piclas.h | 5 + src/restart/restart_field.f90 | 12 +- 9 files changed, 372 insertions(+), 281 deletions(-) diff --git a/src/analyze/analyzefield.f90 b/src/analyze/analyzefield.f90 index 7486743ee..d155a8899 100644 --- a/src/analyze/analyzefield.f90 +++ b/src/analyze/analyzefield.f90 @@ -63,6 +63,7 @@ SUBROUTINE AnalyzeField(Time) USE MOD_Preproc USE MOD_Analyze_Vars ,ONLY: DoFieldAnalyze,CalcEpot,WEl USE MOD_Analyze_Vars ,ONLY: CalcBoundaryFieldOutput,BFO +USE MOD_Mesh_Vars ,ONLY: BoundaryName #if (PP_nVar>=6) USE MOD_Analyze_Vars ,ONLY: CalcPoyntingInt,nPoyntingIntPlanes,PosPoyntingInt #endif /*PP_nVar>=6*/ @@ -75,7 +76,6 @@ SUBROUTINE AnalyzeField(Time) #if USE_HDG USE MOD_HDG_Vars ,ONLY: HDGNorm,iterationTotal,RunTimeTotal,UseFPC,FPC,UseEPC,EPC USE MOD_Analyze_Vars ,ONLY: AverageElectricPotential,CalcAverageElectricPotential,EDC,CalcElectricTimeDerivative -USE MOD_Mesh_Vars ,ONLY: BoundaryName USE MOD_TimeDisc_Vars ,ONLY: dt #endif /*USE_HDG*/ #ifdef PARTICLES @@ -246,8 +246,9 @@ SUBROUTINE AnalyzeField(Time) IF(CalcBoundaryFieldOutput)THEN DO iBoundary=1,BFO%NFieldBoundaries nOutputVarTotal = nOutputVarTotal + 1 - WRITE(StrVarNameTmp,'(A,I0.3)') 'BFO-boundary',iBoundary - WRITE(tmpStr(nOutputVarTotal),'(A,I0.3,A)')delimiter//'"',nOutputVarTotal,'-'//TRIM(StrVarNameTmp)//'"' + WRITE(StrVarNameTmp,'(A,I0.3)') 'BFO-boundary-',iBoundary + WRITE(tmpStr(nOutputVarTotal),'(A,I0.3,A)')delimiter//'"',nOutputVarTotal,'-'//TRIM(StrVarNameTmp)//'-'//& + TRIM(BoundaryName(BFO%FieldBoundaries(iBoundary)))//'"' END DO END IF @@ -1837,6 +1838,8 @@ SUBROUTINE CalculateBoundaryFieldOutput(iBC,Time,BoundaryFieldOutput) CALL ExactFunc( -1 , x , BoundaryFieldOutput , t=Time , iRefState=BCState) CASE(6) ! exact BC = Dirichlet BC !! ExactFunc via RefState (Time is optional) CALL ExactFunc( -2 , x , BoundaryFieldOutput , t=time , iRefState=BCState) + CASE(50) ! exact BC = Dirichlet BC !! ExactFunc via bias voltage DC + CALL ExactFunc( -5 , x , BoundaryFieldOutput ) END SELECT ! BCType END ASSOCIATE #else diff --git a/src/equations/poisson/equation.f90 b/src/equations/poisson/equation.f90 index 6ef2e007f..4c61d4d40 100644 --- a/src/equations/poisson/equation.f90 +++ b/src/equations/poisson/equation.f90 @@ -43,8 +43,16 @@ MODULE MOD_Equation INTERFACE FinalizeEquation MODULE PROCEDURE FinalizeEquation END INTERFACE +#if USE_MPI +INTERFACE SynchronizeCPP + MODULE PROCEDURE SynchronizeCPP +END INTERFACE +#endif /*USE_MPI*/ -PUBLIC::InitEquation,ExactFunc,CalcSource,FinalizeEquation, CalcSourceHDG,DivCleaningDamping +PUBLIC :: InitEquation,ExactFunc,CalcSource,FinalizeEquation, CalcSourceHDG,DivCleaningDamping +#if USE_MPI +PUBLIC :: SynchronizeCPP +#endif /*USE_MPI*/ !=================================================================================================================================== PUBLIC::DefineParametersEquation CONTAINS @@ -129,7 +137,6 @@ SUBROUTINE InitEquation() END IF LBWRITE(UNIT_StdOut,'(132("-"))') LBWRITE(UNIT_stdOut,'(A)') ' INIT POISSON...' -!IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) ! Read in boundary parameters IniExactFunc = GETINT('IniExactFunc') @@ -255,6 +262,7 @@ END SUBROUTINE InitEquation !=================================================================================================================================== !> Initialize containers for coupled power potential (CPP) for adjusting the electric potential on a BC to meet a specific power !> absorbed by the charged particles +!> !> 1.) Get global number of coupled power potential boundaries in [1:nBCs] !> 2.) Get parameters !> 3.) Establish sub-communicator (all BCs directly connected to a coupled power potential boundary); MPIRoot is always in the group) @@ -281,6 +289,8 @@ SUBROUTINE InitCoupledPowerPotential() !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES INTEGER, PARAMETER :: BCTypeCPP(1:2) = (/2,52/) ! BCType which allows coupled power potential control +! ! 2: DC potential adjustment +! ! 52: bias voltage + cos(wt) function + coupled power for AC potential adjustment INTEGER :: iBC,CPPBoundaries,BCType,BCState #if USE_MPI LOGICAL :: BConProc @@ -319,8 +329,6 @@ SUBROUTINE InitCoupledPowerPotential() CoupledPowerRelaxFac = GETREAL('CoupledPowerRelaxFac') IF(CoupledPowerRelaxFac.LE.0.) CALL CollectiveStop(__STAMP__,'CoupledPowerRelaxFac must be > 0.') #if USE_LOADBALANCE -ELSE - CoupledPowerPotential(1:3) = 0. END IF ! .NOT.PerformLoadBalance #endif /*USE_LOADBALANCE*/ @@ -330,6 +338,7 @@ SUBROUTINE InitCoupledPowerPotential() IF(MPIRoot)THEN BConProc = .TRUE. ELSE + CoupledPowerPotential(1:3) = 0. ! Check local sides DO SideID=1,nBCSides iBC = BC(SideID) @@ -463,7 +472,7 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState,iLinState,BCState) USE MOD_Globals_Vars ,ONLY: PI,ElementaryCharge,eps0 USE MOD_Equation_Vars ,ONLY: IniCenter,IniHalfwidth,IniAmplitude,RefState,LinPhi,LinPhiHeight,LinPhiNormal,LinPhiBasePoint #if defined(PARTICLES) -USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,BiasVoltage USE MOD_Particle_Vars ,ONLY: Species,nSpecies #endif /*defined(PARTICLES)*/ USE MOD_Dielectric_Vars ,ONLY: DielectricRatio,Dielectric_E_0,DielectricRadiusValue,DielectricEpsR @@ -488,6 +497,10 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState,iLinState,BCState) REAL :: Omega,r1,r2,r_2D,r_3D,r_bary,cos_theta,eps1,eps2,xi,a(3),b(3),Q !=================================================================================================================================== SELECT CASE (ExactFunction) +#if defined(PARTICLES) +CASE(-5) ! Bias voltage DC boundary + Resu(:) = BiasVoltage%BVData(1) +#endif /*defined(PARTICLES)*/ CASE(-4) ! Electric potential condition (EPC) where charge accumulated over one time step is removed and creates a voltage Resu(:) = EPC%Voltage(EPC%Group(BCState,2)) CASE(-3) ! linear function with base point, normal vector and heigh: Requires BoundaryType = (/7,X/) diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index 113c8af3a..ed1282c95 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -40,7 +40,7 @@ MODULE MOD_HDG PUBLIC :: HDG, RestartHDG PUBLIC :: DefineParametersHDG #if USE_MPI -PUBLIC :: SynchronizeChargeOnFPC,SynchronizeVoltageOnEPC +PUBLIC :: SynchronizeChargeOnFPC,SynchronizeVoltageOnEPC,SynchronizeBV #endif /*USE_MPI*/ #endif /*USE_HDG*/ !=================================================================================================================================== @@ -68,26 +68,25 @@ SUBROUTINE DefineParametersHDG() CALL prms%CreateIntOption( 'AdaptIterNewtonToLinear','Maximum number of iterations where the exact source derivative is used before it is switched to the linearization', '100') CALL prms%CreateIntOption( 'MaxIterNewton' ,'Maximum number of iterations in the Newton solver', '10000') CALL prms%CreateRealOption( 'EpsNonLinear' ,'Abort residual of the Newton solver', '1.0E-6') -CALL prms%CreateIntOption( 'PrecondType' ,'Preconditioner type\n 0: no preconditioner\n 1: Side-block SPD preconditioner'& - //' matrix (already Cholesky decomposed)\n 2: Inverse of diagonal preconditioned', '2') +CALL prms%CreateIntOption( 'PrecondType' ,'Preconditioner type\n 0: no preconditioner\n 1: Side-block SPD preconditioner matrix (already Cholesky decomposed)\n 2: Inverse of diagonal preconditioned', '2') CALL prms%CreateRealOption( 'epsCG' ,'Abort residual of the CG solver', '1.0E-6') CALL prms%CreateIntOption( 'OutIterCG' ,'Number of iteration steps between output of CG solver info to std out', '1') - CALL prms%CreateLogicalOption('useRelativeAbortCrit' ,'Switch between relative and absolute abort criterion', '.FALSE.') CALL prms%CreateIntOption( 'MaxIterCG' ,'Maximum number of iterations in the CG solver', '500') CALL prms%CreateLogicalOption('ExactLambda' ,'Initially set lambda on all sides (volume and boundaries) via a pre-defined function (ExactFunc)', '.FALSE.') - CALL prms%CreateIntOption( 'HDGSkip' ,'Number of time step iterations until the HDG solver is called (i.e. all intermediate calls are skipped)', '0') -CALL prms%CreateIntOption( 'HDGSkipInit' ,'Number of time step iterations until the HDG solver is called (i.e. all intermediate calls are skipped)'& - //' while time < HDGSkip_t0 (if HDGSkip > 0)', '0') +CALL prms%CreateIntOption( 'HDGSkipInit' ,'Number of time step iterations until the HDG solver is called (i.e. all intermediate calls are skipped) while time < HDGSkip_t0 (if HDGSkip > 0)', '0') CALL prms%CreateRealOption( 'HDGSkip_t0' ,'Time during which HDGSkipInit is used instead of HDGSkip (if HDGSkip > 0)', '0.') - CALL prms%CreateLogicalOption('HDGDisplayConvergence' ,'Display divergence criteria: Iterations, RunTime and Residual', '.FALSE.') - -CALL prms%CreateIntOption( 'HDGZeroPotentialDir' ,'Direction in which a Dirichlet condition with phi=0 is superimposed on the boundary conditions'& - //' (1: x, 2: y, 3: z). The default chooses the direction automatically when no other Dirichlet boundary conditions are defined.'& - ,'-1') +CALL prms%CreateIntOption( 'HDGZeroPotentialDir' ,'Direction in which a Dirichlet condition with phi=0 is superimposed on the boundary conditions (1: x, 2: y, 3: z). The default chooses the direction automatically when no other Dirichlet boundary conditions are defined.','-1') CALL prms%CreateRealArrayOption( 'EPC-Resistance' , 'Vector (length corresponds to the number of EPC boundaries) with the resistance for each EPC in Ohm', no=0) +#if defined(PARTICLES) +CALL prms%CreateLogicalOption( 'UseBiasVoltage' , 'Activate usage of bias voltage adjustment (for specific boundaries only)', '.FALSE.') +CALL prms%CreateIntOption( 'BiasVoltage-NPartBoundaries' , 'Number of particle boundaries where the total ion excess is to be calculated for bias voltage model') +CALL prms%CreateIntArrayOption( 'Biasvoltage-PartBoundaries' , 'Particle boundary index of boundaries where the total ion excess is to be calculated for bias voltage model', no=0) +CALL prms%CreateRealOption( 'BiasVoltage-Frequency' , 'Frequency of the sinusoidal field boundary where the bias voltage is applied (a value of 0.0 corresponds to a DC potential BC). The total particle electric current over one cycle is required to converge to zero.') +CALL prms%CreateRealOption( 'BiasVoltage-Delta' , 'Bias voltage difference used for adjusting the DC voltage of the corresponding BC') +#endif /*defined(PARTICLES)*/ ! --- BR electron fluid #if defined(PARTICLES) @@ -246,7 +245,7 @@ SUBROUTINE InitHDG() BCType =BoundaryType(BC(SideID),BC_TYPE) BCState=BoundaryType(BC(SideID),BC_STATE) SELECT CASE(BCType) - CASE(2,4,5,6,7,8,50) ! Dirichlet + CASE(HDGDIRICHLETBCSIDEIDS) ! Dirichlet nDirichletBCsides=nDirichletBCsides+1 CASE(10,11) ! Neumann nNeumannBCsides=nNeumannBCsides+1 @@ -263,13 +262,13 @@ SUBROUTINE InitHDG() ! Conductor: Initialize electric potential condition (resistive decharging of surface and electric potential calculation) CALL InitEPC() +#if defined(PARTICLES) ! Bias Voltage: Initialize containers and sub-communicator -! BCType: 50,X for bias voltage -! BCType: 51,X for bias voltage + cos(wt) function +! BCType: 50,X for bias voltage + DC boundary condition +! BCType: 51,X for bias voltage + cos(wt) function boundary condition ! BCType: 52,X for bias voltage + cos(wt) function + coupled power adjustment (for AC and not DC in this case) -!CALL InitBV() -!write(*,*) "2.) CALL InitBV()" -!IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) +CALL InitBV() +#endif /*defined(PARTICLES)*/ ! Check if zero potential must be set on a boundary (or periodic side) CALL InitZeroPotential() @@ -290,7 +289,7 @@ SUBROUTINE InitHDG() BCType =BoundaryType(BC(SideID),BC_TYPE) BCState=BoundaryType(BC(SideID),BC_STATE) SELECT CASE(BCType) - CASE(2,4,5,6,7,8,50) ! Dirichlet + CASE(HDGDIRICHLETBCSIDEIDS) ! Dirichlet nDirichletBCsides=nDirichletBCsides+1 DirichletBC(nDirichletBCsides)=SideID MaskedSide(SideID)=1 @@ -452,7 +451,6 @@ SUBROUTINE InitHDG() ! nAffectedBlockSides = 22 !END IF ! FPC%nFPCBounds !!IPWRITE(UNIT_StdOut,*) "nAffectedBlockSides =", nAffectedBlockSides -!!IF(myrank.eq.0) read*; CALL MPI_BARRIER(MPI_COMM_WORLD,iError) PetscCallA(MatSEQSBAIJSetPreallocation(Smat_petsc,nGP_face,22,PETSC_NULL_INTEGER,ierr)) PetscCallA(MatMPISBAIJSetPreallocation(Smat_petsc,nGP_face,22,PETSC_NULL_INTEGER,22-1,PETSC_NULL_INTEGER,ierr)) PetscCallA(MatZeroEntries(Smat_petsc,ierr)) @@ -566,7 +564,7 @@ SUBROUTINE InitZeroPotential() SELECT CASE(BCType) CASE(1) ! periodic ! do nothing - CASE(2,4,5,6,7,8,50) ! Dirichlet + CASE(HDGDIRICHLETBCSIDEIDS) ! Dirichlet ZeroPotentialSideID = 0 ! no zero potential required EXIT ! as soon as one Dirichlet BC is found, no zero potential must be used CASE(10,11) ! Neumann @@ -1210,292 +1208,233 @@ SUBROUTINE InitEPC() END SUBROUTINE InitEPC +#if defined(PARTICLES) !=================================================================================================================================== !> Create containers and communicators for each bias-voltage boundary condition where impacting charges are removed and !> subsequently an electric potential is created (the particle communication is part of the BPO analysis and required here). !> -!> 1.) Loop over all field BCs and check if the current processor is either the MPI root or has at least one of the BCs that -!> contribute to the total bias-voltage boundary condition. If yes, then this processor is part of the communicator - -!> 2.) Create Mapping from electric potential boundary condition BC index to field BC index -!> 3.) Create Mapping from field BC index to electric potential boundary condition BC index -!> 4.0) Check if field BC is on current proc (or MPI root) -!> 4.1.) Each processor loops over all of his elements -!> 4.2.) Loop over all compute-node elements (every processor loops over all of these elements) -!> 5.) Create MPI sub-communicators +!> 1.) Activate bias voltage and check number of boundaries +!> 2.) Get bias voltage parameters +!> 3.) Check if actual bias voltage BC is on current process (or MPI root) +!> 4.) Create MPI sub-communicators !=================================================================================================================================== SUBROUTINE InitBV() ! MODULES -USE MOD_Globals ! ,ONLY: MPIRoot,iError,myrank,UNIT_stdOut,MPI_COMM_WORLD -USE MOD_Preproc -USE MOD_Mesh_Vars ,ONLY: nBCs,BoundaryType -USE MOD_Analyze_Vars ,ONLY: DoFieldAnalyze -USE MOD_HDG_Vars ,ONLY: UseEPC,EPC -USE MOD_Mesh_Vars ,ONLY: nBCSides,BC +USE MOD_Globals ,ONLY: CollectiveStop +USE MOD_ReadInTools ,ONLY: GETLOGICAL,GETREAL,GETINT,GETINTARRAY +USE MOD_Particle_Boundary_Vars ,ONLY: PartBound +USE MOD_Mesh_Vars ,ONLY: nBCs,BoundaryType,BoundaryName +USE MOD_HDG_Vars ,ONLY: UseBiasVoltage,BiasVoltage +USE MOD_SurfaceModel_Analyze_Vars ,ONLY: CalcBoundaryParticleOutput,BPO #if USE_LOADBALANCE -USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance +USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ -#if USE_MPI && defined(PARTICLES) -USE MOD_Mesh_Tools ,ONLY: GetGlobalElemID -USE MOD_Globals ,ONLY: ElementOnProc -USE MOD_Particle_Mesh_Vars ,ONLY: ElemInfo_Shared,BoundsOfElem_Shared,SideInfo_Shared -USE MOD_MPI_Shared_Vars ,ONLY: nComputeNodeTotalElems -USE MOD_Particle_MPI_Vars ,ONLY: halo_eps -USE MOD_Mesh_Vars ,ONLY: nElems, offsetElem -#endif /*USE_MPI && defined(PARTICLES)*/ -USE MOD_ReadInTools ,ONLY: GETREALARRAY +#if USE_MPI +USE MOD_Globals ,ONLY: IERROR,MPI_COMM_NULL,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,MPI_INFO_NULL,MPI_UNDEFINED,MPIRoot +USE MOD_Globals ,ONLY: UNIT_StdOut +USE MOD_Mesh_Vars ,ONLY: nBCSides,BC +#endif /*USE_MPI*/ IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! ! INPUT / OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------! -!----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -INTEGER, PARAMETER :: BCTypeEPC = 8 -INTEGER :: BCType,BCState,iUniqueEPCBC -INTEGER :: SideID,iBC +INTEGER, PARAMETER :: BCTypeBV(1:3) = (/50,51,52/) ! BCType which allows bias voltage control +! ! 50: pure DC potential +! ! 51: cos(wt) function with DC bias +! ! 52: cos(wt) function with DC bias + coupled power for AC potential adjustment +INTEGER :: BCType,BVBoundaries,BCState,iBoundary +INTEGER :: SideID,iBC,iPBC #if USE_MPI -INTEGER :: color,WithSides -LOGICAL,ALLOCATABLE :: BConProc(:) -#if defined(PARTICLES) -INTEGER :: iElem,iCNElem -REAL :: iElemCenter(1:3),iGlobElemCenter(1:3) -REAL :: iElemRadius,iGlobElemRadius -INTEGER :: iGlobElem,BCIndex,iSide -#endif /*defined(PARTICLES)*/ +INTEGER :: color +LOGICAL :: BConProc #endif /*USE_MPI*/ -CHARACTER(5) :: hilf,hilf2 !=================================================================================================================================== -! Get global number of EPC boundaries in [1:nBCs], they might belong to the same group (will be reduced to "nUniqueEPCBounds" below) -! EPC boundaries with the same BCState will be in the same group (electrically connected) -UseEPC = .FALSE. -EPC%nEPCBounds = 0 -EPC%nUniqueEPCBounds = 0 -DO iBC=1,nBCs - BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries - BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC - EPC%nEPCBounds=EPC%nEPCBounds+1 - IF(BCState.LE.0) CALL CollectiveStop(__STAMP__,' BCState for EPC must be >0! BCState=',IntInfo=BCState) -END DO +!> 1.) Activate bias voltage and check number of boundaries +! Activate model +UseBiasVoltage = GETLOGICAL('UseBiasVoltage') -IF(EPC%nEPCBounds.EQ.0) RETURN ! Already determined in HDG initialization +! Skip the following if bias voltage is not active +IF(.NOT.UseBiasVoltage)THEN + ! Sanity check before returning: bias voltage BCs cannot be used without activating the bias voltage model + ! Count the number of boundaries that allow bias voltage + BVBoundaries = 0 + DO iBC=1,nBCs + BCType = BoundaryType(iBC,BC_TYPE) + IF(.NOT.ANY(BCType.EQ.BCTypeBV)) CYCLE ! Skip other boundaries + BVBoundaries = BVBoundaries + 1 + END DO + IF(BVBoundaries.GT.0) CALL CollectiveStop(__STAMP__,' Bias voltage BCs require UseBiasVoltage=T!') -UseEPC = .TRUE. + ! Exit this subroutine + RETURN +END IF -ALLOCATE(EPC%Group(1:EPC%nEPCBounds,3)) -EPC%Group = 0 ! Initialize +! CalcBoundaryParticleOutput=T and boundaries must be set correctly +IF(.NOT.CalcBoundaryParticleOutput) CALL CollectiveStop(__STAMP__,' UseBiasVoltage=T requires CalcBoundaryParticleOutput=T!') -! 1.) Loop over all field BCs and check if the current processor is either the MPI root or has at least one of the BCs that -! contribute to the total electric potential boundary condition. If yes, then this processor is part of the communicator +! Count the number of boundaries that allow bias voltage +BVBoundaries = 0 DO iBC=1,nBCs BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries - BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC - WRITE(UNIT=hilf,FMT='(I0)') BCState - WRITE(UNIT=hilf2,FMT='(I0)') EPC%nEPCBounds - IF(BCState.GT.EPC%nEPCBounds) CALL CollectiveStop(__STAMP__,& - 'BCState='//TRIM(hilf)//' must be smaller or equal than the total number of '//TRIM(hilf2)//' EPC boundaries!') - EPC%Group(BCState,1) = EPC%Group(BCState,1) + 1 - IF(EPC%Group(BCState,1).EQ.1) THEN - EPC%nUniqueEPCBounds = EPC%nUniqueEPCBounds +1 ! Only count once - EPC%Group(BCState,2) = EPC%nUniqueEPCBounds - END IF + IF(.NOT.ANY(BCType.EQ.BCTypeBV)) CYCLE ! Skip other boundaries + BVBoundaries = BVBoundaries + 1 END DO -! Automatically activate surface model analyze flag -DoFieldAnalyze = .TRUE. +IF(BVBoundaries.NE.1) CALL CollectiveStop(__STAMP__,' UseBiasVoltage=T requires exactly one boundary with this feature!') -! Read resistances for each unique EPC -EPC%Resistance = GETREALARRAY('EPC-Resistance',EPC%nUniqueEPCBounds) - -! 2.) Create Mapping from electric potential boundary condition BC index to BCState -ALLOCATE(EPC%BCState(EPC%nUniqueEPCBounds)) -EPC%BCState = 0 -DO iBC=1,nBCs - BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.BCTypeEPC) CYCLE - BCState = BoundaryType(iBC,BC_STATE) ! State is iEPC - iUniqueEPCBC = EPC%Group(BCState,2) - EPC%BCState(iUniqueEPCBC) = BCState +!> 2.) Get bias voltage parameters +BiasVoltage%NPartBoundaries = GETINT('BiasVoltage-NPartBoundaries') +BiasVoltage%PartBoundaries = GETINTARRAY('Biasvoltage-PartBoundaries',biasvoltage%npartboundaries) +BiasVoltage%Frequency = GETReal('BiasVoltage-Frequency') +BiasVoltage%Delta = GETReal('BiasVoltage-Delta') +#if USE_LOADBALANCE +! Do not set during load balance in order to keep the old value +IF(.NOT.PerformLoadBalance)& +#endif /*USE_LOADBALANCE*/ + BiasVoltage%BVData = 0. + +IF(BiasVoltage%NPartBoundaries.LT.1) CALL CollectiveStop(__STAMP__,' UseBiasVoltage=T requires one or more particle boundaries!') + +DO iBoundary=1,BiasVoltage%NPartBoundaries + iPBC = BiasVoltage%PartBoundaries(iBoundary) + IF(.NOT.ANY(iPBC.EQ.BPO%PartBoundaries(:))) CALL CollectiveStop(__STAMP__,& + 'One of Biasvoltage-PartBoundaries not defined in any BPO-PartBoundaries') + iBC = PartBound%MapToFieldBC(iPBC) + IF(iBC.GT.SIZE(BoundaryName)) CALL CollectiveStop(__STAMP__,'BiasVoltage-PartBoundaries BC index maps to wrong field BCID= ',& + IntInfo=iBC) + BCType = BoundaryType(iBC,BC_TYPE) + BCState = BoundaryType(iBC,BC_STATE) + LBWRITE(UNIT_stdOut,'(A,I0,A,I0,A)') ' Activated bias voltage by collecting currents from ['//TRIM(BoundaryName(iBC))& + //'] with BCType [',BCType,'] and BCState [',BCState,']' END DO -! Allocate the containers -! This container is not deallocated for MPIRoot when performing load balance (only root needs this info to write it to .csv) -IF(.NOT.ALLOCATED(EPC%Voltage))THEN - ALLOCATE(EPC%Voltage(1:EPC%nUniqueEPCBounds)) - EPC%Voltage = 0. -END IF ! .NOT.ALLOCATED(EPC%Voltage) -ALLOCATE(EPC%VoltageProc(1:EPC%nUniqueEPCBounds)) -EPC%VoltageProc = 0. -! This container is not deallocated for MPIRoot when performing load balance as this process updates the information on the new -! sub-communicator processes during load balance -IF(.NOT.ALLOCATED(EPC%Charge))THEN - ALLOCATE(EPC%Charge(1:EPC%nUniqueEPCBounds)) - EPC%Charge = 0. -END IF ! .NOT.ALLOCATED(EPC%Charge) -ALLOCATE(EPC%ChargeProc(1:EPC%nUniqueEPCBounds)) -EPC%ChargeProc = 0. - -!! 3.) Create Mapping from field BC index to electric potential boundary condition BC index -!ALLOCATE(EPC%BCIDToEPCBCID(nBCs)) -!EPC%BCIDToEPCBCID = -1 -!DO iEPCBC = 1, EPC%NBoundaries -! iBC = EDC%FieldBoundaries(iEDCBC) -! EDC%BCIDToEDCBCID(iBC) = iEDCBC -!END DO ! iEDCBC = 1, EDC%NBoundaries - -! Get processor-local number of EPC sides associated with each i-th EPC boundary -! Check local sides -DO SideID=1,nBCSides - iBC = BC(SideID) - BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries - BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iEPC - EPC%Group(BCState,3) = EPC%Group(BCState,3) + 1 -END DO ! SideID=1,nBCSides - #if USE_MPI -! 4.0) Check if field BC is on current proc (or MPI root) -ALLOCATE(BConProc(EPC%nUniqueEPCBounds)) +!> 3.) Check if actual bias voltage BC is on current process (or MPI root) BConProc = .FALSE. IF(MPIRoot)THEN BConProc = .TRUE. ELSE - + BiasVoltage%BVData = 0. ! Check local sides DO SideID=1,nBCSides iBC = BC(SideID) BCType = BoundaryType(iBC,BC_TYPE) - IF(BCType.NE.BCTypeEPC) CYCLE ! Skip non-EPC boundaries - BCState = BoundaryType(iBC,BC_STATE) ! BCState corresponds to iEPC - iUniqueEPCBC = EPC%Group(BCState,2) - BConProc(iUniqueEPCBC) = .TRUE. + IF(.NOT.ANY(BCType.EQ.BCTypeBV)) CYCLE ! Skip other boundaries + BConProc = .TRUE. END DO ! SideID=1,nBCSides +END IF ! MPIRoot -#if defined(PARTICLES) - ! Check if all FPCs have already been found - IF(.NOT.(ALL(BConProc)))THEN - ! Particles might impact the EPC on another proc/node. Therefore check if a particle can travel from a local element to an - ! element that has at least one side, which is an EPC - ! 4.1.) Each processor loops over all of his elements - iElemLoop: DO iElem = 1+offsetElem, nElems+offsetElem - - iElemCenter(1:3) = (/ SUM(BoundsOfElem_Shared(1:2,1,iElem)),& - SUM(BoundsOfElem_Shared(1:2,2,iElem)),& - SUM(BoundsOfElem_Shared(1:2,3,iElem)) /) / 2. - iElemRadius = VECNORM ((/ BoundsOfElem_Shared(2,1,iElem)-BoundsOfElem_Shared(1,1,iElem),& - BoundsOfElem_Shared(2,2,iElem)-BoundsOfElem_Shared(1,2,iElem),& - BoundsOfElem_Shared(2,3,iElem)-BoundsOfElem_Shared(1,3,iElem) /) / 2.) - - ! 4.2.) Loop over all compute-node elements (every processor loops over all of these elements) - ! Loop ALL compute-node elements (use global element index) - iCNElemLoop: DO iCNElem = 1,nComputeNodeTotalElems - iGlobElem = GetGlobalElemID(iCNElem) +! 4.) Create MPI sub-communicators +! Create new communicator +color = MERGE(BVBoundaries, MPI_UNDEFINED, BConProc) - ! Skip my own elements as they have already been tested when the local sides are checked - IF(ElementOnProc(iGlobElem)) CYCLE iCNElemLoop +! Set communicator id +BiasVoltage%COMM%ID = BVBoundaries - ! Check if one of the six sides of the compute-node element is a EPC - ! Note that iSide is in the range of 1:nNonUniqueGlobalSides - DO iSide = ElemInfo_Shared(ELEM_FIRSTSIDEIND,iGlobElem)+1,ElemInfo_Shared(ELEM_LASTSIDEIND,iGlobElem) - ! Get BC index of the global side index - BCIndex = SideInfo_Shared(SIDE_BCID,iSide) - ! Only check BC sides with BC index > 0 - IF(BCIndex.GT.0)THEN - ! Get boundary type - BCType = BoundaryType(BCIndex,BC_TYPE) - ! Check if EPC has been found - IF(BCType.EQ.BCTypeEPC)THEN - - ! Check if the BC can be reached - iGlobElemCenter(1:3) = (/ SUM(BoundsOfElem_Shared(1:2,1,iGlobElem)),& - SUM(BoundsOfElem_Shared(1:2,2,iGlobElem)),& - SUM(BoundsOfElem_Shared(1:2,3,iGlobElem)) /) / 2. - iGlobElemRadius = VECNORM ((/ BoundsOfElem_Shared(2,1,iGlobElem)-BoundsOfElem_Shared(1,1,iGlobElem),& - BoundsOfElem_Shared(2,2,iGlobElem)-BoundsOfElem_Shared(1,2,iGlobElem),& - BoundsOfElem_Shared(2,3,iGlobElem)-BoundsOfElem_Shared(1,3,iGlobElem) /) / 2.) +! Create new emission communicator for electric potential boundary condition communication. Pass MPI_INFO_NULL as rank to follow the original ordering +CALL MPI_COMM_SPLIT(MPI_COMM_WORLD, color, MPI_INFO_NULL, BiasVoltage%COMM%UNICATOR, iError) - ! check if compute-node element "iGlobElem" is within halo_eps of processor-local element "iElem" - IF (VECNORM( iElemCenter(1:3) - iGlobElemCenter(1:3) ) .LE. ( halo_eps + iElemRadius + iGlobElemRadius ) )THEN - BCState = BoundaryType(BCIndex,BC_STATE) ! BCState corresponds to iEPC - IF(BCState.LT.1) CALL abort(__STAMP__,'BCState cannot be <1',IntInfoOpt=BCState) - iUniqueEPCBC = EPC%Group(BCState,2) - ! Flag the i-th EPC - BConProc(iUniqueEPCBC) = .TRUE. - ! Check if all FPCs have been found -> exit complete loop - IF(ALL(BConProc)) EXIT iElemLoop - ! Go to next element - CYCLE iCNElemLoop - END IF ! VECNORM( ... - END IF ! BCType.EQ.BCTypeEPC - END IF ! BCIndex.GT.0 - END DO ! iSide = ElemInfo_Shared(ELEM_FIRSTSIDEIND,iGlobElem)+1,ElemInfo_Shared(ELEM_LASTSIDEIND,iGlobElem) - END DO iCNElemLoop ! iCNElem = 1,nComputeNodeTotalElems - END DO iElemLoop ! iElem = 1, nElems - END IF ! .NOT.(ALL(BConProc)) -#endif /*defined(PARTICLES)*/ +! Find my rank on the shared communicator, comm size and process name +IF(BConProc)THEN + CALL MPI_COMM_RANK(BiasVoltage%COMM%UNICATOR, BiasVoltage%COMM%MyRank, iError) + CALL MPI_COMM_SIZE(BiasVoltage%COMM%UNICATOR, BiasVoltage%COMM%nProcs, iError) -END IF ! MPIRoot + ! Inform about size of emission communicator + IF (BiasVoltage%COMM%MyRank.EQ.0) THEN +#if USE_LOADBALANCE + IF(.NOT.PerformLoadBalance)& +#endif /*USE_LOADBALANCE*/ + WRITE(UNIT_StdOut,'(A,I0,A,I0)') ' Bias voltage communicator on ',BiasVoltage%COMM%nProcs,' procs for BCState ',BCState + END IF +END IF ! BConProc +#endif /*USE_MPI*/ -! 5.) Create MPI sub-communicators -ALLOCATE(EPC%COMM(EPC%nUniqueEPCBounds)) -DO iUniqueEPCBC = 1, EPC%nUniqueEPCBounds - ! create new communicator - color = MERGE(iUniqueEPCBC, MPI_UNDEFINED, BConProc(iUniqueEPCBC)) +! When restarting, load the deposited charge on each EPC from the .h5 state file +CALL ReadBVDataFromH5() - ! set communicator id - EPC%COMM(iUniqueEPCBC)%ID=iUniqueEPCBC +END SUBROUTINE InitBV - ! create new emission communicator for electric potential boundary condition communication. Pass MPI_INFO_NULL as rank to follow the original ordering - CALL MPI_COMM_SPLIT(MPI_COMM_WORLD, color, MPI_INFO_NULL, EPC%COMM(iUniqueEPCBC)%UNICATOR, iError) - ! Find my rank on the shared communicator, comm size and proc name - IF(BConProc(iUniqueEPCBC))THEN - CALL MPI_COMM_RANK(EPC%COMM(iUniqueEPCBC)%UNICATOR, EPC%COMM(iUniqueEPCBC)%MyRank, iError) - CALL MPI_COMM_SIZE(EPC%COMM(iUniqueEPCBC)%UNICATOR, EPC%COMM(iUniqueEPCBC)%nProcs, iError) +!=================================================================================================================================== +!> Read the bias voltage (BV) data from a .h5 state file. +!> 1. The MPI root process reads the info and checks data consistency +!> 2. The MPI root process distributes the information among the sub-communicator processes connected to the BV boundary. +!=================================================================================================================================== +SUBROUTINE ReadBVDataFromH5() +! MODULES +USE MOD_io_hdf5 +USE MOD_Globals ,ONLY: UNIT_stdOut,MPIRoot,IK,abort +#if USE_LOADBALANCE +USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance,UseH5IOLoadBalance +#endif /*USE_LOADBALANCE*/ +USE MOD_IO_HDF5 ,ONLY: OpenDataFile,CloseDataFile,File_ID +USE MOD_Restart_Vars ,ONLY: DoRestart,RestartFile +USE MOD_HDF5_Input ,ONLY: DatasetExists,ReadArray,GetDataSize +USE MOD_HDG_Vars ,ONLY: BiasVoltage,BVDataLength +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! +! INPUT / OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +CHARACTER(255) :: ContainerName +LOGICAL :: BVExists +REAL :: BVDataHDF5(1:BVDataLength) +!=================================================================================================================================== +! Only required during restart +IF(.NOT.DoRestart) RETURN - ! inform about size of emission communicator - IF (EPC%COMM(iUniqueEPCBC)%MyRank.EQ.0) THEN #if USE_LOADBALANCE - IF(.NOT.PerformLoadBalance)& +! Do not try to read the data from .h5 if load balance is performed without creating a .h5 restart file +IF(PerformLoadBalance.AND..NOT.(UseH5IOLoadBalance)) RETURN #endif /*USE_LOADBALANCE*/ - WRITE(UNIT_StdOut,'(A,I0,A,I0,A,I0)') ' electric potential boundary condition: Emission-Communicator ',iUniqueEPCBC,' on ',& - EPC%COMM(iUniqueEPCBC)%nProcs,' procs for BCState ',EPC%BCState(iUniqueEPCBC) - END IF - END IF ! BConProc(iUniqueEPCBC) -END DO ! iUniqueEPCBC = 1, EPC%nUniqueEPCBounds -DEALLOCATE(BConProc) -! Get the number of procs that actually have a local BC side that is an EPC (required for voltage output to .csv) -! Procs might have zero EPC sides but are in the group because 1.) MPIRoot or 2.) the EPC is in the halo region -! Because only the MPI root process writes the .csv data, the information regarding the voltage on each EPC must be -! communicated with this process even though it might not be connected to each EPC boundary -DO iUniqueEPCBC = 1, EPC%nUniqueEPCBounds - ASSOCIATE( COMM => EPC%COMM(iUniqueEPCBC)%UNICATOR, nProcsWithSides => EPC%COMM(iUniqueEPCBC)%nProcsWithSides ) - IF(EPC%COMM(iUniqueEPCBC)%UNICATOR.NE.MPI_COMM_NULL)THEN - ! Check if the current processor is actually connected to the EPC via a BC side - IF(EPC%Group(EPC%BCState(iUniqueEPCBC),3).EQ.0)THEN - WithSides = 0 - ELSE - WithSides = 1 - END IF ! EPC%Group(EPC%BCState(iUniqueEPCBC),3).EQ.0 - ! Calculate the sum across the sub-communicator. Only the MPI root process needs this information - IF(MPIRoot)THEN - CALL MPI_REDUCE(WithSides, nProcsWithSides, 1 ,MPI_INTEGER, MPI_SUM, 0, COMM, iError) - ! Sanity check - IF(nProcsWithSides.EQ.0) CALL abort(__STAMP__,'Found EPC with no processors connected to it') - ELSE - CALL MPI_REDUCE(WithSides, 0 , 1 ,MPI_INTEGER, MPI_SUM, 0, COMM, IError) - END IF ! MPIRoot - END IF ! EPC%COMM(iUniqueEPCBC)%UNICATOR.NE.MPI_COMM_NULL - END ASSOCIATE -END DO ! iUniqueEPCBC = 1, EPC%nUniqueEPCBounds +! 1. The MPI root process reads the info and checks data consistency +! Only root reads the values and distributes them via MPI Broadcast +IF(MPIRoot)THEN + CALL OpenDataFile(RestartFile,create=.FALSE.,single=.TRUE.,readOnly=.TRUE.) + ! Check old parameter name + ContainerName='BiasVoltage' + CALL DatasetExists(File_ID,TRIM(ContainerName),BVExists) + ! Check for new parameter name + IF(BVExists)THEN + CALL ReadArray(TRIM(ContainerName) , 2 , (/1_IK , INT(BVDataLength,IK)/) , 0_IK , 1 , RealArray=BVDataHDF5) + WRITE(UNIT_stdOut,'(3(A,ES10.2E3))') " Read CoupledPowerPotential from restart file ["//TRIM(RestartFile)//& + "] Bias voltage[V]: ",BVDataHDF5(1),", Ion excess[C]: ",BVDataHDF5(2),", next adjustment time[s]: ",BVDataHDF5(3) + BiasVoltage%BVData = BVDataHDF5 + END IF ! BVExists + CALL CloseDataFile() +END IF ! MPIRoot + +#if USE_MPI +! 2. The MPI root process distributes the information among the sub-communicator processes for each EPC +CALL SynchronizeBV() #endif /*USE_MPI*/ -! When restarting, load the deposited charge on each EPC from the .h5 state file -CALL ReadEPCDataFromH5() +END SUBROUTINE ReadBVDataFromH5 -END SUBROUTINE InitBV + +!=================================================================================================================================== +!> Communicate the bias voltage values from MPIRoot to sub-communicator processes +!=================================================================================================================================== +SUBROUTINE SynchronizeBV() +! MODULES +USE MOD_Globals ,ONLY: IERROR,MPI_COMM_NULL,MPI_DOUBLE_PRECISION +USE MOD_HDG_Vars ,ONLY: BiasVoltage,BVDataLength +! insert modules here +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------! +! INPUT / OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +!=================================================================================================================================== +IF(BiasVoltage%COMM%UNICATOR.NE.MPI_COMM_NULL)THEN + ! Broadcast from root to other processors on the sub-communicator + CALL MPI_BCAST(BiasVoltage%BVData, BVDataLength, MPI_DOUBLE_PRECISION, 0, BiasVoltage%COMM%UNICATOR, IERROR) +END IF +END SUBROUTINE SynchronizeBV +#endif /*defined(PARTICLES)*/ !=================================================================================================================================== @@ -2012,6 +1951,11 @@ SUBROUTINE HDGLinear(time,U_out) r=q*(PP_N+1) + p+1 CALL ExactFunc(-4,Face_xGP(:,p,q,SideID),lambda(PP_nVar,r:r,SideID),t=time,BCState=BCState) END DO; END DO !p,q + CASE(50) ! exact BC = Dirichlet BC !! ExactFunc via DC bias voltage + DO q=0,PP_N; DO p=0,PP_N + r=q*(PP_N+1) + p+1 + CALL ExactFunc(-5,Face_xGP(:,p,q,SideID),lambda(PP_nVar,r:r,SideID),t=time,BCState=BCState) + END DO; END DO !p,q END SELECT ! BCType END DO !BCsideID=1,nDirichletBCSides #if (PP_nVar!=1) @@ -2452,8 +2396,8 @@ SUBROUTINE HDGNewton(time,U_out,td_iter,ForceCGSolverIteration_opt) r=q*(PP_N+1) + p+1 CALL ExactFunc(-3,Face_xGP(:,p,q,SideID),lambda(PP_nVar,r:r,SideID),t=time,iLinState=BCState) END DO; END DO !p,q - CASE(8) ! exact BC = Dirichlet BC !! ExactFunc via electric potential and decharing of a surface - CALL abort(__STAMP__,'Dirichlet BC=8 model not implemented for HDG Newton!') + CASE(8,50) ! exact BC = Dirichlet BC !! ExactFunc via electric potential and decharing of a surface + CALL abort(__STAMP__,'Dirichlet BC=8,50 model not implemented for HDG Newton!') END SELECT ! BCType END DO !BCsideID=1,nDirichletBCSides diff --git a/src/hdg/hdg_vars.f90 b/src/hdg/hdg_vars.f90 index 5a11d1793..9ccc3f015 100644 --- a/src/hdg/hdg_vars.f90 +++ b/src/hdg/hdg_vars.f90 @@ -238,15 +238,37 @@ MODULE MOD_HDG_Vars !-- Special BC with floating potential that is defined by the absorbed power of the charged particles !=================================================================================================================================== -LOGICAL :: CalcPCouplElectricPotential! Switch calculation on/off +LOGICAL :: CalcPCouplElectricPotential ! Switch calculation on/off #if USE_MPI -TYPE(tMPIGROUP) :: CPPCOMM !< communicator and ID for parallel execution +TYPE(tMPIGROUP) :: CPPCOMM !< communicator and ID for parallel execution #endif /*USE_MPI*/ -REAL :: CoupledPowerPotential(3) !< (/min, start, max/) electric potential at all BoundaryType = (/2,2/) -REAL :: CoupledPowerTarget !< Target input power at all BoundaryType = (/2,2/) -REAL :: CoupledPowerRelaxFac !< Relaxation factor for calculation of new electric potential +REAL :: CoupledPowerPotential(3) !< (/min, start, max/) electric potential at all BoundaryType = (/2,2/) +REAL :: CoupledPowerTarget !< Target input power at all BoundaryType = (/2,2/) +REAL :: CoupledPowerRelaxFac !< Relaxation factor for calculation of new electric potential + +!=================================================================================================================================== +!-- Bias Voltage (for calculating a BC voltage bias for certain BCs) +!=================================================================================================================================== + +LOGICAL :: UseBiasVoltage !< Automatic flag when bias voltage is to be used +INTEGER,PARAMETER :: BVDataLength=3 !< Number of variables in BVData + +TYPE tBV +#if USE_MPI + TYPE(tMPIGROUP) :: COMM !< communicator and ID for parallel execution +#endif /*USE_MPI*/ + INTEGER :: NPartBoundaries !< Total number of boundaries where the particles are counted + INTEGER,ALLOCATABLE :: PartBoundaries(:) !< Part-boundary number on which the particles are counted + REAL :: Frequency !< Frequency with which the bias voltage is adjusted (every period T = 1/f the bias voltage is changed) + REAL :: Delta !< Voltage difference used to change the current bias voltage (may also be adjusted over time automatically) + REAL :: BVData(BVDataLength) !< 1: bias voltage +! !< 2: Ion excess +! !< 3: sim. time when next adjustment happens +END TYPE + +TYPE(tBV) :: BiasVoltage #endif /*defined(PARTICLES)*/ !=================================================================================================================================== diff --git a/src/io_hdf5/hdf5_output_state.f90 b/src/io_hdf5/hdf5_output_state.f90 index ec948f4a3..7b5a12852 100644 --- a/src/io_hdf5/hdf5_output_state.f90 +++ b/src/io_hdf5/hdf5_output_state.f90 @@ -72,7 +72,7 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) USE MOD_Equation_Vars ,ONLY: E,Phi #endif /*PP_POIS*/ #if USE_HDG -USE MOD_HDG_Vars ,ONLY: nGP_face,iLocSides,UseFPC,FPC,UseEPC,EPC +USE MOD_HDG_Vars ,ONLY: nGP_face,iLocSides,UseFPC,FPC,UseEPC,EPC,UseBiasVoltage,BiasVoltage,BVDataLength #if PP_nVar==1 USE MOD_Equation_Vars ,ONLY: E,Et #elif PP_nVar==3 @@ -170,8 +170,8 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) #ifdef PARTICLES INTEGER :: i,j,k,iElem #endif /*PARTICLES*/ -REAL,ALLOCATABLE :: FPCDataHDF5(:,:),EPCDataHDF5(:,:) -INTEGER :: nVarFPC,nVarEPC +REAL,ALLOCATABLE :: FPCDataHDF5(:,:),EPCDataHDF5(:,:),BVDataHDF5(:,:) +INTEGER :: nVarFPC,nVarEPC,nVarBV #endif /*USE_HDG*/ !=================================================================================================================================== #ifdef EXTRAE @@ -591,6 +591,20 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) #endif /*PARTICLES*/ #if USE_HDG +! Bias voltage +IF(UseBiasVoltage.AND.MPIRoot)THEN + nVarBV = BVDataLength + ALLOCATE(BVDataHDF5(1:nVarBV,1)) + CALL OpenDataFile(FileName,create=.FALSE.,single=.TRUE.,readOnly=.FALSE.) + BVDataHDF5(1:nVarBV,1) = BiasVoltage%BVData(1:nVarBV) + CALL WriteArrayToHDF5( DataSetName = 'BiasVoltage' , rank = 2 , & + nValGlobal = (/1_IK , INT(nVarBV,IK)/), & + nVal = (/1_IK , INT(nVarBV,IK)/), & + offset = (/0_IK , 0_IK/) , & + collective = .FALSE., RealArray = BVDataHDF5(1:nVarBV,1)) + CALL CloseDataFile() + DEALLOCATE(BVDataHDF5) +END IF ! CalcBulkElectronTempi.AND.MPIRoot ! Floating boundary condition: Store charge on each FPC IF(UseFPC.AND.MPIRoot)THEN nVarFPC = FPC%nUniqueFPCBounds diff --git a/src/mesh/prepare_mesh.f90 b/src/mesh/prepare_mesh.f90 index 2715f4300..a02668ab6 100644 --- a/src/mesh/prepare_mesh.f90 +++ b/src/mesh/prepare_mesh.f90 @@ -1041,7 +1041,7 @@ SUBROUTINE fillMeshInfo() SELECT CASE(BCType) CASE(1) !periodic CALL abort(__STAMP__,'SideID.LE.nBCSides and SideID is periodic should not happen') - CASE(2,4,5,6,7,8,50) ! Dirichlet + CASE(HDGDIRICHLETBCSIDEIDS) ! Dirichlet ! do not consider this side CASE(10,11) ! Neumann HDGSides = HDGSides + 1 diff --git a/src/particles/surfacemodel/surfacemodel_analyze.f90 b/src/particles/surfacemodel/surfacemodel_analyze.f90 index f85dd0fc0..f8f4a8a40 100644 --- a/src/particles/surfacemodel/surfacemodel_analyze.f90 +++ b/src/particles/surfacemodel/surfacemodel_analyze.f90 @@ -158,6 +158,8 @@ SUBROUTINE AnalyzeSurface(Time) #if USE_HDG USE MOD_Analyze_Vars ,ONLY: EDC USE MOD_Analyze_Vars ,ONLY: CalcElectricTimeDerivative +USE MOD_HDG_Vars ,ONLY: UseBiasVoltage,BiasVoltage +USE MOD_HDG ,ONLY: SynchronizeBV #endif /*USE_HDG*/ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -175,7 +177,7 @@ SUBROUTINE AnalyzeSurface(Time) INTEGER :: iBC,iPartBound,iSEE,iBPO,iSpec REAL :: charge,TotalElectricCharge #if USE_HDG -INTEGER :: iEDCBC +INTEGER :: iEDCBC,i,iBoundary,iPartBound2 #endif /*USE_HDG*/ !=================================================================================================================================== IF((nComputeNodeSurfSides.EQ.0).AND.(.NOT.CalcBoundaryParticleOutput).AND.(.NOT.UseNeutralization).AND.(.NOT.CalcElectronSEE)) RETURN @@ -237,6 +239,20 @@ SUBROUTINE AnalyzeSurface(Time) END DO END IF ! BPO%OutputTotalElectricCurrent END IF +#if USE_HDG + ! Output of bias voltage variables + IF(UseBiasVoltage)THEN + WRITE(unit_index,'(A1)',ADVANCE='NO') ',' + WRITE(unit_index,'(I3.3,A)',ADVANCE='NO') OutputCounter,'-BiasVoltage[V]' + OutputCounter = OutputCounter + 1 + WRITE(unit_index,'(A1)',ADVANCE='NO') ',' + WRITE(unit_index,'(I3.3,A)',ADVANCE='NO') OutputCounter,'-BiasVoltageIonExcess[C]' + OutputCounter = OutputCounter + 1 + WRITE(unit_index,'(A1)',ADVANCE='NO') ',' + WRITE(unit_index,'(I3.3,A)',ADVANCE='NO') OutputCounter,'-BiasVoltageUpdateTime[s]' + OutputCounter = OutputCounter + 1 + END IF ! UseBiasVoltage +#endif /*USE_HDG*/ IF(UseNeutralization)THEN ! Ion thruster neutralization current (virtual cathode electrons) CALL WriteDataHeaderInfo(unit_index,'NeutralizationParticles',OutputCounter) END IF ! UseNeutralization @@ -298,7 +314,7 @@ SUBROUTINE AnalyzeSurface(Time) ! Output of total electric current IF(BPO%OutputTotalElectricCurrent)THEN DO iBPO = 1, BPO%NPartBoundaries - TotalElectricCharge = 0. ! initialize + TotalElectricCharge = 0. ! Initialize total charge for each boundary IF(ABS(SurfModelAnalyzeSampleTime).LE.0.0)THEN CALL WriteDataInfo(unit_index,1,RealArray=(/0.0/)) ELSE @@ -330,18 +346,79 @@ SUBROUTINE AnalyzeSurface(Time) ! Sampling time has already been considered due to the displacement current CALL WriteDataInfo(unit_index,RealScalar=TotalElectricCharge) END IF ! ABS(SurfModelAnalyzeSampleTime).LE.0.0 - END DO + END DO ! iBPO = 1, BPO%NPartBoundaries END IF ! BPO%OutputTotalElectricCurrent - ! Reset +#if USE_HDG + ! Bias voltage + IF(UseBiasVoltage)THEN + TotalElectricCharge = 0. ! Initialize sum over all boundaries + ! Ion excess + DO iBoundary = 1, BiasVoltage%NPartBoundaries + iPartBound = BiasVoltage%PartBoundaries(iBoundary) + iBPO = BPO%BCIDToBPOBCID(iPartBound) + ! Sum over all fluxes (only if the species has a charge) + DO iSpec = 1, BPO%NSpecies + charge = Species(BPO%Species(iSpec))%ChargeIC + ! Impacting charged particles: positive number for positive ions (+) and negative number for electrons (-) + IF(ABS(charge).GT.0.0) TotalElectricCharge = TotalElectricCharge + BPO%RealPartOut(iBPO,iSpec)*charge + END DO + ! Released secondary electrons (always a positive number). SEE%BCIDToSEEBCID(iPartBound) yields the iSEEBCIndex + IF(CalcElectronSEE)THEN + ! Get particle boundary ID + iPartBound2 = BPO%PartBoundaries(iBPO) + ! Sanity Check + IF(iPartBound.NE.iPartBound2) CALL abort(__STAMP__,'AnalyzeSurface(): Wrong particle boundary encountered!') + ! Get SEE ID + iSEE = SEE%BCIDToSEEBCID(iPartBound) + ! Add SEE current if this BC has secondary electron emission + IF(iSEE.GT.0) TotalElectricCharge = TotalElectricCharge + SEE%RealElectronOut(iSEE) + END IF ! CalcElectronSEE + END DO ! iBoundary = 1, BiasVoltage%NPartBoundaries + + ! Add total electric charge (without displacement current!) + BiasVoltage%BVData(2) = BiasVoltage%BVData(2) + TotalElectricCharge + + ! Simulation time threshold + IF(time.GE.BiasVoltage%BVData(3))THEN + ! Update time + IF(BiasVoltage%Frequency.GT.0.0)THEN + BiasVoltage%BVData(3) = BiasVoltage%BVData(3) + 1.0/BiasVoltage%Frequency + END IF ! BiasVoltage%Frequency.GT.0.0 + + ! Update Voltage + IF(BiasVoltage%BVData(2).GT.0.0)THEN + ! Increase voltage + BiasVoltage%BVData(1) = BiasVoltage%BVData(1) + BiasVoltage%Delta + ELSEIF(BiasVoltage%BVData(2).LT.0.0)THEN + ! Decrease voltage + BiasVoltage%BVData(1) = BiasVoltage%BVData(1) - BiasVoltage%Delta + ELSE + ! do nothing + END IF ! BiasVoltage%BVData(2).LT.0.0 + + ! Reset ion excess counter + BiasVoltage%BVData(2) = 0. + END IF ! time.GE.BiasVoltage%BVData(3) + + ! Write: Voltage, Ion excess and simulation update time + DO i = 1, 3 + CALL WriteDataInfo(unit_index, RealScalar=BiasVoltage%BVData(i)) + END DO ! i = 1, 3 + END IF ! UseBiasVoltage +#endif /*USE_HDG*/ + + ! Reset BPO containers DO iPartBound = 1, BPO%NPartBoundaries DO iSpec = 1, BPO%NSpecies ! Reset PartMPI%MPIRoot counters after writing the data to the file, ! non-PartMPI%MPIRoot are reset in SyncBoundaryParticleOutput() BPO%RealPartOut(iPartBound,iSpec) = 0. - END DO - END DO - END IF + END DO ! iSpec = 1, BPO%NSpecies + END DO ! iPartBound = 1, BPO%NPartBoundaries + END IF ! CalcBoundaryParticleOutput + IF(UseNeutralization) CALL WriteDataInfo(unit_index,RealScalar=REAL(NeutralizationBalanceGlobal)) + IF(CalcElectronSEE)THEN DO iPartBound = 1, SEE%NPartBoundaries IF(ABS(SurfModelAnalyzeSampleTime).LE.0.0)THEN @@ -357,6 +434,11 @@ SUBROUTINE AnalyzeSurface(Time) WRITE(unit_index,'(A)') '' #if USE_MPI END IF +!----------------------------------------------------------------------------------------------------------------------------------- +#if USE_HDG + ! Bias voltage: Update values of sub-communicator processes (broadcast from MPIRoot to all sub-communicator processes) + IF(UseBiasVoltage) CALL SynchronizeBV() +#endif /*USE_HDG*/ #endif /*USE_MPI*/ !----------------------------------------------------------------------------------------------------------------------------------- SurfModelAnalyzeSampleTime = Time ! Backup "old" time value for next output diff --git a/src/piclas.h b/src/piclas.h index e194f97e4..e10bc2548 100644 --- a/src/piclas.h +++ b/src/piclas.h @@ -321,3 +321,8 @@ ! Secondary electron emission #define SEE_MODELS_ID 5,6,7,8,9,10,11 + +#if USE_HDG +! HDG Dirichlet BC Side IDs +#define HDGDIRICHLETBCSIDEIDS 2,4,5,6,7,8,50 +#endif diff --git a/src/restart/restart_field.f90 b/src/restart/restart_field.f90 index 9611c9fdc..ccb205d67 100644 --- a/src/restart/restart_field.f90 +++ b/src/restart/restart_field.f90 @@ -74,8 +74,9 @@ SUBROUTINE FieldRestart() USE MOD_MPI ,ONLY: StartReceiveMPIData,StartSendMPIData,FinishExchangeMPIData #endif /*USE_MPI*/ #if USE_LOADBALANCE -USE MOD_HDG ,ONLY: SynchronizeVoltageOnEPC -USE MOD_HDG_Vars ,ONLY: UseEPC +USE MOD_Equation ,ONLY: SynchronizeCPP +USE MOD_HDG ,ONLY: SynchronizeVoltageOnEPC,SynchronizeBV +USE MOD_HDG_Vars ,ONLY: UseEPC,UseBiasVoltage,CalcPCouplElectricPotential #if defined(PARTICLES) ! TODO: make ElemInfo available with PARTICLES=OFF and remove this preprocessor if/else as soon as possible USE MOD_Mesh_Vars ,ONLY: SideToNonUniqueGlobalSide,nElems @@ -152,6 +153,13 @@ SUBROUTINE FieldRestart() #if USE_LOADBALANCE IF(PerformLoadBalance.AND.(.NOT.UseH5IOLoadBalance))THEN #if USE_HDG + ! Coupled Bias voltage (BV): The MPI root process distributes the information among the sub-communicator processes + ! (before and after load balancing, the root process is always part of each sub-communicator group) + IF(UseBiasVoltage) CALL SynchronizeBV() + + ! Coupled Power Potential (CPP): The MPI root process distributes the information among the sub-communicator processes + ! (before and after load balancing, the root process is always part of each sub-communicator group) + IF(CalcPCouplElectricPotential) CALL SynchronizeCPP() #if USE_PETSC ! FPC: The MPI root process distributes the information among the sub-communicator processes for each FPC From 4f845a07b289ed7344299d7862a5e9342928d53f Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 19 Apr 2023 15:53:32 +0200 Subject: [PATCH 04/11] Added reggie for bias voltage calculation based on the 1D "turner" setup but with DC potential instead of the usual AC boundary condition. --- REGGIE.md | 17 +- .../turner_bias-voltage_DC/DSMC.ini | 50 +++++ .../DSMCSpecies_electronic_state_full_Data.h5 | Bin 0 -> 12448 bytes .../PartAnalyze_ref.csv | 24 +++ .../turner_bias-voltage_DC/analyze.ini | 14 ++ .../turner_bias-voltage_DC/command_line.ini | 2 + .../turner_bias-voltage_DC/excludeBuild.ini | 2 + .../turner_bias-voltage_DC/externals.ini | 6 + .../turner_bias-voltage_DC/parameter.ini | 191 ++++++++++++++++++ .../turner_bias-voltage_DC/pre-hopr/hopr.ini | 45 +++++ .../turner_bias-voltage_DC/readme.md | 5 + src/hdg/hdg.f90 | 34 ++-- src/hdg/hdg_vars.f90 | 2 +- src/io_hdf5/hdf5_output_state.f90 | 13 +- .../surfacemodel/surfacemodel_analyze.f90 | 48 ++--- 15 files changed, 395 insertions(+), 58 deletions(-) create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMC.ini create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMCSpecies_electronic_state_full_Data.h5 create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/PartAnalyze_ref.csv create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/command_line.ini create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/excludeBuild.ini create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr/hopr.ini create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md diff --git a/REGGIE.md b/REGGIE.md index fa60d8074..7d7f72055 100644 --- a/REGGIE.md +++ b/REGGIE.md @@ -328,16 +328,17 @@ Testing PIC compiled with Boris-Leapfrog integration (poisson,Boris-Leapfrog), s Testing PIC compiled with Runge-Kutta 3 integration, solving Poisson's equation: [Link to build](regressioncheck/NIG_PIC_poisson_RK3/builds.ini). -| **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | +| **No.** | **Case** | **CMAKE-CONFIG** | **Feature** | **Execution** | **Comparing** | **Readme** | | :-----: | :-----------------------------------------------------------: | :--------------: | :-----------------------------------------------------------------------------------------------------------------------------: | :-------------: | :--------------------------------------: | :--------------------------------------------------------------------------------------------------------------------------------------: | -| 1 | parallel_plates | | CalcCoupledPower | nProcs=1 | PartAnalyzeRK3_ref.csv | [Link](regressioncheck/NIG_PIC_poisson_RK3/parallel_plates/readme.md) | -| 2 | parallel_plates_AC | | CalcCoupledPower | nProcs=1 | PartAnalyzeRK3_ref.csv | [Link](regressioncheck/NIG_PIC_poisson_RK3/parallel_plates_AC/readme.md) | -| 3 | plasma_sheath_BR-electrons_conforming | | non-linear HDG (BR electrons) | nProcs=2 | TimeAvg | [Link](regressioncheck/NIG_PIC_poisson_RK3/plasma_sheath_BR-electrons_conforming/readme.md) | -| 4 | plasma_sheath_BR-electrons_conforming_auto-switch | | non-linear HDG (BR electrons), automatic switching BR/kinetic | nProcs=1,2,4 | - | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_conforming_auto-switch/readme.md) | -| 5 | plasma_sheath_BR-electrons_conforming_auto-switch_auto-ref | | non-linear HDG (BR electrons), automatic switching BR/kinetic, automatic ref. values, change nSkipAnalyze during the simulation | nProcs=1,2,4,11 | integrate Te over time (PartAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_conforming_auto-switch_auto-ref/readme.md) | +| 1 | parallel_plates | | CalcCoupledPower | nProcs=1 | PartAnalyzeRK3_ref.csv | [Link](regressioncheck/NIG_PIC_poisson_RK3/parallel_plates/readme.md) | +| 2 | parallel_plates_AC | | CalcCoupledPower | nProcs=1 | PartAnalyzeRK3_ref.csv | [Link](regressioncheck/NIG_PIC_poisson_RK3/parallel_plates_AC/readme.md) | +| 3 | plasma_sheath_BR-electrons_conforming | | non-linear HDG (BR electrons) | nProcs=2 | TimeAvg | [Link](regressioncheck/NIG_PIC_poisson_RK3/plasma_sheath_BR-electrons_conforming/readme.md) | +| 4 | plasma_sheath_BR-electrons_conforming_auto-switch | | non-linear HDG (BR electrons), automatic switching BR/kinetic | nProcs=1,2,4 | - | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_conforming_auto-switch/readme.md) | +| 5 | plasma_sheath_BR-electrons_conforming_auto-switch_auto-ref | | non-linear HDG (BR electrons), automatic switching BR/kinetic, automatic ref. values, change nSkipAnalyze during the simulation | nProcs=1,2,4,11 | integrate Te over time (PartAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_conforming_auto-switch_auto-ref/readme.md) | | 6 | plasma_sheath_BR-electrons_conforming_auto-switch_variable_Te | | non-linear HDG (BR electrons), automatic switching BR/kinetic, variable Te, change nSkipAnalyze during the simulation | nProcs=1,2,4,11 | integrate Te over time (PartAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_conforming_auto-switch_variable_Te/readme.md) | -| 7 | plasma_sheath_BR-electrons_mortar | | non-linear HDG (BR electrons), Mortars | nProcs=2 | TimeAvg | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_mortar/readme.md) | -| 8 | turner | | | nProcs=4 | L2 error, PartAnalyze.csv | | +| 7 | plasma_sheath_BR-electrons_mortar | | non-linear HDG (BR electrons), Mortars | nProcs=2 | TimeAvg | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_mortar/readme.md) | +| 8 | turner | | | nProcs=4 | L2 error, PartAnalyze.csv | | +| 9 | turner_bias-voltage_DC | | bias voltage for DC potential boundaries (BCType=50) | nProcs=1,2,4,10 | PartAnalyze.csv, SurfaceAnalyze.csv | [Link](regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md) | ### NIG_PIC_maxwell_RK4 diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMC.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMC.ini new file mode 100644 index 000000000..5794e1ff9 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMC.ini @@ -0,0 +1,50 @@ +! =============================================================================== ! +! Species1, He +! =============================================================================== ! +Part-Species1-SpeciesName = He +Part-Species1-InteractionID = 1 +Part-Species1-Tref = 300 +Part-Species1-dref = 3.0E-10 +Part-Species1-omega=0.2 +Part-Species1-HeatOfFormation_K=0.0 +! =============================================================================== ! +! Species2, e +! =============================================================================== ! +Part-Species2-SpeciesName = electron +Part-Species2-InteractionID = 4 +Part-Species2-Tref = 300 +Part-Species2-dref = 2.817920E-15 +Part-Species2-omega=0.2 +! =============================================================================== ! +! Species3, He+ +! =============================================================================== ! +Part-Species3-SpeciesName = HeIon1 +Part-Species3-InteractionID = 10 ! atomar ion +Part-Species3-Tref = 300 +Part-Species3-dref = 3.0E-10 +Part-Species3-omega=0.2 +Part-Species3-PreviousState=1 +! =============================================================================== ! +! Species4, He+2 +! =============================================================================== ! +Part-Species4-SpeciesName = HeIon2 +Part-Species4-InteractionID = 4 ! es verhalte sich wie 4 +Part-Species4-Tref = 300 +Part-Species4-dref = 3.0E-10 +Part-Species4-omega=0.2 +Part-Species4-PreviousState=3 +! =============================================================================== ! +! Data for chemical reactions +! =============================================================================== ! +DSMC-NumOfReactions=2 +! =============================================================================== ! +! Ionization +! =============================================================================== ! +! Reaction 1 | He + e --> HeIon1 + e + e +DSMC-Reaction1-ReactionModel=QK +DSMC-Reaction1-Reactants=(/1,2,0/) +DSMC-Reaction1-Products=(/3,2,2,0/) +! Reaction 2 | HeIon1 + e --> HeIon2 + e + e +DSMC-Reaction2-ReactionModel=QK +DSMC-Reaction2-Reactants=(/3,2,0/) +DSMC-Reaction2-Products=(/4,2,2,0/) diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMCSpecies_electronic_state_full_Data.h5 b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMCSpecies_electronic_state_full_Data.h5 new file mode 100644 index 0000000000000000000000000000000000000000..d4b92d1500d35cf2ce0b16d5b5fbebc11020ab5c GIT binary patch literal 12448 zcmeI&3s4o;835qrAsh)1z!eFs^s>sLn|S46Ktm_{^9mpdPm79!-o$EEMlgnIFbx$G zEoy_YNHtc4f+E$RMxz-)AfyV@gI$pZdS8Mfo}bG_77-hZ8JM2!j<75G1@ zfT}!rvILg$aWQj;xJKEIt4j06B@F@N+%)rWsqQ!XUfxShAnVq_bqjT?O6Q*X?$>n_ zh+JG9Q0(rmjVJI}+p9-2hI#1oUrNoI)D=n1gIl7*^qM4U0;jTP7NS)}-WA}sLyrqE z!ReX9b(3j64_C{(+Ry~9opE)kXW}1nd-FfnP2kG0g7}~;I-bDJbLIlC1JtTgb)w{Q zQ42K>Gb5U~`Xkf?(lR)Aw5qiBuik?eAYtOUYp1`U`8=fLRdu5Y%x}Kp-p<7LZk>^G zLALB|$3IZ_jtlsq;Sttv$+=l_TnU|Z4k++1gx&%>JZ&O9`M?2k{pj^^R) z*@C*FP8E3u3~3-HDU z?;x}YXTB*LjV9m&Uq8VSdj13!*(L1H^ru#p_6gqsWcp#@kRKPDS3Ldf~@>1YYE!<(DXGNia9HKVQT zRkY_Lv;@ftSsN1vm89w_Xpqfq_Y3v>vLGwdrtQ>(kUcB-akFITub%;PclC zcc4{8T(qwhEkNm0V{(s2(0m?-dlfUbuIDe_d4-M7{?OufO3xQytI2oh*SbYG-&p94 z=HYZrhYwnSeVH!^Xc2yu_gQdNX2!&9E& zzeJ0$@5y0{&=PnBHZhjLd-VL37&pT)!EY5>g0L+I8Ovb*+lbW|CoukRj+@DNT0a7p zs}F9i0BTiv`;+vmanaP)?;DTCYkj&2NWaVIMytxI{x+w6JuX0PXPIwh433WrKf#zl zbw`9h8=smv_)?XgZ@oSn6K`zSErLU=xd+Wdt*PZET7V(jLT{l(=v_8qe~ms~f|8-) z0kiKxpRo+Jl^NY^{?Gl7euoZ$&WQoNXc_v(FJNqjCF|Z}OyFk3l64>J@1Hw9`3H_o|%wDP2k9o$7i8cWyH*#3(*2x98`LtPS5ATGpUg=0mmTI zB{m*nxh)^k^967&4RC7GEy69ifiVw0AzwLT+`3-7%Lk)H5TCnJf|lTR|CTjq0ZKZe z-bRa%e{VNq35vh$-i>h?BA3nCgBBqo=Qv{t{PF|d!?+CBMy{+wn<3}kSBxbHUK?1C zaT)Tzv@kY<$yD8l@et_S_E*L-JU`U83FBs%dpDjjfxD?k*PqnuN8sDxjlWczkIY;&51l6i^Uwms1Qskni?F{o<2c&7URO^swq6h0pFe?d z5k}qF#8?9V=v{xuxC}35IQ#=`{k#{aGM3;aX%%A`GCo~#9rLZVcYdpH2()&I1Qy%ee@{j%B-(U5>$rvyyU{t`UfKdUX0zc>qJlqdd Om1QrqE_^XMEdHPWK*b0E literal 0 HcmV?d00001 diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/PartAnalyze_ref.csv b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/PartAnalyze_ref.csv new file mode 100644 index 000000000..fbf3ffa2d --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/PartAnalyze_ref.csv @@ -0,0 +1,24 @@ +001-TIME,002-nPart-Spec-001,003-nPart-Spec-002,004-nPart-Spec-003,005-nPart-Spec-004,006-nPart-Spec-005,007-E-Vib001,008-E-Vib002,009-E-Vib003,010-E-Vib004,011-E-Vib005,012-E-Rot001,013-E-Rot002,014-E-Rot003,015-E-Rot004,016-E-Rot005,017-E-Elec001,018-E-Elec002,019-E-Elec003,020-E-Elec004,021-E-Elec005,022-E-TotalPart +0.0000000000000000E+000,0.0000000000000000E+000,0.4012000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8024000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2517594416481664E-013 +0.4580000000000000E-010,0.0000000000000000E+000,0.4011000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8023000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2518698200783414E-013 +0.9160000000000000E-010,0.0000000000000000E+000,0.4010000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8022000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2522462715369796E-013 +0.1374000000000000E-009,0.0000000000000000E+000,0.4008000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8020000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2528153121183654E-013 +0.1832000000000000E-009,0.0000000000000000E+000,0.4008000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8020000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2537316774570108E-013 +0.2290000000000000E-009,0.0000000000000000E+000,0.4007000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8019000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2547895782038987E-013 +0.2748000000000000E-009,0.0000000000000000E+000,0.4007000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8019000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2561151429942829E-013 +0.3206000000000000E-009,0.0000000000000000E+000,0.4005000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8017000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2575341266174985E-013 +0.3663999999999999E-009,0.0000000000000000E+000,0.4004000000000000E+004,0.4012000000000000E+004,0.0000000000000000E+000,0.8016000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2592967043329345E-013 +0.4121999999999999E-009,0.0000000000000000E+000,0.4002000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8015000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2606812409604058E-013 +0.4579999999999999E-009,0.0000000000000000E+000,0.4000000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8013000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2624571920908253E-013 +0.5000000000000000E-009,0.0000000000000000E+000,0.3999000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8012000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2641262427096647E-013 +0.5458000000000000E-009,0.0000000000000000E+000,0.3997000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8010000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2657709585005035E-013 +0.5916000000000000E-009,0.0000000000000000E+000,0.3997000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8010000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2680186010280453E-013 +0.6374000000000000E-009,0.0000000000000000E+000,0.3996000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8009000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2703333139144175E-013 +0.6832000000000000E-009,0.0000000000000000E+000,0.3993000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8006000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2722071667924760E-013 +0.7290000000000000E-009,0.0000000000000000E+000,0.3991000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8004000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2736908745529201E-013 +0.7748000000000000E-009,0.0000000000000000E+000,0.3989000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.8002000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2756365047732255E-013 +0.8205999999999999E-009,0.0000000000000000E+000,0.3985000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.7998000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2774018444632331E-013 +0.8663999999999999E-009,0.0000000000000000E+000,0.3982000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.7995000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2784862968426264E-013 +0.9121999999999999E-009,0.0000000000000000E+000,0.3980000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.7993000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2799497523940786E-013 +0.9580000000000000E-009,0.0000000000000000E+000,0.3978000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.7991000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2813652599466704E-013 +0.1000000000000000E-008,0.0000000000000000E+000,0.3978000000000000E+004,0.4013000000000000E+004,0.0000000000000000E+000,0.7991000000000000E+004,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.0000000000000000E+000,0.2828197009080410E-013 diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini new file mode 100644 index 000000000..757b7efac --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini @@ -0,0 +1,14 @@ +! compare the last row in PartAnalyze.csv with a reference file +compare_data_file_name = PartAnalyze.csv +compare_data_file_reference = PartAnalyze_ref.csv +compare_data_file_tolerance = 20.0e-2 +compare_data_file_tolerance_type = relative +compare_data_file_max_differences = 2 + +! integrate columns x:y in a data file as integral(y(x), x, x(1), x(end)) +integrate_line_file = SurfaceAnalyze.csv ! data file name +integrate_line_columns = 0:9 ! columns x:y +integrate_line_integral_value = -1.0172 ! bias voltage over time +integrate_line_tolerance_value = 1.0e-3 ! tolerance +integrate_line_tolerance_type = relative ! absolute or relative comparison +integrate_line_multiplier = 1e9 ! 1/tEnd diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/command_line.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/command_line.ini new file mode 100644 index 000000000..bbe3efc24 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/command_line.ini @@ -0,0 +1,2 @@ +MPI = 1,2,4,10 +cmd_suffix = DSMC.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/excludeBuild.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/excludeBuild.ini new file mode 100644 index 000000000..f8a3a0b81 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/excludeBuild.ini @@ -0,0 +1,2 @@ +! maxwell RK +PICLAS_EQNSYSNAME=maxwell diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini new file mode 100644 index 000000000..aa8f4729f --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini @@ -0,0 +1,6 @@ +! --- Externals Tool Reggie +MPI = 1 +externalbinary = ./hopr/build/bin/hopr +externaldirectory = pre-hopr +externalruntime = pre +cmd_pre_execute = ln\s-s\s./pre-hopr/turner2013_mesh.h5\s../. diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini new file mode 100644 index 000000000..c2190c375 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini @@ -0,0 +1,191 @@ +! =============================================================================== ! +! EQUATION (linearscalaradvection) +! =============================================================================== ! +IniExactFunc = 0 ! empty +! =============================================================================== ! +! DISCRETIZATION +! =============================================================================== ! +N = 1 ! Polynomial degree +NAnalyze = 10 ! Number of analyze points +! =============================================================================== ! +! MESH +! =============================================================================== ! +MeshFile = ./turner2013_mesh.h5 +useCurveds = F +! =============================================================================== ! +! OUTPUT / VISUALIZATION +! =============================================================================== ! +ProjectName = turner2013 +TrackingMethod = triatracking +Part-SafetyFactor = 1 +Logging = F +WriteErrorFiles = F +printRandomSeeds = F +DoCalcErrorNorms = T +FlushInitialState = T +! =============================================================================== ! +! ANALYZE +! =============================================================================== ! +!CalcKineticEnergy = T +CalcPotentialEnergy = T +CalcNumSpec = T +CalcInternalEnergy = T +!CalcTemp = T +CalcPartBalance = F +CalcVelos = F +VelocityDirections = (/1,0,0,1/) ! x,y,z,abs + +PIC-OutputSource = T + +CalcCollRates = T ! piclas +CalcReacRates = T ! piclas +Particles-DSMC-CalcQualityFactors = F ! piclas: Pmax/Pmean + +CalcPointsPerDebyeLength = T +CalcPICCFLCondition = T +CalcMaxPartDisplacement = T +! =============================================================================== ! +! CALCULATION +! =============================================================================== ! +DoLoadBalance = T +DoInitialAutoRestart = T +Load-DeviationThreshold = 1e-3 +LoadBalanceMaxSteps = 10 + +ManualTimestep = 4.58E-11 +IterDisplayStep = 10 +tend = 1.0E-9 +Analyze_dt = 5.0E-10 +! =============================================================================== ! +! Particle Boundaries +! =============================================================================== ! +BoundaryName = BC_left +BoundaryType = (/50,0/) ! Dirichlet with 0V initial BC + +BoundaryName = BC_right +BoundaryType = (/4,0/) ! 4: Dirichlet with zero potential + +! Automatic adjustment of the DC potential +UseBiasVoltage = T +BiasVoltage-NPartBoundaries = 2 ! Nbr of particle boundaries where the ion excess is measured +Biasvoltage-PartBoundaries = (/1,2/) ! Particle boundary index of where the ion excess is measured +BiasVoltage-Frequency = 3e9 ! Frequency with which the bias voltage is adapted +BiasVoltage-Delta = 1.0 ! Voltage difference used for adapting the bias voltage +! =============================================================================== ! +! PARTICLES +! =============================================================================== ! +Part-maxParticleNumber = 500000 +Part-nSpecies = 4 +PIC-externalField = (/0.,0.,0.,0.,0.,0./) + +Part-FIBGMdeltas = (/0.00013,3.42e-5,3.42e-5/) +PIC-Deposition-Type = cell_volweight + +Part-nBounds = 6 + +Part-Boundary1-SourceName = BC_left +Part-Boundary1-Condition = open + +Part-Boundary2-SourceName = BC_right +Part-Boundary2-Condition = open + +Part-Boundary3-SourceName = BC_periodicy+ +Part-Boundary3-Condition = periodic + +Part-Boundary4-SourceName = BC_periodicy- +Part-Boundary4-Condition = periodic + +Part-Boundary5-SourceName = BC_periodicz+ +Part-Boundary5-Condition = periodic + +Part-Boundary6-SourceName = BC_periodicz- +Part-Boundary6-Condition = periodic + +Part-nPeriodicVectors = 2 + +! Output of integral particle flux on specific BCs +CalcBoundaryParticleOutput = T +BPO-NPartBoundaries = 2 +BPO-PartBoundaries = (/1,2/) +BPO-NSpecies = 3 +BPO-Species = (/2,3,4/) + +! =============================================================================== ! +! DSMC +! =============================================================================== ! +UseDSMC = T +Particles-DSMC-CollisMode = 3 !(1:elast coll, 2: elast + rela, 3:chem) +Particles-DSMC-ElectronicModel = 1 +Particles-DSMCElectronicDatabase = DSMCSpecies_electronic_state_full_Data.h5 ! when supplied: doQK = true +EpsMergeElectronicState = 1.e-2 ! merge QK levels when difference falls below eps + +Part-NumberOfRandomSeeds = 2 +Particles-RandomSeed1 = 1 +Particles-RandomSeed2 = 2 + +Particles-HaloEpsVelo = 3E5 ! 300E6 + +! HDG +epsCG = 1e-6 +maxIterCG = 10000 !'500' + +! =============================================================================== ! +! Species1 | He +! =============================================================================== ! +Part-Species1-ChargeIC = 0 +Part-Species1-MassIC = 6.64647640919434E-027 +Part-Species1-MacroParticleFactor = 10 +Part-Species1-nInits = 1 + +Part-Species1-Init1-SpaceIC = background +Part-Species1-Init1-PartDensity = 96.4E+20 +Part-Species1-Init1-velocityDistribution = maxwell_lpn +Part-Species1-Init1-MWTemperatureIC = 300.0 +Part-Species1-Init1-VeloIC = 0 +Part-Species1-Init1-VeloVecIC = (/0.,0.,1./) +Part-Species1-Init1-Tempelec = 300.0 +! =============================================================================== ! +! Species2 | e +! =============================================================================== ! +Part-Species2-ChargeIC = -1.60217653E-19 +Part-Species2-MassIC = 9.1093826E-31 +Part-Species2-MacroParticleFactor = 10 +Part-Species2-nInits = 1 + +Part-Species2-Init1-SpaceIC = cuboid +Part-Species2-Init1-velocityDistribution = maxwell +Part-Species2-Init1-MWTemperatureIC = 30000.0 +Part-Species2-Init1-PartDensity = 5.12E14 +Part-Species2-Init1-BasePointIC = (/0.,0.,0./) +Part-Species2-Init1-BaseVector1IC = (/0.,3.42e-5,0./) +Part-Species2-Init1-BaseVector2IC = (/0.,0.,3.42e-5/) +Part-Species2-Init1-NormalIC = (/1.,0.,0./) +Part-Species2-Init1-CuboidHeightIC = 0.067 +Part-Species2-Init1-VeloIC = 0 +Part-Species2-Init1-VeloVecIC = (/0.,0.,1./) +! =============================================================================== ! +! Species3 | HeIon +! =============================================================================== ! +Part-Species3-ChargeIC = 1.60217653E-19 +Part-Species3-MassIC = 6.645565470903E-027 +Part-Species3-MacroParticleFactor = 10 +Part-Species3-nInits = 1 + +Part-Species3-Init1-SpaceIC = cuboid +Part-Species3-Init1-velocityDistribution = maxwell +Part-Species3-Init1-MWTemperatureIC = 300.0 +Part-Species3-Init1-PartDensity = 5.12E14 +Part-Species3-Init1-BasePointIC = (/0.,0.,0./) +Part-Species3-Init1-BaseVector1IC = (/0.,3.42e-5,0./) +Part-Species3-Init1-BaseVector2IC = (/0.,0.,3.42e-5/) +Part-Species3-Init1-NormalIC = (/1.,0.,0./) +Part-Species3-Init1-CuboidHeightIC = 0.067 +Part-Species3-Init1-VeloIC = 0 +Part-Species3-Init1-VeloVecIC = (/0.,0.,1./) +Part-Species3-Init1-Tempelec = 300.0 +! =============================================================================== ! +! Species3 | HeIon2 +! =============================================================================== ! +Part-Species4-ChargeIC = 3.20435306E-019 +Part-Species4-MassIC = 6.64465453261166E-027 +Part-Species4-MacroParticleFactor = 10 diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr/hopr.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr/hopr.ini new file mode 100644 index 000000000..8c120e058 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr/hopr.ini @@ -0,0 +1,45 @@ +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! +ProjectName = turner2013 ! name of the project (used for filenames) +Debugvisu =T ! Write debug mesh to tecplot file +Logging =F ! Write log files + +!=============================================================================== ! +! MESH +!=============================================================================== ! +Mode =1 ! 1 Cartesian 2 gambit file 3 CGNS +nZones =1 ! number of zones +Corner =(/0.,0.,0.0,,6.7,0.,0.0,,6.7,3.42e-3,0.0,,0.,3.42e-3,0.0 ,,0.,0.,3.42e-3,,6.7,0.,3.42e-3,,6.7,3.42e-3,3.42e-3,,0.,3.42e-3,3.42e-3/) ! [-3,3]x[-3,3]x[-3,3] +nElems =(/512,1,1/) ! Anzahl der Elemente in jede Richtung +BCIndex =(/1,3,6,4,5,2/) ! Indices of Boundary Conditions for six Boundary Faces (z-,y-,x+,y+,x-,z+) +elemtype =108 ! Elementform (108: Hexaeder) +useCurveds =F ! T if curved boundaries defined +SpaceQuandt =1. ! characteristic length of the mesh +ConformConnect=T + +postScaleMesh = T +MeshScale = 1e-2 +!=============================================================================== ! +! BOUNDARY CONDITIONS +!=============================================================================== ! + nUserDefinedBoundaries=6 + + +BoundaryName=BC_periodicz- ! BC index 1 (from position in parameterfile) +BoundaryType=(/1,0,0,1/) ! (/ Type, curveIndex, State, alpha /) +BoundaryName=BC_periodicz+ ! BC index 2 +BoundaryType=(/1,0,0,-1/) ! here the direction of the vector 1 is changed, because it is the opposite side +vv=(/0.,0.,3.42e-3/) ! vector for periodic BC in z direction (zminus,zplus), index=1 + +BoundaryName=BC_periodicy- ! BC index 3 +BoundaryType=(/1,0,0,2/) +BoundaryName=BC_periodicy+ ! BC index 4 +BoundaryType=(/1,0,0,-2/) ! (/ BCType=1: periodic, 0, 0, Index of second vector vv in parameter file /) +vv=(/0.,3.42e-3,0./) ! vector for periodic BC in y direction (yminus,yplus), index=2 + + +BoundaryName=BC_left +BoundaryType=(/4,0,0,0/) ! ideal conductor +BoundaryName=BC_right +BoundaryType=(/4,0,0,0/) ! ideal conductor diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md new file mode 100644 index 000000000..1fb6ef241 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md @@ -0,0 +1,5 @@ +# 1D turner setup and bias voltage feature +- He plasma between two electrodes one of which is usually AC, but for this bias-voltage test case it is set to DC starting at 0V +- bias voltage for DC only (field BCType=50) +- ion excess on both electrodes generates a positive bias voltage and electron excess a negative potential +- comparison of PartAnalyze.csv with reference file and integration of the bias voltage over time from SurfaceAnalyze.csv diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index ed1282c95..cd4961eac 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -1255,16 +1255,17 @@ SUBROUTINE InitBV() ! Activate model UseBiasVoltage = GETLOGICAL('UseBiasVoltage') +! Count the number of boundaries that allow bias voltage +BVBoundaries = 0 +DO iBC=1,nBCs + BCType = BoundaryType(iBC,BC_TYPE) + IF(.NOT.ANY(BCType.EQ.BCTypeBV)) CYCLE ! Skip other boundaries + BVBoundaries = BVBoundaries + 1 +END DO + ! Skip the following if bias voltage is not active IF(.NOT.UseBiasVoltage)THEN ! Sanity check before returning: bias voltage BCs cannot be used without activating the bias voltage model - ! Count the number of boundaries that allow bias voltage - BVBoundaries = 0 - DO iBC=1,nBCs - BCType = BoundaryType(iBC,BC_TYPE) - IF(.NOT.ANY(BCType.EQ.BCTypeBV)) CYCLE ! Skip other boundaries - BVBoundaries = BVBoundaries + 1 - END DO IF(BVBoundaries.GT.0) CALL CollectiveStop(__STAMP__,' Bias voltage BCs require UseBiasVoltage=T!') ! Exit this subroutine @@ -1274,14 +1275,7 @@ SUBROUTINE InitBV() ! CalcBoundaryParticleOutput=T and boundaries must be set correctly IF(.NOT.CalcBoundaryParticleOutput) CALL CollectiveStop(__STAMP__,' UseBiasVoltage=T requires CalcBoundaryParticleOutput=T!') -! Count the number of boundaries that allow bias voltage -BVBoundaries = 0 -DO iBC=1,nBCs - BCType = BoundaryType(iBC,BC_TYPE) - IF(.NOT.ANY(BCType.EQ.BCTypeBV)) CYCLE ! Skip other boundaries - BVBoundaries = BVBoundaries + 1 -END DO - +! Check the number of boundaries that allow bias voltage: Must be exactly 1 IF(BVBoundaries.NE.1) CALL CollectiveStop(__STAMP__,' UseBiasVoltage=T requires exactly one boundary with this feature!') !> 2.) Get bias voltage parameters @@ -1290,10 +1284,15 @@ SUBROUTINE InitBV() BiasVoltage%Frequency = GETReal('BiasVoltage-Frequency') BiasVoltage%Delta = GETReal('BiasVoltage-Delta') #if USE_LOADBALANCE -! Do not set during load balance in order to keep the old value -IF(.NOT.PerformLoadBalance)& +! Do not nullify during load balance in order to keep the old value on the MPIRoot +IF((.NOT.PerformLoadBalance).OR.(.NOT.MPIRoot))THEN #endif /*USE_LOADBALANCE*/ BiasVoltage%BVData = 0. + ! Update time + IF(BiasVoltage%Frequency.GT.0.0) BiasVoltage%BVData(3) = 1.0/BiasVoltage%Frequency +#if USE_LOADBALANCE +END IF +#endif /*USE_LOADBALANCE*/ IF(BiasVoltage%NPartBoundaries.LT.1) CALL CollectiveStop(__STAMP__,' UseBiasVoltage=T requires one or more particle boundaries!') @@ -1316,7 +1315,6 @@ SUBROUTINE InitBV() IF(MPIRoot)THEN BConProc = .TRUE. ELSE - BiasVoltage%BVData = 0. ! Check local sides DO SideID=1,nBCSides iBC = BC(SideID) diff --git a/src/hdg/hdg_vars.f90 b/src/hdg/hdg_vars.f90 index 9ccc3f015..74e98b1be 100644 --- a/src/hdg/hdg_vars.f90 +++ b/src/hdg/hdg_vars.f90 @@ -261,7 +261,7 @@ MODULE MOD_HDG_Vars #endif /*USE_MPI*/ INTEGER :: NPartBoundaries !< Total number of boundaries where the particles are counted INTEGER,ALLOCATABLE :: PartBoundaries(:) !< Part-boundary number on which the particles are counted - REAL :: Frequency !< Frequency with which the bias voltage is adjusted (every period T = 1/f the bias voltage is changed) + REAL :: Frequency !< Adaption nrequency with which the bias voltage is adjusted (every period T = 1/f the bias voltage is changed) REAL :: Delta !< Voltage difference used to change the current bias voltage (may also be adjusted over time automatically) REAL :: BVData(BVDataLength) !< 1: bias voltage ! !< 2: Ion excess diff --git a/src/io_hdf5/hdf5_output_state.f90 b/src/io_hdf5/hdf5_output_state.f90 index 7b5a12852..dc6b85ee7 100644 --- a/src/io_hdf5/hdf5_output_state.f90 +++ b/src/io_hdf5/hdf5_output_state.f90 @@ -171,7 +171,7 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) INTEGER :: i,j,k,iElem #endif /*PARTICLES*/ REAL,ALLOCATABLE :: FPCDataHDF5(:,:),EPCDataHDF5(:,:),BVDataHDF5(:,:) -INTEGER :: nVarFPC,nVarEPC,nVarBV +INTEGER :: nVarFPC,nVarEPC #endif /*USE_HDG*/ !=================================================================================================================================== #ifdef EXTRAE @@ -593,15 +593,14 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) #if USE_HDG ! Bias voltage IF(UseBiasVoltage.AND.MPIRoot)THEN - nVarBV = BVDataLength - ALLOCATE(BVDataHDF5(1:nVarBV,1)) + ALLOCATE(BVDataHDF5(BVDataLength,1)) CALL OpenDataFile(FileName,create=.FALSE.,single=.TRUE.,readOnly=.FALSE.) - BVDataHDF5(1:nVarBV,1) = BiasVoltage%BVData(1:nVarBV) + BVDataHDF5(BVDataLength,1) = BiasVoltage%BVData(BVDataLength) CALL WriteArrayToHDF5( DataSetName = 'BiasVoltage' , rank = 2 , & - nValGlobal = (/1_IK , INT(nVarBV,IK)/), & - nVal = (/1_IK , INT(nVarBV,IK)/), & + nValGlobal = (/1_IK , INT(BVDataLength,IK)/), & + nVal = (/1_IK , INT(BVDataLength,IK)/), & offset = (/0_IK , 0_IK/) , & - collective = .FALSE., RealArray = BVDataHDF5(1:nVarBV,1)) + collective = .FALSE., RealArray = BVDataHDF5(BVDataLength,1)) CALL CloseDataFile() DEALLOCATE(BVDataHDF5) END IF ! CalcBulkElectronTempi.AND.MPIRoot diff --git a/src/particles/surfacemodel/surfacemodel_analyze.f90 b/src/particles/surfacemodel/surfacemodel_analyze.f90 index f8f4a8a40..673405bb8 100644 --- a/src/particles/surfacemodel/surfacemodel_analyze.f90 +++ b/src/particles/surfacemodel/surfacemodel_analyze.f90 @@ -375,30 +375,30 @@ SUBROUTINE AnalyzeSurface(Time) END IF ! CalcElectronSEE END DO ! iBoundary = 1, BiasVoltage%NPartBoundaries - ! Add total electric charge (without displacement current!) - BiasVoltage%BVData(2) = BiasVoltage%BVData(2) + TotalElectricCharge - - ! Simulation time threshold - IF(time.GE.BiasVoltage%BVData(3))THEN - ! Update time - IF(BiasVoltage%Frequency.GT.0.0)THEN - BiasVoltage%BVData(3) = BiasVoltage%BVData(3) + 1.0/BiasVoltage%Frequency - END IF ! BiasVoltage%Frequency.GT.0.0 - - ! Update Voltage - IF(BiasVoltage%BVData(2).GT.0.0)THEN - ! Increase voltage - BiasVoltage%BVData(1) = BiasVoltage%BVData(1) + BiasVoltage%Delta - ELSEIF(BiasVoltage%BVData(2).LT.0.0)THEN - ! Decrease voltage - BiasVoltage%BVData(1) = BiasVoltage%BVData(1) - BiasVoltage%Delta - ELSE - ! do nothing - END IF ! BiasVoltage%BVData(2).LT.0.0 - - ! Reset ion excess counter - BiasVoltage%BVData(2) = 0. - END IF ! time.GE.BiasVoltage%BVData(3) + ASSOCIATE( V => BiasVoltage%BVData(1), Q => BiasVoltage%BVData(2), tBV => BiasVoltage%BVData(3) ) + ! Add total electric charge (without displacement current!) + Q = Q + TotalElectricCharge + + ! Simulation time threshold + IF(time.GE.tBV)THEN + ! Update time + IF(BiasVoltage%Frequency.GT.0.0) tBV = tBV + 1.0/BiasVoltage%Frequency + + ! Update Voltage + IF(Q.GT.0.0)THEN + ! Increase voltage + V = V + BiasVoltage%Delta + ELSEIF(Q.LT.0.0)THEN + ! Decrease voltage + V = V - BiasVoltage%Delta + ELSE + ! do nothing + END IF ! Q.LT.0.0 + + ! Reset ion excess counter + Q = 0. + END IF ! time.GE.tBV + END ASSOCIATE ! Write: Voltage, Ion excess and simulation update time DO i = 1, 3 From 922594805cbbf05d269f6fda986c1223be535d6e Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 19 Apr 2023 16:15:44 +0200 Subject: [PATCH 05/11] Added user guide description for bias voltage (currently only pure DC adjustment implemented). --- .../features-and-models/BC-field-solver.md | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/docs/documentation/userguide/features-and-models/BC-field-solver.md b/docs/documentation/userguide/features-and-models/BC-field-solver.md index 6385fcd3b..b5bdc1c7d 100644 --- a/docs/documentation/userguide/features-and-models/BC-field-solver.md +++ b/docs/documentation/userguide/features-and-models/BC-field-solver.md @@ -96,6 +96,9 @@ as detailed in the following table. | (/11,0/) | Neumann | q*n=1 | | | | | | (/20,1/) | FPC | 1: Assign BC to FPC group nbr. 1 (different BCs can be assigned the same FPC), see {ref}`sec:floating-boundary-condition` | +| | | | +| (/50,0/) | Dirichlet | Bias voltage for DC BCs (0: this number has no meaning). For more details, see {ref}`sec:bias-voltage-for-dc` | +| | | | ### RefState boundaries {-} @@ -233,6 +236,42 @@ To selected the direction by hand, simply supply the desired direction via with 1: x-, 2: y-, 3: z-direction. +(sec:bias-voltage-for-dc)= +### Bias Voltage: DC boundary +A bias voltage develops if charged particles impact on an electrode that is connected via capacitor to a matching box. Hence, if +there is an overhead of electrons impacting on the electrode, a negative bias voltage arises, which counters the flow of electrons +to that boundary. This feature only works when `PARTICLES=ON` is enabled and Poisson's equation is solved `PICLAS_EQNSYSNAME=poisson`. +The following parameters are required + + UseBiasVoltage = T ! Activate bias voltage model + BiasVoltage-NPartBoundaries = 2 ! Nbr of particle boundaries where the ion excess is measured + Biasvoltage-PartBoundaries = (/1,2/) ! Particle boundary index of where the ion excess is measured + BiasVoltage-Frequency = 3e9 ! Frequency with which the bias voltage is adapted + BiasVoltage-Delta = 1.0 ! Voltage difference used for adapting the bias voltage + +where `UseBiasVoltage` switches the modell on/off, `BiasVoltage-NPartBoundaries` defines the number of boundaries over which the +total electric charge is summarized, `Biasvoltage-PartBoundaries` specifies the indices of the particle boundaries where the total +electric charge is summarized, e.g., the powered and the grounded electrode., `BiasVoltage-Frequency` is the frequency with which +the bias voltage is updated by the voltage difference `BiasVoltage-Delta`. + +Additionally, the field boundary must be defined + + BoundaryName = BC_left + BoundaryType = (/50,0/) ! Dirichlet with 0V initial BC + +where the value *50* is used for for pure DC boundaries. +Furthermore, the bias voltage model relies on the functionality of the integral particle flux measurement on the corresponding +boundaries, see *BoundaryParticleOutput (BPO)* in Section {ref}`sec:integral-surface-variables`. +The definition of BPO variables must include at least the ones above, e.g. + + CalcBoundaryParticleOutput = T + BPO-NPartBoundaries = 2 + BPO-PartBoundaries = (/1,2/) + BPO-NSpecies = 3 + BPO-Species = (/2,3,4/) + +where the defined species information must consider all relevant charged particle species that affect the bias voltage. + ## Dielectric Materials Dielectric material properties can be considered by defining regions (or specific elements) From 4be09a45b35411d9215399d92d32ae17f7debdab Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 19 Apr 2023 18:58:14 +0200 Subject: [PATCH 06/11] Added bias voltage functionality for sinusoidal AC excitation boundaries. --- .../turner_bias-voltage_DC/analyze.ini | 4 +- .../turner_bias-voltage_DC/externals.ini | 12 +++--- .../turner_bias-voltage_DC/parameter.ini | 20 ++++++--- .../{pre-hopr => pre-hopr-AC}/hopr.ini | 12 +++--- .../pre-hopr-DC/hopr.ini | 43 +++++++++++++++++++ src/analyze/analyzefield.f90 | 8 +++- src/equations/poisson/equation.f90 | 13 ++++-- src/hdg/hdg.f90 | 9 +++- src/piclas.h | 2 +- 9 files changed, 95 insertions(+), 28 deletions(-) rename regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/{pre-hopr => pre-hopr-AC}/hopr.ini (87%) create mode 100644 regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-DC/hopr.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini index 757b7efac..205f85175 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini @@ -8,7 +8,7 @@ compare_data_file_max_differences = 2 ! integrate columns x:y in a data file as integral(y(x), x, x(1), x(end)) integrate_line_file = SurfaceAnalyze.csv ! data file name integrate_line_columns = 0:9 ! columns x:y -integrate_line_integral_value = -1.0172 ! bias voltage over time -integrate_line_tolerance_value = 1.0e-3 ! tolerance +integrate_line_integral_value = -1.0 ! integrated bias voltage over time (time is normalised to [0,1]) +integrate_line_tolerance_value = 1.0e-4 ! tolerance integrate_line_tolerance_type = relative ! absolute or relative comparison integrate_line_multiplier = 1e9 ! 1/tEnd diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini index aa8f4729f..e9d673cb3 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini @@ -1,6 +1,8 @@ ! --- Externals Tool Reggie -MPI = 1 -externalbinary = ./hopr/build/bin/hopr -externaldirectory = pre-hopr -externalruntime = pre -cmd_pre_execute = ln\s-s\s./pre-hopr/turner2013_mesh.h5\s../. +MPI = 1 , 1 +externalbinary = ./hopr/build/bin/hopr , ./hopr/build/bin/hopr +externaldirectory = pre-hopr-DC , pre-hopr-AC +externalruntime = pre , pre +cmd_pre_execute = ln\s-s\s./pre-hopr-DC/turner2013_BCType50_DC_mesh.h5\s../. , ln\s-s\s./pre-hopr-AC/turner2013_BCType51_AC_mesh.h5\s../. + +nocrosscombination:MPI,externalbinary,externaldirectory,externalruntime,cmd_pre_execute diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini index c2190c375..5a659123e 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini @@ -10,7 +10,7 @@ NAnalyze = 10 ! Number of analyze points ! =============================================================================== ! ! MESH ! =============================================================================== ! -MeshFile = ./turner2013_mesh.h5 +MeshFile = turner2013_BCType50_DC_mesh.h5, turner2013_BCType51_AC_mesh.h5 useCurveds = F ! =============================================================================== ! ! OUTPUT / VISUALIZATION @@ -52,15 +52,25 @@ DoInitialAutoRestart = T Load-DeviationThreshold = 1e-3 LoadBalanceMaxSteps = 10 -ManualTimestep = 4.58E-11 +ManualTimestep = 3.00E-11 IterDisplayStep = 10 tend = 1.0E-9 Analyze_dt = 5.0E-10 ! =============================================================================== ! ! Particle Boundaries ! =============================================================================== ! -BoundaryName = BC_left -BoundaryType = (/50,0/) ! Dirichlet with 0V initial BC +! This BC is defined in hopr.ini +!BoundaryName = BC_left +!BoundaryType = (/50,0/),(/51,0/) ! Dirichlet with 0V initial BC, either DC or AC + +! The refstate for the AC case must be defined here +!RefState = (/150.0 , 13.56E6 , -1.57079632679/) ! RefState Nbr 1: Voltage, Frequency and Phase shift +RefState = (/10.0 , 3e9 , -1.57079632679/) ! RefState Nbr 1: Voltage, Frequency and Phase shift + +! Boundary Field Output: Write current voltage to FieldAnalyze.csv +CalcBoundaryFieldOutput = T +BFO-NFieldBoundaries = 1 ! Nbr of boundaries +BFO-FieldBoundaries = (/5/) ! Field-Boundary1 BoundaryName = BC_right BoundaryType = (/4,0/) ! 4: Dirichlet with zero potential @@ -69,7 +79,7 @@ BoundaryType = (/4,0/) ! 4: Dirichlet with zero potential UseBiasVoltage = T BiasVoltage-NPartBoundaries = 2 ! Nbr of particle boundaries where the ion excess is measured Biasvoltage-PartBoundaries = (/1,2/) ! Particle boundary index of where the ion excess is measured -BiasVoltage-Frequency = 3e9 ! Frequency with which the bias voltage is adapted +BiasVoltage-Frequency = 3e9 ! Frequency with which the bias voltage is adapted BiasVoltage-Delta = 1.0 ! Voltage difference used for adapting the bias voltage ! =============================================================================== ! ! PARTICLES diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr/hopr.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-AC/hopr.ini similarity index 87% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr/hopr.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-AC/hopr.ini index 8c120e058..7346f6f95 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr/hopr.ini +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-AC/hopr.ini @@ -1,9 +1,9 @@ !=============================================================================== ! ! OUTPUT !=============================================================================== ! -ProjectName = turner2013 ! name of the project (used for filenames) -Debugvisu =T ! Write debug mesh to tecplot file -Logging =F ! Write log files +ProjectName = turner2013_BCType51_AC ! name of the project (used for filenames) +Debugvisu = T ! Write debug mesh to tecplot file +Logging = F ! Write log files !=============================================================================== ! ! MESH @@ -23,8 +23,7 @@ MeshScale = 1e-2 !=============================================================================== ! ! BOUNDARY CONDITIONS !=============================================================================== ! - nUserDefinedBoundaries=6 - +nUserDefinedBoundaries=6 BoundaryName=BC_periodicz- ! BC index 1 (from position in parameterfile) BoundaryType=(/1,0,0,1/) ! (/ Type, curveIndex, State, alpha /) @@ -37,9 +36,8 @@ BoundaryType=(/1,0,0,2/) BoundaryName=BC_periodicy+ ! BC index 4 BoundaryType=(/1,0,0,-2/) ! (/ BCType=1: periodic, 0, 0, Index of second vector vv in parameter file /) vv=(/0.,3.42e-3,0./) ! vector for periodic BC in y direction (yminus,yplus), index=2 - BoundaryName=BC_left -BoundaryType=(/4,0,0,0/) ! ideal conductor +BoundaryType=(/51,0,1,0/) ! ideal conductor BoundaryName=BC_right BoundaryType=(/4,0,0,0/) ! ideal conductor diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-DC/hopr.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-DC/hopr.ini new file mode 100644 index 000000000..449e22b3d --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-DC/hopr.ini @@ -0,0 +1,43 @@ +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! +ProjectName = turner2013_BCType50_DC ! name of the project (used for filenames) +Debugvisu = T ! Write debug mesh to tecplot file +Logging = F ! Write log files + +!=============================================================================== ! +! MESH +!=============================================================================== ! +Mode =1 ! 1 Cartesian 2 gambit file 3 CGNS +nZones =1 ! number of zones +Corner =(/0.,0.,0.0,,6.7,0.,0.0,,6.7,3.42e-3,0.0,,0.,3.42e-3,0.0 ,,0.,0.,3.42e-3,,6.7,0.,3.42e-3,,6.7,3.42e-3,3.42e-3,,0.,3.42e-3,3.42e-3/) ! [-3,3]x[-3,3]x[-3,3] +nElems =(/512,1,1/) ! Anzahl der Elemente in jede Richtung +BCIndex =(/1,3,6,4,5,2/) ! Indices of Boundary Conditions for six Boundary Faces (z-,y-,x+,y+,x-,z+) +elemtype =108 ! Elementform (108: Hexaeder) +useCurveds =F ! T if curved boundaries defined +SpaceQuandt =1. ! characteristic length of the mesh +ConformConnect=T + +postScaleMesh = T +MeshScale = 1e-2 +!=============================================================================== ! +! BOUNDARY CONDITIONS +!=============================================================================== ! +nUserDefinedBoundaries=6 + +BoundaryName=BC_periodicz- ! BC index 1 (from position in parameterfile) +BoundaryType=(/1,0,0,1/) ! (/ Type, curveIndex, State, alpha /) +BoundaryName=BC_periodicz+ ! BC index 2 +BoundaryType=(/1,0,0,-1/) ! here the direction of the vector 1 is changed, because it is the opposite side +vv=(/0.,0.,3.42e-3/) ! vector for periodic BC in z direction (zminus,zplus), index=1 + +BoundaryName=BC_periodicy- ! BC index 3 +BoundaryType=(/1,0,0,2/) +BoundaryName=BC_periodicy+ ! BC index 4 +BoundaryType=(/1,0,0,-2/) ! (/ BCType=1: periodic, 0, 0, Index of second vector vv in parameter file /) +vv=(/0.,3.42e-3,0./) ! vector for periodic BC in y direction (yminus,yplus), index=2 + +BoundaryName=BC_left +BoundaryType=(/50,0,0,0/) ! ideal conductor +BoundaryName=BC_right +BoundaryType=(/4,0,0,0/) ! ideal conductor diff --git a/src/analyze/analyzefield.f90 b/src/analyze/analyzefield.f90 index d155a8899..497d31b62 100644 --- a/src/analyze/analyzefield.f90 +++ b/src/analyze/analyzefield.f90 @@ -1834,8 +1834,12 @@ SUBROUTINE CalculateBoundaryFieldOutput(iBC,Time,BoundaryFieldOutput) CALL ExactFunc(BCState , x , BoundaryFieldOutput , t=Time) CASE(4) ! exact BC = Dirichlet BC !! Zero potential BoundaryFieldOutput = 0. - CASE(5) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional) - CALL ExactFunc( -1 , x , BoundaryFieldOutput , t=Time , iRefState=BCState) + CASE(5,51) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional) + IF(BCType.EQ.51)THEN + CALL ExactFunc( 51 , x , BoundaryFieldOutput , t=Time , iRefState=BCState) + ELSE + CALL ExactFunc( -1 , x , BoundaryFieldOutput , t=Time , iRefState=BCState) + END IF ! BCType.EQ.51 CASE(6) ! exact BC = Dirichlet BC !! ExactFunc via RefState (Time is optional) CALL ExactFunc( -2 , x , BoundaryFieldOutput , t=time , iRefState=BCState) CASE(50) ! exact BC = Dirichlet BC !! ExactFunc via bias voltage DC diff --git a/src/equations/poisson/equation.f90 b/src/equations/poisson/equation.f90 index 4c61d4d40..cfec62b02 100644 --- a/src/equations/poisson/equation.f90 +++ b/src/equations/poisson/equation.f90 @@ -130,6 +130,8 @@ SUBROUTINE InitEquation() CHARACTER(LEN=255) :: BCName INTEGER :: nRefStateMax INTEGER :: nLinState,nLinStateMax +INTEGER,PARAMETER :: BYTypeRefstate(1:2)=(/5,51/) +CHARACTER(LEN=32) :: hilf !=================================================================================================================================== IF((.NOT.InterpolationInitIsDone).OR.EquationInitIsDone)THEN LBWRITE(*,*) "InitPoisson not ready to be called or already called." @@ -149,14 +151,15 @@ SUBROUTINE InitEquation() BCType = BoundaryType(i,BC_TYPE) BCState = BoundaryType(i,BC_STATE) BCName = BoundaryName(i) - IF(BCType.EQ.5.AND.BCState.LE.0) THEN + IF(ANY(BCType.EQ.BYTypeRefstate).AND.BCState.LE.0) THEN SWRITE(*,'(A)') "Error found for the following boundary condition" SWRITE(*,'(A,I0)') " BC: ",i SWRITE(*,'(A,I0)') " Type: ",BCType SWRITE(*,'(A,I0)') "State: ",BCState SWRITE(*,'(A)') " Name: "//TRIM(BCName) - CALL abort(__STAMP__,'BCState is <= 0 for BCType=5 is not allowed! Set a positive integer for the n-th RefState') - ELSEIF(BCType.EQ.5.AND.BCState.GT.0)THEN + WRITE(UNIT=hilf,FMT='(I0)') BCType + CALL abort(__STAMP__,'BCState is <= 0 for BCType='//TRIM(hilf)//' is not allowed! Set a positive integer for the n-th RefState') + ELSEIF(ANY(BCType.EQ.BYTypeRefstate).AND.BCState.GT.0)THEN nRefStateMax = MAX(nRefStateMax,BCState) ELSEIF(BCType.EQ.7.AND.BCState.GT.0)THEN nLinStateMax = MAX(nLinStateMax,BCState) @@ -527,12 +530,14 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState,iLinState,BCState) Omega = 2.*PI*RefState(2,iRefState) r1 = RefState(1,iRefState) / 2.0 Resu(:) = r1*(COS(Omega*t+RefState(3,iRefState)) + 1.0) -CASE(-1) ! Signal with zero-crossing: Amplitude, Frequency and Phase Shift supplied by RefState +CASE(-1,51) ! Signal with zero-crossing: Amplitude, Frequency and Phase Shift supplied by RefState ! RefState(1,iRefState): amplitude ! RefState(2,iRefState): frequency ! RefState(3,iRefState): phase shift Omega = 2.*PI*RefState(2,iRefState) Resu(:) = RefState(1,iRefState)*COS(Omega*t+RefState(3,iRefState)) + ! Add bias potential (only if bias voltage model is activated) + IF(ExactFunction.EQ.51) Resu(:) = Resu(:) + BiasVoltage%BVData(1) CASE(0) ! constant 0. Resu(:)=0. #if defined(PARTICLES) diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index cd4961eac..26176df72 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -1954,6 +1954,11 @@ SUBROUTINE HDGLinear(time,U_out) r=q*(PP_N+1) + p+1 CALL ExactFunc(-5,Face_xGP(:,p,q,SideID),lambda(PP_nVar,r:r,SideID),t=time,BCState=BCState) END DO; END DO !p,q + CASE(51) ! exact BC = Dirichlet BC !! ExactFunc via DC bias voltage + DO q=0,PP_N; DO p=0,PP_N + r=q*(PP_N+1) + p+1 + CALL ExactFunc(51,Face_xGP(:,p,q,SideID),lambda(PP_nVar,r:r,SideID),t=time,iRefState=BCState) + END DO; END DO !p,q END SELECT ! BCType END DO !BCsideID=1,nDirichletBCSides #if (PP_nVar!=1) @@ -2394,8 +2399,8 @@ SUBROUTINE HDGNewton(time,U_out,td_iter,ForceCGSolverIteration_opt) r=q*(PP_N+1) + p+1 CALL ExactFunc(-3,Face_xGP(:,p,q,SideID),lambda(PP_nVar,r:r,SideID),t=time,iLinState=BCState) END DO; END DO !p,q - CASE(8,50) ! exact BC = Dirichlet BC !! ExactFunc via electric potential and decharing of a surface - CALL abort(__STAMP__,'Dirichlet BC=8,50 model not implemented for HDG Newton!') + CASE(8,50,51) ! exact BC = Dirichlet BC !! ExactFunc via electric potential and decharing of a surface + CALL abort(__STAMP__,'Dirichlet BC=8,50,51 model not implemented for HDG Newton!') END SELECT ! BCType END DO !BCsideID=1,nDirichletBCSides diff --git a/src/piclas.h b/src/piclas.h index e10bc2548..4c6a7ed94 100644 --- a/src/piclas.h +++ b/src/piclas.h @@ -324,5 +324,5 @@ #if USE_HDG ! HDG Dirichlet BC Side IDs -#define HDGDIRICHLETBCSIDEIDS 2,4,5,6,7,8,50 +#define HDGDIRICHLETBCSIDEIDS 2,4,5,6,7,8,50,51 #endif From 7185457f7fc6bceb7fd45443ba23f6d4a28a37fe Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 19 Apr 2023 19:01:31 +0200 Subject: [PATCH 07/11] Renamed bias voltage reggie, whcih now also tests bias voltage in combination with AC (BCType=51) --- REGGIE.md | 2 +- .../DSMC.ini | 0 .../DSMCSpecies_electronic_state_full_Data.h5 | Bin .../PartAnalyze_ref.csv | 0 .../analyze.ini | 0 .../command_line.ini | 0 .../excludeBuild.ini | 0 .../externals.ini | 0 .../parameter.ini | 0 .../pre-hopr-AC/hopr.ini | 0 .../pre-hopr-DC/hopr.ini | 0 .../readme.md | 4 ++-- 12 files changed, 3 insertions(+), 3 deletions(-) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/DSMC.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/DSMCSpecies_electronic_state_full_Data.h5 (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/PartAnalyze_ref.csv (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/analyze.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/command_line.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/excludeBuild.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/externals.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/parameter.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/pre-hopr-AC/hopr.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/pre-hopr-DC/hopr.ini (100%) rename regressioncheck/NIG_PIC_poisson_RK3/{turner_bias-voltage_DC => turner_bias-voltage_AC-DC}/readme.md (63%) diff --git a/REGGIE.md b/REGGIE.md index 7d7f72055..8097e2de8 100644 --- a/REGGIE.md +++ b/REGGIE.md @@ -338,7 +338,7 @@ Testing PIC compiled with Runge-Kutta 3 integration, solving Poisson's equation: | 6 | plasma_sheath_BR-electrons_conforming_auto-switch_variable_Te | | non-linear HDG (BR electrons), automatic switching BR/kinetic, variable Te, change nSkipAnalyze during the simulation | nProcs=1,2,4,11 | integrate Te over time (PartAnalyze.csv) | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_conforming_auto-switch_variable_Te/readme.md) | | 7 | plasma_sheath_BR-electrons_mortar | | non-linear HDG (BR electrons), Mortars | nProcs=2 | TimeAvg | [Link](regressioncheck/NIG_PIC_poisson_RK3-electrons_conforming/plasma_sheath_BR-electrons_mortar/readme.md) | | 8 | turner | | | nProcs=4 | L2 error, PartAnalyze.csv | | -| 9 | turner_bias-voltage_DC | | bias voltage for DC potential boundaries (BCType=50) | nProcs=1,2,4,10 | PartAnalyze.csv, SurfaceAnalyze.csv | [Link](regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md) | +| 9 | turner_bias-voltage_AC-DC | | bias voltage for AC and DC potential boundaries (BCType=50) | nProcs=1,2,4,10 | PartAnalyze.csv, SurfaceAnalyze.csv | [Link](regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/readme.md) | ### NIG_PIC_maxwell_RK4 diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMC.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/DSMC.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMC.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/DSMC.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMCSpecies_electronic_state_full_Data.h5 b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/DSMCSpecies_electronic_state_full_Data.h5 similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/DSMCSpecies_electronic_state_full_Data.h5 rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/DSMCSpecies_electronic_state_full_Data.h5 diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/PartAnalyze_ref.csv b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/PartAnalyze_ref.csv similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/PartAnalyze_ref.csv rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/PartAnalyze_ref.csv diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/analyze.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/analyze.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/analyze.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/command_line.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/command_line.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/command_line.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/command_line.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/excludeBuild.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/excludeBuild.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/excludeBuild.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/excludeBuild.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/externals.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/externals.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/externals.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/parameter.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/parameter.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/parameter.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-AC/hopr.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/pre-hopr-AC/hopr.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-AC/hopr.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/pre-hopr-AC/hopr.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-DC/hopr.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/pre-hopr-DC/hopr.ini similarity index 100% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/pre-hopr-DC/hopr.ini rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/pre-hopr-DC/hopr.ini diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/readme.md similarity index 63% rename from regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md rename to regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/readme.md index 1fb6ef241..6a6300fb4 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_DC/readme.md +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/readme.md @@ -1,5 +1,5 @@ # 1D turner setup and bias voltage feature -- He plasma between two electrodes one of which is usually AC, but for this bias-voltage test case it is set to DC starting at 0V -- bias voltage for DC only (field BCType=50) +- He plasma between two electrodes one of which is usually AC, but for this bias-voltage test case it is set to AC or DC starting at 0V +- bias voltage for AC or DC only (field BCType=50 or 51, the latter requiring a refstate to define the sinusoidal function) - ion excess on both electrodes generates a positive bias voltage and electron excess a negative potential - comparison of PartAnalyze.csv with reference file and integration of the bias voltage over time from SurfaceAnalyze.csv From 83f371bfc449468c33c8a96e490589bf72a2daac Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Wed, 19 Apr 2023 23:43:16 +0200 Subject: [PATCH 08/11] Fixed bias voltage compilation with PARTICLES=OFF, MPI=OFF etc. --- src/equations/poisson/equation.f90 | 12 +++++----- src/hdg/hdg.f90 | 12 ++++++---- src/io_hdf5/hdf5_output_state.f90 | 11 +++++---- .../surfacemodel/surfacemodel_analyze.f90 | 2 ++ src/restart/restart_field.f90 | 24 ++++++++++--------- 5 files changed, 35 insertions(+), 26 deletions(-) diff --git a/src/equations/poisson/equation.f90 b/src/equations/poisson/equation.f90 index cfec62b02..14fafc547 100644 --- a/src/equations/poisson/equation.f90 +++ b/src/equations/poisson/equation.f90 @@ -43,16 +43,14 @@ MODULE MOD_Equation INTERFACE FinalizeEquation MODULE PROCEDURE FinalizeEquation END INTERFACE -#if USE_MPI + +PUBLIC :: InitEquation,ExactFunc,CalcSource,FinalizeEquation, CalcSourceHDG,DivCleaningDamping +#if USE_MPI && defined(PARTICLES) INTERFACE SynchronizeCPP MODULE PROCEDURE SynchronizeCPP END INTERFACE -#endif /*USE_MPI*/ - -PUBLIC :: InitEquation,ExactFunc,CalcSource,FinalizeEquation, CalcSourceHDG,DivCleaningDamping -#if USE_MPI PUBLIC :: SynchronizeCPP -#endif /*USE_MPI*/ +#endif /*USE_MPI && defined(PARTICLES)*/ !=================================================================================================================================== PUBLIC::DefineParametersEquation CONTAINS @@ -536,8 +534,10 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState,iLinState,BCState) ! RefState(3,iRefState): phase shift Omega = 2.*PI*RefState(2,iRefState) Resu(:) = RefState(1,iRefState)*COS(Omega*t+RefState(3,iRefState)) +#if defined(PARTICLES) ! Add bias potential (only if bias voltage model is activated) IF(ExactFunction.EQ.51) Resu(:) = Resu(:) + BiasVoltage%BVData(1) +#endif /*defined(PARTICLES)*/ CASE(0) ! constant 0. Resu(:)=0. #if defined(PARTICLES) diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index 26176df72..0b1c74e6e 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -40,8 +40,11 @@ MODULE MOD_HDG PUBLIC :: HDG, RestartHDG PUBLIC :: DefineParametersHDG #if USE_MPI -PUBLIC :: SynchronizeChargeOnFPC,SynchronizeVoltageOnEPC,SynchronizeBV -#endif /*USE_MPI*/ +PUBLIC :: SynchronizeChargeOnFPC,SynchronizeVoltageOnEPC +#if defined(PARTICLES) + PUBLIC :: SynchronizeBV +#endif /*defined(PARTICLES)*/ +#endif /*USE_MPI */ #endif /*USE_HDG*/ !=================================================================================================================================== @@ -1220,7 +1223,7 @@ END SUBROUTINE InitEPC !=================================================================================================================================== SUBROUTINE InitBV() ! MODULES -USE MOD_Globals ,ONLY: CollectiveStop +USE MOD_Globals ,ONLY: CollectiveStop,UNIT_StdOut USE MOD_ReadInTools ,ONLY: GETLOGICAL,GETREAL,GETINT,GETINTARRAY USE MOD_Particle_Boundary_Vars ,ONLY: PartBound USE MOD_Mesh_Vars ,ONLY: nBCs,BoundaryType,BoundaryName @@ -1231,7 +1234,6 @@ SUBROUTINE InitBV() #endif /*USE_LOADBALANCE*/ #if USE_MPI USE MOD_Globals ,ONLY: IERROR,MPI_COMM_NULL,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,MPI_INFO_NULL,MPI_UNDEFINED,MPIRoot -USE MOD_Globals ,ONLY: UNIT_StdOut USE MOD_Mesh_Vars ,ONLY: nBCSides,BC #endif /*USE_MPI*/ IMPLICIT NONE @@ -1413,6 +1415,7 @@ SUBROUTINE ReadBVDataFromH5() END SUBROUTINE ReadBVDataFromH5 +#if USE_MPI !=================================================================================================================================== !> Communicate the bias voltage values from MPIRoot to sub-communicator processes !=================================================================================================================================== @@ -1432,6 +1435,7 @@ SUBROUTINE SynchronizeBV() CALL MPI_BCAST(BiasVoltage%BVData, BVDataLength, MPI_DOUBLE_PRECISION, 0, BiasVoltage%COMM%UNICATOR, IERROR) END IF END SUBROUTINE SynchronizeBV +#endif /*USE_MPI*/ #endif /*defined(PARTICLES)*/ diff --git a/src/io_hdf5/hdf5_output_state.f90 b/src/io_hdf5/hdf5_output_state.f90 index dc6b85ee7..1663e56ad 100644 --- a/src/io_hdf5/hdf5_output_state.f90 +++ b/src/io_hdf5/hdf5_output_state.f90 @@ -72,7 +72,7 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) USE MOD_Equation_Vars ,ONLY: E,Phi #endif /*PP_POIS*/ #if USE_HDG -USE MOD_HDG_Vars ,ONLY: nGP_face,iLocSides,UseFPC,FPC,UseEPC,EPC,UseBiasVoltage,BiasVoltage,BVDataLength +USE MOD_HDG_Vars ,ONLY: nGP_face,iLocSides,UseFPC,FPC,UseEPC,EPC #if PP_nVar==1 USE MOD_Equation_Vars ,ONLY: E,Et #elif PP_nVar==3 @@ -93,6 +93,7 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) USE MOD_Mesh_Vars ,ONLY: GlobalUniqueSideID USE MOD_Analyze_Vars ,ONLY: CalcElectricTimeDerivative #ifdef PARTICLES +USE MOD_HDG_Vars ,ONLY: UseBiasVoltage,BiasVoltage,BVDataLength USE MOD_PICInterpolation_Vars ,ONLY: useAlgebraicExternalField,AlgebraicExternalField USE MOD_Analyze_Vars ,ONLY: AverageElectricPotential USE MOD_Mesh_Vars ,ONLY: Elem_xGP @@ -587,10 +588,6 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) collective = .FALSE., RealArray = TmpArray2(1:3,1)) CALL CloseDataFile() END IF ! CalcBulkElectronTempi.AND.MPIRoot -#endif /*USE_HDG*/ -#endif /*PARTICLES*/ - -#if USE_HDG ! Bias voltage IF(UseBiasVoltage.AND.MPIRoot)THEN ALLOCATE(BVDataHDF5(BVDataLength,1)) @@ -604,6 +601,10 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) CALL CloseDataFile() DEALLOCATE(BVDataHDF5) END IF ! CalcBulkElectronTempi.AND.MPIRoot +#endif /*USE_HDG*/ +#endif /*PARTICLES*/ + +#if USE_HDG ! Floating boundary condition: Store charge on each FPC IF(UseFPC.AND.MPIRoot)THEN nVarFPC = FPC%nUniqueFPCBounds diff --git a/src/particles/surfacemodel/surfacemodel_analyze.f90 b/src/particles/surfacemodel/surfacemodel_analyze.f90 index 673405bb8..5e5f15276 100644 --- a/src/particles/surfacemodel/surfacemodel_analyze.f90 +++ b/src/particles/surfacemodel/surfacemodel_analyze.f90 @@ -159,7 +159,9 @@ SUBROUTINE AnalyzeSurface(Time) USE MOD_Analyze_Vars ,ONLY: EDC USE MOD_Analyze_Vars ,ONLY: CalcElectricTimeDerivative USE MOD_HDG_Vars ,ONLY: UseBiasVoltage,BiasVoltage +#if USE_MPI USE MOD_HDG ,ONLY: SynchronizeBV +#endif /*USE_MPI*/ #endif /*USE_HDG*/ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE diff --git a/src/restart/restart_field.f90 b/src/restart/restart_field.f90 index ccb205d67..a60c016b2 100644 --- a/src/restart/restart_field.f90 +++ b/src/restart/restart_field.f90 @@ -74,10 +74,12 @@ SUBROUTINE FieldRestart() USE MOD_MPI ,ONLY: StartReceiveMPIData,StartSendMPIData,FinishExchangeMPIData #endif /*USE_MPI*/ #if USE_LOADBALANCE -USE MOD_Equation ,ONLY: SynchronizeCPP -USE MOD_HDG ,ONLY: SynchronizeVoltageOnEPC,SynchronizeBV -USE MOD_HDG_Vars ,ONLY: UseEPC,UseBiasVoltage,CalcPCouplElectricPotential +USE MOD_HDG ,ONLY: SynchronizeVoltageOnEPC +USE MOD_HDG_Vars ,ONLY: UseEPC #if defined(PARTICLES) +USE MOD_Equation ,ONLY: SynchronizeCPP +USE MOD_HDG ,ONLY: SynchronizeBV +USE MOD_HDG_Vars ,ONLY: UseBiasVoltage,CalcPCouplElectricPotential ! TODO: make ElemInfo available with PARTICLES=OFF and remove this preprocessor if/else as soon as possible USE MOD_Mesh_Vars ,ONLY: SideToNonUniqueGlobalSide,nElems USE MOD_LoadBalance_Vars ,ONLY: MPInSideSend,MPInSideRecv,MPIoffsetSideSend,MPIoffsetSideRecv @@ -153,14 +155,6 @@ SUBROUTINE FieldRestart() #if USE_LOADBALANCE IF(PerformLoadBalance.AND.(.NOT.UseH5IOLoadBalance))THEN #if USE_HDG - ! Coupled Bias voltage (BV): The MPI root process distributes the information among the sub-communicator processes - ! (before and after load balancing, the root process is always part of each sub-communicator group) - IF(UseBiasVoltage) CALL SynchronizeBV() - - ! Coupled Power Potential (CPP): The MPI root process distributes the information among the sub-communicator processes - ! (before and after load balancing, the root process is always part of each sub-communicator group) - IF(CalcPCouplElectricPotential) CALL SynchronizeCPP() - #if USE_PETSC ! FPC: The MPI root process distributes the information among the sub-communicator processes for each FPC ! (before and after load balancing, the root process is always part of each sub-communicator group) @@ -171,6 +165,14 @@ SUBROUTINE FieldRestart() IF(UseEPC) CALL SynchronizeVoltageOnEPC() #if defined(PARTICLES) + ! Coupled Bias voltage (BV): The MPI root process distributes the information among the sub-communicator processes + ! (before and after load balancing, the root process is always part of each sub-communicator group) + IF(UseBiasVoltage) CALL SynchronizeBV() + + ! Coupled Power Potential (CPP): The MPI root process distributes the information among the sub-communicator processes + ! (before and after load balancing, the root process is always part of each sub-communicator group) + IF(CalcPCouplElectricPotential) CALL SynchronizeCPP() + !checkRank=MIN(3,nProcessors-1) ! Store lambda solution on global non-unique array for MPI communication ASSOCIATE( firstSide => ElemInfo_Shared(ELEM_FIRSTSIDEIND,offsetElem+1) + 1 ,& From 70506de7187d6fc7e98920d5ee008ada6b36dc4b Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Thu, 20 Apr 2023 00:00:23 +0200 Subject: [PATCH 09/11] Added user guide description for bias voltage model + AC potential field boundary. --- .../features-and-models/BC-field-solver.md | 25 ++++++++++++++----- 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/docs/documentation/userguide/features-and-models/BC-field-solver.md b/docs/documentation/userguide/features-and-models/BC-field-solver.md index b5bdc1c7d..628d6e603 100644 --- a/docs/documentation/userguide/features-and-models/BC-field-solver.md +++ b/docs/documentation/userguide/features-and-models/BC-field-solver.md @@ -83,9 +83,9 @@ as detailed in the following table. | | | | | (/4,0/) | Dirichlet | zero-potential (Phi=0) | | | | | -| (/5,1/) | Dirichlet | 1: use RefState Nbr 1 and $\cos(\omega t)$ function (see details below) | +| (/5,1/) | Dirichlet | 1: use RefState Nbr 1 and $\cos(\omega t)$ function (for details see {ref}`sec:ref-state-bcs`) | | | | | -| (/6,1/) | Dirichlet | 1: use RefState Nbr 1 and $\cos(\omega t)+1$ function that does not switch sign (see details below) | +| (/6,1/) | Dirichlet | 1: use RefState Nbr 1 and $\cos(\omega t)+1$ function that does not switch sign (for details see {ref}`sec:ref-state-bcs`) | | | | | | (/7,1/) | Dirichlet | 1: use LinState Nbr 1, linear function for Phi, see {ref}`sec:linear-potential` | | | | | @@ -97,10 +97,12 @@ as detailed in the following table. | | | | | (/20,1/) | FPC | 1: Assign BC to FPC group nbr. 1 (different BCs can be assigned the same FPC), see {ref}`sec:floating-boundary-condition` | | | | | -| (/50,0/) | Dirichlet | Bias voltage for DC BCs (0: this number has no meaning). For more details, see {ref}`sec:bias-voltage-for-dc` | +| (/50,0/) | Dirichlet | Bias voltage for ***DC*** BCs (0: this number has no meaning). For more details, see {ref}`sec:bias-voltage-for-dc` | +| (/51,1/) | Dirichlet | Bias voltage for **AC** BCs (1: use RefState Nbr 1). For more details, see {ref}`sec:bias-voltage-for-ac` | | | | | -### RefState boundaries {-} +(sec:ref-state-bcs)= +### RefState boundaries For each boundary of type *5* (reference state boundary *RefState*), e.g., by setting the boundary in the *parameter.ini* file @@ -116,14 +118,14 @@ function (a frequency of 0 results in a fixed potential over time) and phase shi This yields the three parameters used in the cosine function - Phi(t) = A*COS(2*pi*f*t + psi) +$$\Phi(t)=A\cos (2\pi f t + \psi)$$ where *A=-0.18011* is the amplitude, *t* is the time, *f=1* is the frequency and *psi=0* is the phase shift. Similar to boundary type *5* is type *6*, which simply uses a cosine function that always has the same sign, depending on the amplitude *A* - Phi(t) = (A/2) * (COS(2*pi*f*t + psi) + 1) +$$\Phi(t)=\frac{A}{2}(\cos (2\pi f t + \psi)+1)$$ (sec:linear-potential)= ### Linear potential function @@ -272,6 +274,17 @@ The definition of BPO variables must include at least the ones above, e.g. where the defined species information must consider all relevant charged particle species that affect the bias voltage. +(sec:bias-voltage-for-ac)= +### Bias Voltage: AC boundary +If an AC potential boundary is coupled with a bias potential, the `BoundaryType` has to be changed as compared with the previous +section + + BoundaryName = BC_left + BoundaryType = (/51,1/) ! Dirichlet with 0V initial BC + +Furthermore, a RefState must be defined, which specifies the parameters for the $cos(\omega t)$ function, for details, see +{ref}`sec:ref-state-bcs`. Note that this BC is only implemented with zero crossing. + ## Dielectric Materials Dielectric material properties can be considered by defining regions (or specific elements) From af1087ee952985196f4c9b2a7afae141233aa653 Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Thu, 20 Apr 2023 11:11:48 +0200 Subject: [PATCH 10/11] Fixed small typo in read-in of bias voltage --- src/hdg/hdg.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index 0b1c74e6e..1f4a4b10c 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -1400,7 +1400,7 @@ SUBROUTINE ReadBVDataFromH5() ! Check for new parameter name IF(BVExists)THEN CALL ReadArray(TRIM(ContainerName) , 2 , (/1_IK , INT(BVDataLength,IK)/) , 0_IK , 1 , RealArray=BVDataHDF5) - WRITE(UNIT_stdOut,'(3(A,ES10.2E3))') " Read CoupledPowerPotential from restart file ["//TRIM(RestartFile)//& + WRITE(UNIT_stdOut,'(3(A,ES10.2E3))') " Read bias voltage from restart file ["//TRIM(RestartFile)//& "] Bias voltage[V]: ",BVDataHDF5(1),", Ion excess[C]: ",BVDataHDF5(2),", next adjustment time[s]: ",BVDataHDF5(3) BiasVoltage%BVData = BVDataHDF5 END IF ! BVExists From 4a353f2f28281089ad5b446b073e858f86194f6f Mon Sep 17 00:00:00 2001 From: Stephen Copplestone Date: Fri, 21 Apr 2023 08:20:57 +0200 Subject: [PATCH 11/11] Fixed bug in output of bias voltage data got .h5 state file. --- src/io_hdf5/hdf5_output_state.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io_hdf5/hdf5_output_state.f90 b/src/io_hdf5/hdf5_output_state.f90 index 1663e56ad..102393f2b 100644 --- a/src/io_hdf5/hdf5_output_state.f90 +++ b/src/io_hdf5/hdf5_output_state.f90 @@ -592,12 +592,12 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) IF(UseBiasVoltage.AND.MPIRoot)THEN ALLOCATE(BVDataHDF5(BVDataLength,1)) CALL OpenDataFile(FileName,create=.FALSE.,single=.TRUE.,readOnly=.FALSE.) - BVDataHDF5(BVDataLength,1) = BiasVoltage%BVData(BVDataLength) + BVDataHDF5(1:BVDataLength,1) = BiasVoltage%BVData(1:BVDataLength) CALL WriteArrayToHDF5( DataSetName = 'BiasVoltage' , rank = 2 , & nValGlobal = (/1_IK , INT(BVDataLength,IK)/), & nVal = (/1_IK , INT(BVDataLength,IK)/), & offset = (/0_IK , 0_IK/) , & - collective = .FALSE., RealArray = BVDataHDF5(BVDataLength,1)) + collective = .FALSE., RealArray = BVDataHDF5(1:BVDataLength,1)) CALL CloseDataFile() DEALLOCATE(BVDataHDF5) END IF ! CalcBulkElectronTempi.AND.MPIRoot