diff --git a/REGGIE.md b/REGGIE.md index 8097e2de8..05c85f360 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_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) | +| 9 | turner_bias-voltage_AC-DC | | bias voltage for AC with BCType=51 and 52 (power control) and DC with BCType=50 potential boundaries | 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/docs/documentation/userguide/features-and-models/BC-field-solver.md b/docs/documentation/userguide/features-and-models/BC-field-solver.md index 628d6e603..76affc808 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. -| (/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 | +| (/Type,State/) | 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,45 +61,47 @@ 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. -| (/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 (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 (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` | -| | | | -| (/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` | -| | | | -| (/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` | -| | | | +| (/Type,State/) | 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 (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 (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` | +| | | | +| (/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` | +| | | | +| (/50,0/) | Dirichlet | {ref}`sec:bias-voltage-for-dc` (0: this number has no meaning) | +| (/51,1/) | Dirichlet | {ref}`sec:bias-voltage-for-ac` (1: use RefState Nbr 1, see {ref}`sec:ref-state-bcs`) | +| (/52,1/) | Dirichlet | {ref}`sec:bias-voltage-for-ac-and-cpp` (1: use RefState Nbr 1, see {ref}`sec:ref-state-bcs`) and {ref}`sec:fixed-coupled-power`| +| | | | +| (/60,1/) | Dirichlet | {ref}`sec:fixed-coupled-power` for $\Phi=A\cos(\omega t)$ (1: use RefState Nbr 1, see {ref}`sec:ref-state-bcs`) | (sec:ref-state-bcs)= ### RefState boundaries @@ -203,16 +205,21 @@ The resulting current $I$ and voltage $U$ are automatically written to *FieldAna "008-EPC-Voltage-BCState-001". (sec:fixed-coupled-power)= -### Fixed coupled power (const. input power) +### Power control -An automatic adjustment of the electric potential to ensure that a fixed power input to the system is achieved requires the -following parameters +An automatic adjustment of a *DC* electric potential to ensure that a fixed power absorbed by the charged particles of the system is +achieved, requires the following parameters BoundaryName = BC_WALL ! BC name in the mesh.h5 file BoundaryType = (/2,2/) ! all BCs with this type will be adjusted to the same electric potential that is adjusted over time -Additionally, a starting value for the potential, lower and upper boundaries and a relaxation factor are required -as well as the target input power, which is set via +For an *AC* boundary the parameters are the following + + BoundaryName = BC_WALL ! BC name in the mesh.h5 file + BoundaryType = (/60,1/) ! applicable to one BC with this BCType. The BCState=1 corresponds to the RefState number + +Additionally, a starting value for the potential $\Phi$ (or amplitude $A$ in the case of *AC*), lower and upper boundaries and a +relaxation factor are required as well as the target input power, which is set via CoupledPowerPotential = (/10. , 1000. , 2000./) ! lower, starting and maximum values for the electric potential CoupledPowerRelaxFac = 0.05 ! the new potential is updated by 5% in each time step @@ -221,6 +228,19 @@ as well as the target input power, which is set via The values in `CoupledPowerPotential` correspond to the lower boundary, the starting value and the upper boundary, respectively. When a simulation is restarted from a state file, the last known value of the BC will be used instead of the starting value, which is only applied when starting a fresh simulation from scratch. +Note that in the case of an AC boundary, the amplitude defined in the `RefState` will be overwritten with the initial value given in +`CoupledPowerPotential`. + +There are three models implemented by which the power control adjusts the electric potential (or peak potential in case of AC) + + CoupledPowerMode = 3 ! 1: instantaneous power (default value), 2: moving average, 3: integrated and averaged over 1 cycle + CoupledPowerFrequency = 13.56e6 ! Frequency with which the coupled power voltage is adapte + +where `CoupledPowerMode` selects the model and `CoupledPowerFrequency` defines the adaption frequency and corresponding period over +which the power is integrated but only when model 3 is used. Model 1 uses the instantaneously determined absorbed power that inherently +oscillates, model 2 uses an averaged value (that is reset during every load balance or normal restart event) and model 3 integrated +the power linearly between particle analysis steps and adjusts the power after each period of the user-defined frequency +`CoupledPowerFrequency`. ### Zero potential enforcement @@ -239,7 +259,7 @@ 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 +### Bias Voltage DC 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`. @@ -258,7 +278,7 @@ the bias voltage is updated by the voltage difference `BiasVoltage-Delta`. Additionally, the field boundary must be defined - BoundaryName = BC_left + BoundaryName = BC_left ! BC name in the mesh.h5 file BoundaryType = (/50,0/) ! Dirichlet with 0V initial BC where the value *50* is used for for pure DC boundaries. @@ -275,16 +295,28 @@ 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 +### Bias Voltage AC 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 + BoundaryName = BC_left ! BC name in the mesh.h5 file 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. +(sec:bias-voltage-for-ac-and-cpp)= +### Bias Voltage AC and Fixed coupled power +If an **AC potential boundary** is coupled with a **bias potential** and an **automatic adjustment of the AC amplitude** to ensure +that a fixed power to the system is achieved, the `BoundaryType` has to be changed as compared with the previous section + + BoundaryName = BC_left ! BC name in the mesh.h5 file + BoundaryType = (/52,1/) ! Dirichlet with 0V initial BC + +where 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. +For details on the power coupling, see {ref}`sec:fixed-coupled-power`. + ## Dielectric Materials Dielectric material properties can be considered by defining regions (or specific elements) diff --git a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/externals.ini b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/externals.ini new file mode 100644 index 000000000..c69adb0a3 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/externals.ini @@ -0,0 +1,8 @@ +! --- Externals Tool Reggie +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/parallel_plates_BCType2_DC_mesh.h5\s../. , ln\s-s\s./pre-hopr-AC/parallel_plates_BCType60_AC_mesh.h5\s../. + +nocrosscombination:MPI,externalbinary,externaldirectory,externalruntime,cmd_pre_execute diff --git a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/hopr.ini b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/hopr.ini deleted file mode 100644 index 982651ec8..000000000 --- a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/hopr.ini +++ /dev/null @@ -1,65 +0,0 @@ -DEFVAR=(INT): iz = 1 ! no. of nFineHexa -!DEFVAR=(REAL): li = 6.7 ! length -DEFVAR=(REAL): li = 1 ! length -DEFVAR=(REAL): lx = 1 ! length -! MAKEFILE PARAMETER (put a "#" in front, NO blanks!) -!=============================================================================== ! -! This is only a dummy parameter needed for the regression check -#MPI= - -!=============================================================================== ! -! OUTPUT -!=============================================================================== ! - ProjectName = parallel_plates ! 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,,li,0.,0.0,,li,lx,0.0,,0.,lx,0.0 ,,0.,0.,lx,,li,0.,lx,,li,lx,lx,,0.,lx,lx/) ! [-3,3]x[-3,3]x[-3,3] - nElems =(/5,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.,lx/) ! 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.,lx,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 - -!=============================================================================== ! -! BASIS -!=============================================================================== ! - NVisu = 4 - -!=============================================================================== ! -! SEARCH -!=============================================================================== ! -! nElemsNodeSearch=50 -! RefineSideSearch=50 diff --git a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/parallel_plates_mesh.h5 b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/parallel_plates_mesh.h5 deleted file mode 100644 index c6de56e2c..000000000 Binary files a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/parallel_plates_mesh.h5 and /dev/null differ diff --git a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/parameter.ini b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/parameter.ini index 59c8e8215..38aea157d 100644 --- a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/parameter.ini +++ b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/parameter.ini @@ -10,7 +10,7 @@ GeometricNGeo = 1 ! Degree of mesh representation ! =============================================================================== ! ! MESH ! =============================================================================== ! -MeshFile = parallel_plates_mesh.h5 +MeshFile = parallel_plates_BCType2_DC_mesh.h5,parallel_plates_BCType60_AC_mesh.h5 useCurveds = F ! =============================================================================== ! ! OUTPUT / VISUALIZATION @@ -46,8 +46,13 @@ VisuParticles = T ! =============================================================================== ! ! Field Boundaries ! =============================================================================== ! -roundaryName = BC_right ! Adjust the right BC automatically depending on the actually input power -BoundaryType = (/2,2/) ! 2: Dirichlet with automatically adjusted electric potential +! This BC is set directly in the two hopr.ini files +!roundaryName = BC_right ! Adjust the right BC automatically depending on the actually input power +!BoundaryType = (/2,2/) ! 2: Dirichlet with automatically adjusted electric potential + +! The refstate for the AC case must be defined here +!RefState = (/10.0 , 3e9 , -1.57079632679/) ! RefState Nbr 1: Voltage, Frequency and Phase shift +RefState = (/1000.0 , 0.0 , 0.0/) ! RefState Nbr 1: Voltage, Frequency and Phase shift CoupledPowerPotential = (/10. , 1000. , 2000./) ! lower, starting and maximum values for the electric potential at all BoundaryType = (/2,2/) BCs CoupledPowerTarget = 1e-10 ! target power of 1e-10 Watt diff --git a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/pre-hopr-AC/hopr.ini b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/pre-hopr-AC/hopr.ini new file mode 100644 index 000000000..6e6fc67d9 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/pre-hopr-AC/hopr.ini @@ -0,0 +1,45 @@ +DEFVAR=(INT): iz = 1 ! no. of nFineHexa +DEFVAR=(REAL): li = 1 ! length +DEFVAR=(REAL): lx = 1 ! length +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! +ProjectName = parallel_plates_BCType60_AC ! 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,,li,0.,0.0,,li,lx,0.0,,0.,lx,0.0 ,,0.,0.,lx,,li,0.,lx,,li,lx,lx,,0.,lx,lx/) ! [-3,3]x[-3,3]x[-3,3] +nElems =(/5,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 + +!=============================================================================== ! +! 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.,lx/) ! 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.,lx,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=(/2,0,2,0/) ! ideal conductor diff --git a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/pre-hopr-DC/hopr.ini b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/pre-hopr-DC/hopr.ini new file mode 100644 index 000000000..6a5e468fc --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/pre-hopr-DC/hopr.ini @@ -0,0 +1,45 @@ +DEFVAR=(INT): iz = 1 ! no. of nFineHexa +DEFVAR=(REAL): li = 1 ! length +DEFVAR=(REAL): lx = 1 ! length +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! +ProjectName = parallel_plates_BCType2_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,,li,0.,0.0,,li,lx,0.0,,0.,lx,0.0 ,,0.,0.,lx,,li,0.,lx,,li,lx,lx,,0.,lx,lx/) ! [-3,3]x[-3,3]x[-3,3] +nElems =(/5,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 + +!=============================================================================== ! +! 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.,lx/) ! 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.,lx,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=(/60,0,1,0/) ! ideal conductor diff --git a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/readme.md b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/readme.md index e885e9661..25146562a 100644 --- a/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/readme.md +++ b/regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_fixed_power_input/readme.md @@ -9,6 +9,13 @@ CoupledPowerTarget = 1e-10 ! target power of 1e-10 Watt CoupledPowerRelaxFac = 0.5 ! Relaxation factor: Adjust the potential slowly to avoid oscillations +* the second test utilizes the AC boundary, but with f=0 (hence resulting in a DC potential) to test if the AC boundary fallback works + + BoundaryName = BC_right ! Adjust the right BC automatically depending on the actually input power + BoundaryType = (/60,1/) ! 1: RefState + + * adjusting an actual AC boundary with a single particle does not work! + * Coupled Power output and adjusted electric potential are compared with reference solution * parameter.ini for activating coupled Power output diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/analyze.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/analyze.ini index 205f85175..91b4f21f3 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/analyze.ini +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/analyze.ini @@ -1,9 +1,9 @@ ! 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 +!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 diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/externals.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/externals.ini index e9d673cb3..28e2ebafc 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/externals.ini +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/externals.ini @@ -1,8 +1,8 @@ ! --- Externals Tool Reggie -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../. +MPI = 1 , 1 , 1 +externalbinary = ./hopr/build/bin/hopr , ./hopr/build/bin/hopr , ./hopr/build/bin/hopr +externaldirectory = pre-hopr-DC , pre-hopr-AC , pre-hopr-AC-52 +externalruntime = pre , 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../. , ln\s-s\s./pre-hopr-AC-52/turner2013_BCType52_AC_mesh.h5\s../. nocrosscombination:MPI,externalbinary,externaldirectory,externalruntime,cmd_pre_execute diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/parameter.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/parameter.ini index 5a659123e..43ec0bdea 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/parameter.ini +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/parameter.ini @@ -10,7 +10,7 @@ NAnalyze = 10 ! Number of analyze points ! =============================================================================== ! ! MESH ! =============================================================================== ! -MeshFile = turner2013_BCType50_DC_mesh.h5, turner2013_BCType51_AC_mesh.h5 +MeshFile = turner2013_BCType50_DC_mesh.h5, turner2013_BCType51_AC_mesh.h5, turner2013_BCType52_AC_mesh.h5 useCurveds = F ! =============================================================================== ! ! OUTPUT / VISUALIZATION @@ -44,6 +44,8 @@ Particles-DSMC-CalcQualityFactors = F ! piclas: Pmax/Pmean CalcPointsPerDebyeLength = T CalcPICCFLCondition = T CalcMaxPartDisplacement = T + +CalcCoupledPower = T ! =============================================================================== ! ! CALCULATION ! =============================================================================== ! @@ -61,12 +63,22 @@ Analyze_dt = 5.0E-10 ! =============================================================================== ! ! This BC is defined in hopr.ini !BoundaryName = BC_left -!BoundaryType = (/50,0/),(/51,0/) ! Dirichlet with 0V initial BC, either DC or AC +!BoundaryType = (/50,0/),(/51,1/),(/51,2/) ! Dirichlet with 0V initial BC, either 1.) DC or 2.) AC or 3.) AC+power control ! 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 +! Power control for AC, only required for BoundaryType = (/52,1/) +CoupledPowerPotential = (/1.0 , 2.0 , 20.0/) ! lower, starting and maximum values for the electric potential at all BoundaryType = (/2,2/) BCs +CoupledPowerTarget = 10e-6 ! target power of 1e-10 Watt +CoupledPowerRelaxFac = 0.01 + +CoupledPowerMode = 1,2,3 +CoupledPowerFrequency = 3.333e9, 3.333e9,3.333e9 + +nocrosscombination:MeshFile,CoupledPowerMode,CoupledPowerFrequency + ! Boundary Field Output: Write current voltage to FieldAnalyze.csv CalcBoundaryFieldOutput = T BFO-NFieldBoundaries = 1 ! Nbr of boundaries diff --git a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/pre-hopr-AC-52/hopr.ini b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/pre-hopr-AC-52/hopr.ini new file mode 100644 index 000000000..9273693a2 --- /dev/null +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/pre-hopr-AC-52/hopr.ini @@ -0,0 +1,43 @@ +!=============================================================================== ! +! OUTPUT +!=============================================================================== ! +ProjectName = turner2013_BCType52_AC ! 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=(/52,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_AC-DC/readme.md b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/readme.md index 6a6300fb4..0713af470 100644 --- a/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/readme.md +++ b/regressioncheck/NIG_PIC_poisson_RK3/turner_bias-voltage_AC-DC/readme.md @@ -1,5 +1,8 @@ -# 1D turner setup and bias voltage feature +# 1D turner setup and bias voltage feature and power control - 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) +- bias voltage for + - DC: BCType=50 (starting from 0V) + - AC: BCType=51 (requiring a refstate to define the sinusoidal function, starting from 10.) + - AC+power control: BCType=52 which utilizes a bias voltage offset and cos function amplitude via power control (starting from 1V) - 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/analyze/analyzefield.f90 b/src/analyze/analyzefield.f90 index 497d31b62..3e9080bc2 100644 --- a/src/analyze/analyzefield.f90 +++ b/src/analyze/analyzefield.f90 @@ -1798,7 +1798,7 @@ END SUBROUTINE CalculateElectricDisplacementCurrentSurface #endif /*USE_HDG*/ !=================================================================================================================================== -!> Determine the field boundary condition values for the given iBC +!> Determine the field boundary output (BFO) values for the given iBC !=================================================================================================================================== SUBROUTINE CalculateBoundaryFieldOutput(iBC,Time,BoundaryFieldOutput) ! MODULES @@ -1834,12 +1834,8 @@ SUBROUTINE CalculateBoundaryFieldOutput(iBC,Time,BoundaryFieldOutput) CALL ExactFunc(BCState , x , BoundaryFieldOutput , t=Time) CASE(4) ! exact BC = Dirichlet BC !! Zero potential BoundaryFieldOutput = 0. - 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(5,51,52,60) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional) + 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 diff --git a/src/equations/poisson/equation.f90 b/src/equations/poisson/equation.f90 index 14fafc547..d04c7aba2 100644 --- a/src/equations/poisson/equation.f90 +++ b/src/equations/poisson/equation.f90 @@ -46,9 +46,6 @@ MODULE MOD_Equation PUBLIC :: InitEquation,ExactFunc,CalcSource,FinalizeEquation, CalcSourceHDG,DivCleaningDamping #if USE_MPI && defined(PARTICLES) -INTERFACE SynchronizeCPP - MODULE PROCEDURE SynchronizeCPP -END INTERFACE PUBLIC :: SynchronizeCPP #endif /*USE_MPI && defined(PARTICLES)*/ !=================================================================================================================================== @@ -62,6 +59,10 @@ SUBROUTINE DefineParametersEquation() ! MODULES USE MOD_Globals USE MOD_ReadInTools ,ONLY: prms +#if defined(PARTICLES) +USE MOD_HDG_Vars ,ONLY: CPPDataLength +USE MOD_ReadInTools ,ONLY: addStrListEntry +#endif /*defined(PARTICLES)*/ IMPLICIT NONE !================================================================================================================================== CALL prms%SetSection("Equation") @@ -88,9 +89,18 @@ SUBROUTINE DefineParametersEquation() #if defined(PARTICLES) ! 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%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=CPPDataLength ) 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 (which is 5%)', '0.05' ) +CALL prms%CreateRealOption( 'CoupledPowerFrequency' , 'Adaption frequency for coupled power using the integrated power for adjustment', '0.0' ) + +CALL prms%CreateIntFromStringOption('CoupledPowerMode', 'Define mode used for coupled power adjustment:\n'//& + ' instantaneous (1): xxx\n'//& + ' moving-average (2): xxx\n'//& + ' integrated (3): xxx\n','instantaneous') +CALL addStrListEntry('CoupledPowerMode' , 'instantaneous' , 1) +CALL addStrListEntry('CoupledPowerMode' , 'moving-average', 2) +CALL addStrListEntry('CoupledPowerMode' , 'integrated' , 3) #endif /*defined(PARTICLES)*/ END SUBROUTINE DefineParametersEquation @@ -128,7 +138,7 @@ SUBROUTINE InitEquation() CHARACTER(LEN=255) :: BCName INTEGER :: nRefStateMax INTEGER :: nLinState,nLinStateMax -INTEGER,PARAMETER :: BYTypeRefstate(1:2)=(/5,51/) +INTEGER,PARAMETER :: BYTypeRefstate(1:4)=(/5,51,52,60/) CHARACTER(LEN=32) :: hilf !=================================================================================================================================== IF((.NOT.InterpolationInitIsDone).OR.EquationInitIsDone)THEN @@ -275,8 +285,9 @@ SUBROUTINE InitCoupledPowerPotential() 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_HDG_Vars ,ONLY: UseCoupledPowerPotential,CoupledPowerPotential,CoupledPowerTarget,CoupledPowerRelaxFac +USE MOD_HDG_Vars ,ONLY: CPPDataLength,CoupledPowerMode,CoupledPowerFrequency +USE MOD_ReadInTools ,ONLY: GETREALARRAY,GETREAL,GETINTFROMSTR,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 @@ -289,9 +300,10 @@ SUBROUTINE InitCoupledPowerPotential() ! INPUT / OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! 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, PARAMETER :: BCTypeCPP(1:3) = (/2,52,60/) ! BCType which allows coupled power potential control +! ! 2: DC + coupled power for DC potential adjustment +! ! 52: AC with cos(wt) + coupled power for AC potential adjustment + bias voltage +! ! 60: AC with cos(wt) + coupled power for AC potential adjustment INTEGER :: iBC,CPPBoundaries,BCType,BCState #if USE_MPI LOGICAL :: BConProc @@ -300,7 +312,7 @@ SUBROUTINE InitCoupledPowerPotential() !=================================================================================================================================== ! 1.) Get global number of coupled power potential boundaries in [1:nBCs] -CalcPCouplElectricPotential = .FALSE. +UseCoupledPowerPotential = .FALSE. CPPBoundaries = 0 ! Loop over global BC list DO iBC=1,nBCs @@ -314,7 +326,7 @@ SUBROUTINE InitCoupledPowerPotential() IF(CPPBoundaries.EQ.0) RETURN ! Already determined in HDG initialization -CalcPCouplElectricPotential = .TRUE. +UseCoupledPowerPotential = .TRUE. ! 2.) Get parameters #if USE_LOADBALANCE @@ -322,13 +334,24 @@ SUBROUTINE InitCoupledPowerPotential() 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) + CoupledPowerPotential = 0. + CoupledPowerPotential(1:3) = GETREALARRAY('CoupledPowerPotential',3) ! Get target power value CoupledPowerTarget = GETREAL('CoupledPowerTarget') IF(CoupledPowerTarget.LE.0.) CALL CollectiveStop(__STAMP__,'CoupledPowerTarget must be > 0.') ! Get relaxation factor CoupledPowerRelaxFac = GETREAL('CoupledPowerRelaxFac') IF(CoupledPowerRelaxFac.LE.0.) CALL CollectiveStop(__STAMP__,'CoupledPowerRelaxFac must be > 0.') + ! Get frequency within which the integrated power is calculated + CoupledPowerFrequency = GETREAL('CoupledPowerFrequency') + IF(CoupledPowerFrequency.LT.0.) CALL CollectiveStop(__STAMP__,'CoupledPowerFrequency must be >= 0.') + ! Sanity check + CoupledPowerMode = GETINTFROMSTR('CoupledPowerMode') + IF(.NOT.ANY(CoupledPowerMode.EQ.(/1,2,3/))) CALL CollectiveStop(__STAMP__,'CoupledPowerMode with unknown value!') + IF((CoupledPowerMode.EQ.3).AND.(ABS(CoupledPowerFrequency).LE.0.))& + CALL CollectiveStop(__STAMP__,'CoupledPowerMode=3 requires a positive value for CoupledPowerFrequency!') + ! Update time + IF(CoupledPowerFrequency.GT.0.0) CoupledPowerPotential(5) = 1.0 / CoupledPowerFrequency #if USE_LOADBALANCE END IF ! .NOT.PerformLoadBalance #endif /*USE_LOADBALANCE*/ @@ -339,7 +362,7 @@ SUBROUTINE InitCoupledPowerPotential() IF(MPIRoot)THEN BConProc = .TRUE. ELSE - CoupledPowerPotential(1:3) = 0. + CoupledPowerPotential(1:CPPDataLength) = 0. ! Check local sides DO SideID=1,nBCSides iBC = BC(SideID) @@ -397,7 +420,7 @@ SUBROUTINE ReadCPPDataFromH5() 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 +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,CPPDataLength IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------! ! INPUT / OUTPUT VARIABLES @@ -405,7 +428,7 @@ SUBROUTINE ReadCPPDataFromH5() ! LOCAL VARIABLES CHARACTER(255) :: ContainerName LOGICAL :: CPPExists -REAL :: CPPDataHDF5(1:3) +REAL :: CPPDataHDF5(1:CPPDataLength) !=================================================================================================================================== ! Only required during restart IF(.NOT.DoRestart) RETURN @@ -424,9 +447,14 @@ SUBROUTINE ReadCPPDataFromH5() 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) + CALL ReadArray(TRIM(ContainerName) , 2 , (/1_IK , INT(CPPDataLength,IK)/) , 0_IK , 1 , RealArray=CPPDataHDF5) + WRITE(UNIT_stdOut,'(6(A,ES10.2E3))') " Read CoupledPowerPotential from restart file ["//TRIM(RestartFile)//& + "] min[V]: " ,CPPDataHDF5(1),& + ", current[V]: " ,CPPDataHDF5(2),& + ", max[V]: " ,CPPDataHDF5(3),& + ", intPower[Ws]: " ,CPPDataHDF5(4),& + ", next adjustment time[s]: ",CPPDataHDF5(5),& + ", last-energy[J]: " ,CPPDataHDF5(6) CoupledPowerPotential = CPPDataHDF5 END IF ! CPPExists CALL CloseDataFile() @@ -448,7 +476,7 @@ 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 +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,CPPDataLength IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------! ! INPUT / OUTPUT VARIABLES @@ -457,7 +485,7 @@ SUBROUTINE SynchronizeCPP() !=================================================================================================================================== 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) + CALL MPI_BCAST(CoupledPowerPotential, CPPDataLength, MPI_DOUBLE_PRECISION, 0, CPPCOMM%UNICATOR, IERROR) END IF END SUBROUTINE SynchronizeCPP #endif /*USE_MPI*/ @@ -473,7 +501,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,BiasVoltage +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,UseCoupledPowerPotential,BiasVoltage,UseBiasVoltage USE MOD_Particle_Vars ,ONLY: Species,nSpecies #endif /*defined(PARTICLES)*/ USE MOD_Dielectric_Vars ,ONLY: DielectricRatio,Dielectric_E_0,DielectricRadiusValue,DielectricEpsR @@ -528,15 +556,22 @@ 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,51) ! Signal with zero-crossing: Amplitude, Frequency and Phase Shift supplied by RefState +CASE(-1) ! 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)) #if defined(PARTICLES) - ! Add bias potential (only if bias voltage model is activated) - IF(ExactFunction.EQ.51) Resu(:) = Resu(:) + BiasVoltage%BVData(1) + ! Match coupled power by adjusting the coefficient of the cos function + IF(UseCoupledPowerPotential)THEN + Resu(:) = CoupledPowerPotential(2)*COS(Omega*t+RefState(3,iRefState)) + ELSE +#endif /*defined(PARTICLES)*/ + Resu(:) = RefState(1,iRefState)*COS(Omega*t+RefState(3,iRefState)) +#if defined(PARTICLES) + END IF ! UseCoupledPowerPotential + ! Add bias potential (only if bias voltage model is activated, BYType is 51 for DC or 52 for AC) + IF(UseBiasVoltage) Resu(:) = Resu(:) + BiasVoltage%BVData(1) #endif /*defined(PARTICLES)*/ CASE(0) ! constant 0. Resu(:)=0. diff --git a/src/hdg/hdg.f90 b/src/hdg/hdg.f90 index 1f4a4b10c..394f078c0 100644 --- a/src/hdg/hdg.f90 +++ b/src/hdg/hdg.f90 @@ -1933,7 +1933,7 @@ SUBROUTINE HDGLinear(time,U_out) r=q*(PP_N+1) + p+1 lambda(iVar,r:r,SideID)=0. END DO; END DO !p,q - CASE(5) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional) for reference state (with zero crossing) + CASE(5,51,52,60) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional) for reference state (with zero crossing) DO q=0,PP_N; DO p=0,PP_N r=q*(PP_N+1) + p+1 CALL ExactFunc(-1,Face_xGP(:,p,q,SideID),lambda(iVar,r:r,SideID),t=time,iRefState=BCState) @@ -1958,11 +1958,6 @@ 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) @@ -2403,8 +2398,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,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!') + CASE(8,50,51,52,60) ! exact BC = Dirichlet BC !! ExactFunc via electric potential and decharing of a surface + CALL abort(__STAMP__,'Dirichlet BC=8,50,51,52,60 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 74e98b1be..1c1d83384 100644 --- a/src/hdg/hdg_vars.f90 +++ b/src/hdg/hdg_vars.f90 @@ -238,15 +238,18 @@ 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 :: UseCoupledPowerPotential !< Switch calculation on/off +INTEGER,PARAMETER :: CPPDataLength=6 !< Number of variables in BVData #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 +REAL :: CoupledPowerPotential(CPPDataLength) !< (/min, start, max/) electric potential, e.g., used 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 :: CoupledPowerFrequency !< Frequency with which the integrated power is calculated (must be consistent Part-AnalyzeStep, i.e., that one cycle with period T=1/f must be larger than Part-AnalyzeStep * dt) +INTEGER :: CoupledPowerMode !< Method for power adjustment with 1: instantaneous power, 2: moving average power, 3: integrated power !=================================================================================================================================== !-- Bias Voltage (for calculating a BC voltage bias for certain BCs) diff --git a/src/io_hdf5/hdf5_output_state.f90 b/src/io_hdf5/hdf5_output_state.f90 index 102393f2b..e63dd87e2 100644 --- a/src/io_hdf5/hdf5_output_state.f90 +++ b/src/io_hdf5/hdf5_output_state.f90 @@ -102,7 +102,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_HDG_Vars ,ONLY: CoupledPowerPotential,CalcPCouplElectricPotential +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,UseCoupledPowerPotential,CPPDataLength #endif /*PARTICLES*/ #endif /*USE_HDG*/ USE MOD_Analyze_Vars ,ONLY: OutputTimeFixed @@ -134,7 +134,7 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) REAL :: NumSpec(nSpecAnalyze),TmpArray(1,1) INTEGER(KIND=IK) :: SimNumSpec(nSpecAnalyze) #if USE_HDG -REAL :: TmpArray2(3,1) +REAL,ALLOCATABLE :: CPPDataHDF5(:,:) #endif /*USE_HDG*/ #endif /*PARTICLES*/ REAL :: StartT,EndT @@ -578,19 +578,21 @@ SUBROUTINE WriteStateToHDF5(MeshFileName,OutputTime,PreviousTime) CALL CloseDataFile() END IF ! CalcBulkElectronTempi.AND.MPIRoot #if USE_HDG -IF(CalcPCouplElectricPotential.AND.MPIRoot)THEN +IF(UseCoupledPowerPotential.AND.MPIRoot)THEN + ALLOCATE(CPPDataHDF5(1:CPPDataLength,1)) CALL OpenDataFile(FileName,create=.FALSE.,single=.TRUE.,readOnly=.FALSE.) - TmpArray2(1:3,1) = CoupledPowerPotential(1:3) + CPPDataHDF5(1:CPPDataLength,1) = CoupledPowerPotential(1:CPPDataLength) CALL WriteArrayToHDF5( DataSetName = 'CoupledPowerPotential' , rank = 2 , & - nValGlobal = (/1_IK , 3_IK/) , & - nVal = (/1_IK , 3_IK/) , & + nValGlobal = (/1_IK , INT(CPPDataLength,IK)/) , & + nVal = (/1_IK , INT(CPPDataLength,IK)/) , & offset = (/0_IK , 0_IK/) , & - collective = .FALSE., RealArray = TmpArray2(1:3,1)) + collective = .FALSE., RealArray = CPPDataHDF5(1:CPPDataLength,1)) CALL CloseDataFile() + DEALLOCATE(CPPDataHDF5) END IF ! CalcBulkElectronTempi.AND.MPIRoot ! Bias voltage IF(UseBiasVoltage.AND.MPIRoot)THEN - ALLOCATE(BVDataHDF5(BVDataLength,1)) + ALLOCATE(BVDataHDF5(1:BVDataLength,1)) CALL OpenDataFile(FileName,create=.FALSE.,single=.TRUE.,readOnly=.FALSE.) BVDataHDF5(1:BVDataLength,1) = BiasVoltage%BVData(1:BVDataLength) CALL WriteArrayToHDF5( DataSetName = 'BiasVoltage' , rank = 2 , & diff --git a/src/particles/analyze/particle_analyze.f90 b/src/particles/analyze/particle_analyze.f90 index 97a309536..70ada040b 100644 --- a/src/particles/analyze/particle_analyze.f90 +++ b/src/particles/analyze/particle_analyze.f90 @@ -138,12 +138,13 @@ SUBROUTINE InitParticleAnalyze() USE MOD_Particle_Analyze_Tools,ONLY: AllocateElectronIonDensityCell,AllocateElectronTemperatureCell,AllocateCalcElectronEnergy #if USE_HDG USE MOD_HDG_Vars ,ONLY: CalcBRVariableElectronTemp,BRAutomaticElectronRef -USE MOD_HDG_Vars ,ONLY: CalcPCouplElectricPotential +USE MOD_HDG_Vars ,ONLY: UseCoupledPowerPotential #endif /*USE_HDG*/ USE MOD_Particle_Analyze_Tools,ONLY: CalcNumberDensityBGGasDistri #if USE_LOADBALANCE USE MOD_LoadBalance_Vars ,ONLY: PerformLoadBalance #endif /*USE_LOADBALANCE*/ +USE MOD_Restart_Vars ,ONLY: DoRestart,RestartTime ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- @@ -167,6 +168,13 @@ SUBROUTINE InitParticleAnalyze() LBWRITE(UNIT_StdOut,'(132("-"))') LBWRITE(UNIT_stdOut,'(A)') ' INIT PARTICLE ANALYZE...' +! Initialize with restart time or zero +IF(DoRestart)THEN + ParticleAnalyzeSampleTime = RestartTime +ELSE + ParticleAnalyzeSampleTime = 0. +END IF ! DoRestart + ! Average number of points per shape function: max. number allowed is (PP_N+1)^3 IF(StringBeginsWith(DepositionType,'shape_function'))THEN CalcPointsPerShapeFunction = GETLOGICAL('CalcPointsPerShapeFunction') @@ -563,13 +571,14 @@ SUBROUTINE InitParticleAnalyze() !-- Coupled Power CalcCoupledPower = GETLOGICAL('CalcCoupledPower') #if USE_HDG -IF(CalcPCouplElectricPotential.AND.(.NOT.CalcCoupledPower)) CALL abort(__STAMP__,'BoundaryType = (/2,2/) requires CalcCoupledPower=T') +IF(UseCoupledPowerPotential.AND.(.NOT.CalcCoupledPower)) CALL abort(__STAMP__,'Coupled power potential requires CalcCoupledPower=T') #endif /*USE_HDG*/ IF(CalcCoupledPower) THEN DisplayCoupledPower = GETLOGICAL('DisplayCoupledPower') DoPartAnalyze = .TRUE. PCouplAverage = 0.0 + PCouplAverageOld = 0.0 #if !((PP_TimeDiscMethod==1) || (PP_TimeDiscMethod==2) || (PP_TimeDiscMethod==6) || (PP_TimeDiscMethod==500) || (PP_TimeDiscMethod==501) || (PP_TimeDiscMethod==502) || (PP_TimeDiscMethod==506) || (PP_TimeDiscMethod==507) || (PP_TimeDiscMethod==508) || (PP_TimeDiscMethod==509)) CALL abort(__STAMP__,'ERROR: CalcCoupledPower is not implemented yet with the chosen time discretization method!') #endif @@ -852,7 +861,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_HDG_Vars ,ONLY: CalcPCouplElectricPotential,CoupledPowerPotential +USE MOD_HDG_Vars ,ONLY: UseCoupledPowerPotential,CoupledPowerPotential,CoupledPowerFrequency,CoupledPowerMode USE MOD_Particle_Analyze_Tools ,ONLY: CalculatePCouplElectricPotential #endif /*USE_HDG*/ USE MOD_Globals_Vars ,ONLY: eV2Kelvin @@ -893,11 +902,12 @@ SUBROUTINE AnalyzeParticles(Time) #if USE_MPI REAL :: tmpArray(1:2) #endif /*USE_MPI*/ +REAL :: PCouplDelta +REAL :: TimeDelta !=================================================================================================================================== - IF ( DoRestart ) THEN - isRestart = .true. - END IF - IF (.NOT.DoPartAnalyze) RETURN +IF(DoRestart) isRestart = .true. +IF(.NOT.DoPartAnalyze) RETURN +ParticleAnalyzeSampleTime = Time - ParticleAnalyzeSampleTime ! Set ParticleAnalyzeSampleTime=Time at the end of this routine #if (PP_TimeDiscMethod==42) IF (useDSMC) THEN IF (CollisMode.NE.0) THEN @@ -1003,6 +1013,9 @@ SUBROUTINE AnalyzeParticles(Time) WRITE(unit_index,'(A1)',ADVANCE='NO') ',' WRITE(unit_index,'(I3.3,A)',ADVANCE='NO') OutputCounter,'-PCoupledMoAv' OutputCounter = OutputCounter + 1 + WRITE(unit_index,'(A1)',ADVANCE='NO') ',' + WRITE(unit_index,'(I3.3,A)',ADVANCE='NO') OutputCounter,'-PCoupledIntAv' + OutputCounter = OutputCounter + 1 END IF IF (CalcLaserInteraction) THEN ! computer laser-plasma interaction DO iSpec=1, nSpecies @@ -1253,9 +1266,9 @@ SUBROUTINE AnalyzeParticles(Time) OutputCounter = OutputCounter + 1 END DO END IF ! CalcBRVariableElectronTemp.OR.BRAutomaticElectronRef - IF(CalcPCouplElectricPotential)THEN + IF(UseCoupledPowerPotential)THEN WRITE(unit_index,'(A1,I3.3,A)',ADVANCE='NO') ',',OutputCounter,'-CoupledPowerPotential-[V]' - END IF ! CalcPCouplElectricPotential + END IF ! UseCoupledPowerPotential #endif /*USE_HDG*/ IF(CalcBulkElectronTemp)THEN WRITE(unit_index,'(A1,I3.3,A,I3.3,A)',ADVANCE='NO') ',',OutputCounter,'-BulkElectronTemp-[K]' @@ -1353,14 +1366,40 @@ SUBROUTINE AnalyzeParticles(Time) ! MPI Communication for values which are not YET communicated ! All routines ABOVE contain the required MPI-Communication !=================================================================================================================================== -#if USE_MPI IF(CalcCoupledPower) THEN + PCouplIntAverage = 0.0 ! Default +#if USE_MPI + ! Collect sum on MPIRoot tmpArray = (/PCoupl, PCouplAverage/) - CALL MPI_ALLREDUCE(MPI_IN_PLACE , tmpArray , 2 , MPI_DOUBLE_PRECISION , MPI_SUM , MPI_COMM_WORLD , IERROR) - PCoupl = tmpArray(1) - PCouplAverage = tmpArray(2) - ! Only the root process keeps PCouplAverage, but all processes require PCoupl if CalcPCouplElectricPotential=T - IF(.NOT.MPIRoot) PCouplAverage = 0. + IF(PartMPI%MPIRoot)THEN + CALL MPI_REDUCE(MPI_IN_PLACE, tmpArray, 2, MPI_DOUBLE_PRECISION, MPI_SUM, 0, PartMPI%COMM, IERROR) + PCoupl = tmpArray(1) + PCouplAverage = tmpArray(2) +#endif /*USE_MPI*/ +#if USE_HDG + ! Only required for integrated average over one cycle + IF(CoupledPowerMode.EQ.3)THEN + PCouplDelta = PCouplAverage - PCouplAverageOld ! y2 - new value + PCouplAverageOld = PCouplAverage + CoupledPowerPotential(4) = CoupledPowerPotential(4) + 0.5*(CoupledPowerPotential(6) + PCouplDelta) ! 0.5*dx*(y1+y2)/dx + CoupledPowerPotential(6) = PCouplDelta ! y1 = y2 - store old value ("last energy") + ! Check sampling frequency + IF(CoupledPowerFrequency.GT.0)THEN + TimeDelta = (1.0 / CoupledPowerFrequency) + Time - CoupledPowerPotential(5) ! = 1/f + t - tCPP + ELSE + TimeDelta = ParticleAnalyzeSampleTime + END IF ! CoupledPowerFrequency.GT.0 + ! Check sampling time + IF(ABS(TimeDelta).GT.0.) PCouplIntAverage = CoupledPowerPotential(4) / TimeDelta ! Calculate average from integral + END IF ! CoupledPowerMode.EQ.3 +#endif /*USE_HDG*/ +#if USE_MPI + ELSE + CALL MPI_REDUCE(tmpArray, 0, 2, MPI_DOUBLE_PRECISION, MPI_SUM, 0, PartMPI%COMM, IERROR) + ! Reset for all processes execpt the MPIRoot (this process keeps the old value, which is therefore considered only once in the + ! next MPI reduce call) + PCouplAverage = 0. + END IF ! PartMPI%MPIRoot END IF ! Switch between root and non-root processes IF (PartMPI%MPIRoot) THEN @@ -1377,20 +1416,21 @@ SUBROUTINE AnalyzeParticles(Time) CALL MPI_REDUCE(PartEkinIn ,0,nSpecAnalyze,MPI_DOUBLE_PRECISION,MPI_SUM,0,PartMPI%COMM,IERROR) CALL MPI_REDUCE(PartEkinOut,0,nSpecAnalyze,MPI_DOUBLE_PRECISION,MPI_SUM,0,PartMPI%COMM,IERROR) END IF - END IF #endif /*USE_MPI*/ + END IF !----------------------------------------------------------------------------------------------------------------------------------- ! Analyze Routines that require MPI_REDUCE of other variables ! Moving Average of PCoupl: -IF(CalcCoupledPower) THEN +IF(CalcCoupledPower.AND.PartMPI%MPIRoot) THEN ! Moving Average of PCoupl: - IF(ABS(Time-RestartTime).GT.0.0) PCouplAverage = PCouplAverage / (Time-RestartTime) + TimeDelta = Time-RestartTime + IF(ABS(TimeDelta).GT.0.0) PCouplAverage = PCouplAverage / TimeDelta ! current PCoupl (Delta_E / Timestep) PCoupl = PCoupl / dt END IF #if USE_HDG ! Calculate electric potential for special BCs BoundaryType = (/2,2/) to meet a specific input power -IF((iter.GT.0).AND.CalcPCouplElectricPotential) CALL CalculatePCouplElectricPotential() +IF((iter.GT.0).AND.UseCoupledPowerPotential) CALL CalculatePCouplElectricPotential() #endif /*USE_HDG*/ !----------------------------------------------------------------------------------------------------------------------------------- ! Perform averaging/summation of the MPI communicated variables on the root only (and for the non-MPI case, MPIRoot is set to true) @@ -1488,6 +1528,7 @@ SUBROUTINE AnalyzeParticles(Time) IF (CalcCoupledPower) THEN WRITE(unit_index,CSVFORMAT,ADVANCE='NO') ',', PCoupl WRITE(unit_index,CSVFORMAT,ADVANCE='NO') ',', PCouplAverage + WRITE(unit_index,CSVFORMAT,ADVANCE='NO') ',', PCouplIntAverage END IF IF (CalcLaserInteraction) THEN DO iSpec=1, nSpecies @@ -1637,9 +1678,9 @@ SUBROUTINE AnalyzeParticles(Time) WRITE(unit_index,CSVFORMAT,ADVANCE='NO') ',', RegionElectronRef(3,iRegions)*ElementaryCharge/BoltzmannConst ! convert eV to K END DO END IF ! CalcBRVariableElectronTemp.OR.BRAutomaticElectronRef - IF(CalcPCouplElectricPotential)THEN + IF(UseCoupledPowerPotential)THEN WRITE(unit_index,CSVFORMAT,ADVANCE='NO') ',', CoupledPowerPotential(2) ! Electric potential in V - END IF ! CalcPCouplElectricPotential + END IF ! UseCoupledPowerPotential #endif /*USE_HDG*/ IF(CalcBulkElectronTemp)THEN WRITE(unit_index,CSVFORMAT,ADVANCE='NO') ',', BulkElectronTemp*eV2Kelvin ! Temperature in Kelvin @@ -1651,8 +1692,8 @@ SUBROUTINE AnalyzeParticles(Time) #endif /*USE_MPI*/ !----------------------------------------------------------------------------------------------------------------------------------- ! Reset coupled power to particles if output of coupled power is active -IF (CalcCoupledPower) THEN - PCouplAverage = PCouplAverage * (Time-RestartTime) ! PCouplAverage is reseted +IF (CalcCoupledPower.AND.PartMPI%MPIRoot) THEN + IF(ABS(TimeDelta).GT.0.0) PCouplAverage = PCouplAverage * TimeDelta ! PCouplAverage is reset END IF !----------------------------------------------------------------------------------------------------------------------------------- ! Reset the particle counter @@ -1661,6 +1702,7 @@ SUBROUTINE AnalyzeParticles(Time) END IF !----------------------------------------------------------------------------------------------------------------------------------- +ParticleAnalyzeSampleTime = Time ! Backup "old" time value for next output END SUBROUTINE AnalyzeParticles diff --git a/src/particles/analyze/particle_analyze_tools.f90 b/src/particles/analyze/particle_analyze_tools.f90 index 62bc3b295..30f8faf25 100644 --- a/src/particles/analyze/particle_analyze_tools.f90 +++ b/src/particles/analyze/particle_analyze_tools.f90 @@ -3093,25 +3093,64 @@ END SUBROUTINE CalcCoupledPowerPart !=================================================================================================================================== SUBROUTINE CalculatePCouplElectricPotential() ! MODULES -USE MOD_Globals -USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,CoupledPowerTarget,CoupledPowerRelaxFac -USE MOD_Particle_Analyze_Vars ,ONLY: PCoupl +USE MOD_Globals ,ONLY: MPIRoot +USE MOD_HDG_Vars ,ONLY: CoupledPowerPotential,CoupledPowerTarget,CoupledPowerRelaxFac,CoupledPowerFrequency +USE MOD_HDG_Vars ,ONLY: CoupledPowerMode +USE MOD_Particle_Analyze_Vars ,ONLY: PCoupl,PCouplAverage,PCouplIntAverage +USE MOD_Restart_Vars ,ONLY: RestartTime +USE MOD_TimeDisc_Vars ,ONLY: Time +#if USE_MPI +USE MOD_Equation ,ONLY: SynchronizeCPP +#endif /*USE_MPI*/ IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------! ! INPUT / OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! LOCAL VARIABLES -REAL :: PowerRatio +REAL :: PowerRatio !=================================================================================================================================== -! Adjust electric potential depending on the instantaneous coupled power -PowerRatio = CoupledPowerTarget / PCoupl +IF(MPIRoot)THEN + ASSOCIATE(& + Vmin => CoupledPowerPotential(1),& + V => CoupledPowerPotential(2),& + Vmax => CoupledPowerPotential(3),& + VInt => CoupledPowerPotential(4),& + tCPP => CoupledPowerPotential(5) & + ) + + PowerRatio = 1.0 ! Default + SELECT CASE(CoupledPowerMode) + CASE(1) + ! Adjust electric potential depending on the instantaneous coupled power + PowerRatio = CoupledPowerTarget / PCoupl + CASE(2) + ! Use moving average power + IF(ABS(Time-RestartTime).GT.0.0) PowerRatio = CoupledPowerTarget / PCouplAverage + CASE(3) + ! Use integrated power (via user-defined frequency) + IF(time.GE.tCPP)THEN ! Simulation time threshold reached + ! Update time + IF(CoupledPowerFrequency.GT.0.0) tCPP = tCPP + 1.0 / CoupledPowerFrequency + ! Update Voltage + PowerRatio = CoupledPowerTarget / PCouplIntAverage ! PCouplIntAverage is the running integrated average + ! Reset integrated power value + VInt = 0. + END IF ! time.GE.tBV + END SELECT -! Relaxation factor -CoupledPowerPotential(2) = CoupledPowerPotential(2) * (1.0 + CoupledPowerRelaxFac * (PowerRatio - 1.0)) + ! Relaxation factor + V = V * (1.0 + CoupledPowerRelaxFac * (PowerRatio - 1.0)) -! Keep boundaries -IF(CoupledPowerPotential(2).GT.CoupledPowerPotential(3)) CoupledPowerPotential(2) = CoupledPowerPotential(3) -IF(CoupledPowerPotential(2).LT.CoupledPowerPotential(1)) CoupledPowerPotential(2) = CoupledPowerPotential(1) + ! Keep boundaries + IF(V.GT.Vmax) V = Vmax + IF(V.LT.Vmin) V = Vmin + END ASSOCIATE +END IF ! MPIRoot + +#if USE_MPI +! Synchronize CoupledPowerPotential from MPIRoot to all processors on the sub-communicator +CALL SynchronizeCPP() +#endif /*USE_MPI*/ END SUBROUTINE CalculatePCouplElectricPotential #endif /*USE_HDG*/ diff --git a/src/particles/analyze/particle_analyze_vars.f90 b/src/particles/analyze/particle_analyze_vars.f90 index 2a165805a..9373a004b 100644 --- a/src/particles/analyze/particle_analyze_vars.f90 +++ b/src/particles/analyze/particle_analyze_vars.f90 @@ -23,6 +23,7 @@ MODULE MOD_Particle_Analyze_Vars ! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- LOGICAL :: ParticleAnalyzeInitIsDone = .FALSE. +REAL :: ParticleAnalyzeSampleTime !< Accumulated simulation time between two outputs to ParticleAnalyze.csv LOGICAL :: CalcSimNumSpec !< Calculate the number of simulated particles per species LOGICAL :: CalcNumDens !< Calculate the number density per species within the domain LOGICAL :: CalcSurfFluxInfo !< Calculate the current/mass flow through or pressure (adaptive/subsonic BC) at the surface flux boundaries @@ -39,6 +40,8 @@ MODULE MOD_Particle_Analyze_Vars !< push (only charged particles) REAL :: PCoupl !< Power that is coupled into plasma in [W] REAL :: PCouplAverage !< Power that is coupled into plasma (moving average) in [W] +REAL :: PCouplIntAverage !< Power that is coupled into plasma (moving integrated average) in [W] +REAL :: PCouplAverageOld !< Power that is coupled into plasma (moving integrated average) in [W] - old value from last call to PartAnalyze() TYPE tPCoupl REAL,ALLOCATABLE :: DensityAvgElem(:) !< Power per volume that is coupled into plasma (moving average !< for each element) in [W/m^3] diff --git a/src/particles/particle_init.f90 b/src/particles/particle_init.f90 index ed1a162be..492959dc3 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_HDG_Vars ,ONLY: CalcPCouplElectricPotential +USE MOD_HDG_Vars ,ONLY: UseCoupledPowerPotential #endif /*defined(PARTICLES) && USE_HDG*/ ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE @@ -463,7 +463,7 @@ SUBROUTINE InitializeVariables() ! init DSMC determines whether DSMC%UseOctree is true or false) DoInterpolation = GETLOGICAL('PIC-DoInterpolation') #if defined(PARTICLES) && USE_HDG -IF(CalcPCouplElectricPotential.AND.(.NOT.DoInterpolation)) CALL abort(__STAMP__,'BoundaryType = (/2,2/) requires DoInterpolation=T') +IF(UseCoupledPowerPotential.AND.(.NOT.DoInterpolation)) CALL abort(__STAMP__,'Coupled power potential requires DoInterpolation=T') #endif /*defined(PARTICLES) && USE_HDG*/ #ifdef CODE_ANALYZE ! Check if an analytic function is to be used for interpolation diff --git a/src/particles/surfacemodel/surfacemodel_analyze.f90 b/src/particles/surfacemodel/surfacemodel_analyze.f90 index 5e5f15276..24ab955c0 100644 --- a/src/particles/surfacemodel/surfacemodel_analyze.f90 +++ b/src/particles/surfacemodel/surfacemodel_analyze.f90 @@ -158,7 +158,7 @@ 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_Vars ,ONLY: UseBiasVoltage,BiasVoltage,BVDataLength #if USE_MPI USE MOD_HDG ,ONLY: SynchronizeBV #endif /*USE_MPI*/ @@ -403,7 +403,7 @@ SUBROUTINE AnalyzeSurface(Time) END ASSOCIATE ! Write: Voltage, Ion excess and simulation update time - DO i = 1, 3 + DO i = 1, BVDataLength CALL WriteDataInfo(unit_index, RealScalar=BiasVoltage%BVData(i)) END DO ! i = 1, 3 END IF ! UseBiasVoltage diff --git a/src/piclas.h b/src/piclas.h index 4c6a7ed94..05082b39c 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,51 +#define HDGDIRICHLETBCSIDEIDS 2,4,5,6,7,8,50,51,52,60 #endif diff --git a/src/restart/restart_field.f90 b/src/restart/restart_field.f90 index 67c3c5dab..9a803e109 100644 --- a/src/restart/restart_field.f90 +++ b/src/restart/restart_field.f90 @@ -79,7 +79,7 @@ SUBROUTINE FieldRestart() #if defined(PARTICLES) USE MOD_Equation ,ONLY: SynchronizeCPP USE MOD_HDG ,ONLY: SynchronizeBV -USE MOD_HDG_Vars ,ONLY: UseBiasVoltage,CalcPCouplElectricPotential +USE MOD_HDG_Vars ,ONLY: UseBiasVoltage,UseCoupledPowerPotential ! 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 @@ -175,7 +175,7 @@ SUBROUTINE FieldRestart() ! 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(UseCoupledPowerPotential) CALL SynchronizeCPP() !checkRank=MIN(3,nProcessors-1) ! Store lambda solution on global non-unique array for MPI communication