Merge branch 'improvement.hdg.bc.dirichlet.linear.function' into 'mas…

…ter.dev'

[improvement.hdg.bc.dirichlet.linear.function] Adjusted the HDG linear potential to have the same functionality as RefStates

See merge request piclas/piclas!703

CMakeListsLib.txt

@@ -295,7 +295,7 @@ IF(NOT LIBS_BUILD_MATH_LIB)
   ENDIF()
 
   # Use Lapack/Blas for GNU
-  FIND_PACKAGE(LAPACK QUIET)
+  FIND_PACKAGE(LAPACK REQUIRED)
   IF (LAPACK_FOUND)
     LIST(APPEND linkedlibs ${LAPACK_LIBRARIES})
     MESSAGE(STATUS "Compiling with system [BLAS/Lapack]")

regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates/parameter.ini

@@ -45,20 +45,19 @@ VisuParticles        = T
 ! =============================================================================== !
 ! Field Boundaries
 ! =============================================================================== !
-BoundaryName = BC_right
-!BoundaryType = (/5,1/)              ! 5: Dirichlet, 1: Nbr of RefState
-!RefState     = (/1000.0, 0.0, 0.0/) ! RefState Nbr 1: Voltage and Frequency
-BoundaryType = (/2,1/)              ! 2: Dirichlet with linear ramp
-
-
-BoundaryName = BC_left
-!BoundaryType = (/4,0/) ! 4: Dirichlet with zero potential
-BoundaryType = (/2,1/) ! 2: Dirichlet with linear ramp
-
-LinPhiBasePoint = (/0. , 0. , 0./)
-LinPhiNormal    = (/1. , 0. , 0./)
-LinPhiHeight    = 1.0
-LinPhi          = 1e3
+BoundaryName    = BC_right
+BoundaryType    = (/7,1/)          ! 7: Dirichlet with linear ramp 1st LinState
+LinPhiBasePoint = (/0. , 0. , 0./) ! 1st LinState
+LinPhiNormal    = (/1. , 0. , 0./) ! 1st LinState
+LinPhiHeight    = 1.0              ! 1st LinState
+LinPhi          = 1e3              ! 1st LinState
+
+BoundaryName    = BC_left
+BoundaryType    = (/7,2/)          ! 7: Dirichlet with linear ramp 2nd LinState
+LinPhiBasePoint = (/0. , 0. , 0./) ! 2nd LinState
+LinPhiNormal    = (/1. , 0. , 0./) ! 2nd LinState
+LinPhiHeight    = 1.0              ! 2nd LinState
+LinPhi          = 0.0              ! 2nd LinState
 ! =============================================================================== !
 ! PARTICLES
 ! =============================================================================== !

regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_AC/analyze.ini

@@ -2,3 +2,10 @@ compare_data_file_name      = PartAnalyze.csv
 compare_data_file_reference = PartAnalyzeLeapfrog_ref.csv
 compare_data_file_tolerance = 0.01
 compare_data_file_tolerance_type = relative
+
+! integrate columns x:y in a data file as integral(y(x), x, x(1), x(end))
+integrate_line_file            = PartAnalyze.csv  ! data file name
+integrate_line_columns         = 0:1              ! columns x:y ("002-PCoupled" over time "001-TIME")
+integrate_line_integral_value  = 0.               ! should be exactly zero as a sin(x) is integrated over 2 periods
+integrate_line_tolerance_value = 1e-20            ! Tolerance
+integrate_line_tolerance_type  = absolute         ! relative or absolute

regressioncheck/NIG_PIC_poisson_Leapfrog/parallel_plates_AC/command_line.ini

@@ -1,2 +1,2 @@
-MPI=1
+MPI=1,2,5
 cmd_suffix=DSMC.ini

src/equations/poisson/equation.f90

@@ -59,8 +59,7 @@ SUBROUTINE DefineParametersEquation()
 IMPLICIT NONE
 !==================================================================================================================================
 CALL prms%SetSection("Equation")
-CALL prms%CreateIntOption(      'IniExactFunc'     , 'TODO-DEFINE-PARAMETER\n'//&
-                                                     'Define exact function necessary for linear scalar advection')
+CALL prms%CreateIntOption(      'IniExactFunc'     , 'Define exact function necessary for linear scalar advection')
 CALL prms%CreateRealArrayOption('RefState'         , 'State(s) for electric potential (amplitude, frequency and phase shift).', multiple=.TRUE., no=3)
 CALL prms%CreateRealArrayOption('IniWavenumber'    , 'TODO-DEFINE-PARAMETER' , '1. , 1. , 1.')
 CALL prms%CreateRealArrayOption('IniCenter'        , 'TODO-DEFINE-PARAMETER' , '0. , 0. , 0.')
@@ -76,10 +75,10 @@ SUBROUTINE DefineParametersEquation()
                                                      'Modified for curved and shape-function influence (c*dt*SafetyFactor+r_cutoff)' , '1.0')
 
 ! Special BC with linear potential ramp (constant in time)
-CALL prms%CreateRealArrayOption('LinPhiBasePoint'  , 'Origin of coordinate system for linear potential ramp for BoundaryType = (/2,1/)'       , no=3 )
-CALL prms%CreateRealArrayOption('LinPhiNormal'     , 'Normal vector of coordinate system for linear potential ramp for BoundaryType = (/2,1/)', no=3 )
-CALL prms%CreateRealOption(     'LinPhiHeight'     , 'Interval for ramping from 0 to LinPhi potential ramp for BoundaryType = (/2,1/)' )
-CALL prms%CreateRealOption(     'LinPhi'           , 'Target potential value for ramping from 0 for BoundaryType = (/2,1/)' )
+CALL prms%CreateRealArrayOption('LinPhiBasePoint'  , 'Origin(s) of coordinate system for linear potential ramp for BoundaryType = (/7,1/)'       , multiple=.TRUE. , no=3 )
+CALL prms%CreateRealArrayOption('LinPhiNormal'     , 'Normal(s) vector of coordinate system for linear potential ramp for BoundaryType = (/7,1/)', multiple=.TRUE. , no=3 )
+CALL prms%CreateRealOption(     'LinPhiHeight'     , 'Interval(s) for ramping from 0 to LinPhi potential ramp for BoundaryType = (/7,1/)'        , multiple=.TRUE.        )
+CALL prms%CreateRealOption(     'LinPhi'           , 'Potential(s) for ramping from 0 for BoundaryType = (/7,1/)'                                , multiple=.TRUE.        )
 
 #if defined(PARTICLES)
 ! Special BC with floating potential that is defined by the absorbed power of the charged particles
@@ -127,7 +126,8 @@ SUBROUTINE InitEquation()
 INTEGER            :: i,BCType,BCState
 CHARACTER(LEN=255) :: BCName
 INTEGER            :: nRefStateMax
-LOGICAL            :: LinPhiAlreadySet,PCouplAlreadySet
+INTEGER            :: nLinState,nLinStateMax
+LOGICAL            :: PCouplAlreadySet
 #if defined(PARTICLES)
 CHARACTER(255)     :: ContainerName
 LOGICAL            :: CoupledPowerPotentialExists
@@ -146,14 +146,14 @@ SUBROUTINE InitEquation()
 
 ! Sanity Check BCs
 nRefStateMax = 0
+nLinStateMax = 0
 ! Defaults
 #if defined(PARTICLES)
 #if USE_LOADBALANCE
 IF(.NOT.PerformLoadBalance)&
 #endif /*USE_LOADBALANCE*/
     CalcPCouplElectricPotential=.FALSE.
 #endif /*defined(PARTICLES)*/
-LinPhiAlreadySet=.FALSE.
 PCouplAlreadySet=.FALSE.
 DO i=1,nBCs
   BCType  = BoundaryType(i,BC_TYPE)
@@ -168,20 +168,9 @@ SUBROUTINE InitEquation()
     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
     nRefStateMax = MAX(nRefStateMax,BCState)
+  ELSEIF(BCType.EQ.7.AND.BCState.GT.0)THEN
+    nLinStateMax = MAX(nLinStateMax,BCState)
   ELSEIF(BCType.EQ.2)THEN
-    ! Special BC with linear potential ramp (constant in time)
-    ! Only for BoundaryType = (/2,1/)
-    IF(BCState.EQ.1)THEN
-      ! Do not set twice
-      IF(LinPhiAlreadySet) CYCLE
-      LinPhiAlreadySet=.TRUE.
-      ! Read linear potential parameters
-      LinPhiBasePoint = GETREALARRAY('LinPhiBasePoint',3)
-      LinPhiNormal    = GETREALARRAY('LinPhiNormal',3)
-      LinPhiNormal    = UNITVECTOR(LinPhiNormal)
-      LinPhiHeight    = GETREAL('LinPhiHeight')
-      LinPhi          = GETREAL('LinPhi')
-    END IF ! BCState.EQ.1
 
 #if defined(PARTICLES)
     ! Special BC with that adjusts a floating potential in order to meet a user-specified power input
@@ -233,11 +222,36 @@ SUBROUTINE InitEquation()
   END IF
 END DO
 
+! Read linear potential boundaries
+nLinState = CountOption('LinPhi')
+IF(nLinStateMax.GT.nLinState)THEN
+  SWRITE(*,'(A,I0)') "nLinStateMax: ",nLinStateMax
+  SWRITE(*,'(A,I0,A)') "   nLinState: ",nLinState," (number of times LinPhi = X occurrences in the parameter file)"
+  CALL abort(__STAMP__&
+      ,'nLinStateMax > nLinState: The given LinState number for boundary type 7 is larger than the supplied LinPhi values. '//&
+       'Define the correct number of LinPhi via, e.g., \n\n  LinPhi = 120.0 ! LinPhi Nbr 1: Voltage\n  LinPhi = 500.0 '//&
+       '! LinPhi Nbr 2: Voltage\n and the corresponding LinPhiBasePoint, LinPhiNormal and LinPhiHeight parameters')
+END IF ! nLinStateMax.GT.nLinState
+IF(nLinState .GT. 0)THEN
+  ALLOCATE(LinPhiBasePoint(3,nLinState))
+  ALLOCATE(LinPhiNormal(3,nLinState))
+  ALLOCATE(LinPhiHeight(nLinState))
+  ALLOCATE(LinPhi(nLinState))
+  DO iState=1,nLinState
+    ! Read linear potential parameters
+    LinPhiBasePoint(1:3,iState) = GETREALARRAY('LinPhiBasePoint',3)
+    LinPhiNormal(1:3,iState)    = GETREALARRAY('LinPhiNormal',3)
+    LinPhiNormal(1:3,iState)    = UNITVECTOR(LinPhiNormal(1:3,iState))
+    LinPhiHeight(iState)        = GETREAL('LinPhiHeight')
+    LinPhi(iState)              = GETREAL('LinPhi')
+  END DO
+END IF
+
 ! Read Boundary information / RefStates / perform sanity check
 nRefState=CountOption('RefState')
 IF(nRefStateMax.GT.nRefState)THEN
   SWRITE(*,'(A,I0)') "nRefStateMax: ",nRefStateMax
-  SWRITE(*,'(A,I0,A)') "   nRefState: ",nRefState," (number of times RefState = (/x,x,x/) occurs in the parameter file)"
+  SWRITE(*,'(A,I0,A)') "   nRefState: ",nRefState," (number of times RefState = (/x,x,x/) occurrences in the parameter file)"
   CALL abort(__STAMP__&
       ,'nRefStateMax > nRefState: The given RefState number for boundary type 5 is larger than the supplied RefStates. '//&
        'Define the correct number of RefStates via, e.g., \n\n  RefState = (/100.0 , 13.56E6 , -1.57079632679/) '//&
@@ -298,7 +312,7 @@ SUBROUTINE InitEquation()
 END SUBROUTINE InitEquation
 
 
-SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState)
+SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState,iLinState)
 !===================================================================================================================================
 ! Specifies all the initial conditions. The state in conservative variables is returned.
 !===================================================================================================================================
@@ -319,15 +333,32 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState)
 INTEGER,INTENT(IN)              :: ExactFunction    ! determines the exact function
 INTEGER,INTENT(IN),OPTIONAL     :: ElemID           ! ElemID
 REAL,INTENT(IN),OPTIONAl        :: t                ! time
-INTEGER,INTENT(IN),OPTIONAL     :: iRefState        ! ElemID
+INTEGER,INTENT(IN),OPTIONAL     :: iRefState        ! i-th reference state
+INTEGER,INTENT(IN),OPTIONAL     :: iLinState        ! i-th linear potential state
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! OUTPUT VARIABLES
 REAL,INTENT(OUT)                :: Resu(1:PP_nVar)    ! state in conservative variables
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
-REAL                            :: Omega,r1,r2,r_2D,r_3D,r_bary,cos_theta,eps1,eps2,xi
+REAL                            :: Omega,r1,r2,r_2D,r_3D,r_bary,cos_theta,eps1,eps2,xi,a(3),b(3)
 !===================================================================================================================================
 SELECT CASE (ExactFunction)
+CASE(-3) ! linear function with base point, normal vector and heigh: Requires BoundaryType = (/7,X/)
+  ASSOCIATE( height => LinPhiHeight(iLinState)    ,&
+             phi    => LinPhi(iLinState)          )
+    a = x(1:3) - LinPhiBasePoint(1:3,iLinState)
+    b = LinPhiNormal(1:3,iLinState)
+    xi = DOT_PRODUCT(a,b)
+    IF(xi.GT.0.)THEN
+      IF(xi.GT.height)THEN
+        Resu(:)=phi
+      ELSE
+        Resu(:)=phi*xi/height
+      END IF ! x(3)
+    ELSE
+      Resu(:)=0.0
+    END IF ! xi.GT.0
+  END ASSOCIATE
 CASE(-2) ! Signal without zero-crossing (always positive or negative), otherwise  like CASE(-1):
          ! Amplitude, Frequency and Phase Shift supplied by RefState
   ! RefState(1,iRefState): amplitude
@@ -344,24 +375,6 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState)
   Resu(:) = RefState(1,iRefState)*COS(Omega*t+RefState(3,iRefState))
 CASE(0) ! constant 0.
     Resu(:)=0.
-CASE(1) ! linear function with base point, normal vector and heigh: Requires BoundaryType = (/2,1/)
-  ASSOCIATE( origin => LinPhiBasePoint ,&
-             normal => LinPhiNormal    ,&
-             height => LinPhiHeight    ,&
-             phi    => LinPhi          )
-    ASSOCIATE( xVec => x(1:3) - origin(1:3) )
-      xi = DOT_PRODUCT(xVec,normal)
-      IF(xi.GT.0.)THEN
-        IF(xi.GT.height)THEN
-          Resu(:)=phi
-        ELSE
-          Resu(:)=phi*xi/height
-        END IF ! x(3)
-      ELSE
-        Resu(:)=0.0
-      END IF ! xi.GT.0
-    END ASSOCIATE
-  END ASSOCIATE
 #if defined(PARTICLES)
 CASE(2) ! Floating voltage that depends on the coupled power. User-specified power input as target value sets the potential.
   Resu(:) = CoupledPowerPotential(2)
@@ -523,9 +536,7 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState)
     IF((r1.LE.0.0).OR.(r2.LE.0.0))THEN
       SWRITE(*,*) "r1=",r1
       SWRITE(*,*) "r2=",r2
-      CALL abort(&
-          __STAMP__&
-          ,'ExactFunc=400: Point source in dielectric region. Cannot evaluate the exact function at the singularity!')
+      CALL abort(__STAMP__,'ExactFunc=400: Point source in dielectric region. Cannot evaluate the exact function at the singularity!')
     END IF
     resu(1:PP_nVar) = (1./eps1)*(&
                                    1./r1 + ((eps1-eps2)/(eps1+eps2))*&
@@ -539,7 +550,7 @@ SUBROUTINE ExactFunc(ExactFunction,x,resu,t,ElemID,iRefState)
   END IF
 
 CASE DEFAULT
-  CALL abort(__STAMP__,'Exactfunction not specified!')
+  CALL abort(__STAMP__,'Exactfunction not specified!', IntInfoOpt=ExactFunction)
 END SELECT ! ExactFunction
 
 END SUBROUTINE ExactFunc
@@ -790,6 +801,10 @@ SUBROUTINE FinalizeEquation()
 SDEALLOCATE(E)
 SDEALLOCATE(Et)
 SDEALLOCATE(RefState)
+SDEALLOCATE(LinPhiBasePoint)
+SDEALLOCATE(LinPhiNormal)
+SDEALLOCATE(LinPhiHeight)
+SDEALLOCATE(LinPhi)
 END SUBROUTINE FinalizeEquation
 
 END MODULE MOD_Equation

src/equations/poisson/equation_vars.f90

@@ -61,10 +61,10 @@ MODULE MOD_Equation_Vars
 REAL,ALLOCATABLE  :: RefState(:,:) !< refstates in primitive variables (as read from ini file)
 
 ! Special BC with linear potential ramp (constant in time)
-REAL              :: LinPhiBasePoint(3)
-REAL              :: LinPhiNormal(3)
-REAL              :: LinPhiHeight
-REAL              :: LinPhi
+REAL,ALLOCATABLE :: LinPhiBasePoint(:,:)
+REAL,ALLOCATABLE :: LinPhiNormal(:,:)
+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

src/hdg/hdg.f90

@@ -206,7 +206,7 @@ SUBROUTINE InitHDG()
   BCType =BoundaryType(BC(SideID),BC_TYPE)
   BCState=BoundaryType(BC(SideID),BC_STATE)
   SELECT CASE(BCType)
-  CASE(2,4,5,6) ! Dirichlet
+  CASE(2,4,5,6,7) ! Dirichlet
     nDirichletBCsides=nDirichletBCsides+1
   CASE(10,11) ! Neumann
     nNeumannBCsides=nNeumannBCsides+1
@@ -232,13 +232,15 @@ SUBROUTINE InitHDG()
   BCType =BoundaryType(BC(SideID),BC_TYPE)
   BCState=BoundaryType(BC(SideID),BC_STATE)
   SELECT CASE(BCType)
-  CASE(2,4,5,6) ! Dirichlet
+  CASE(2,4,5,6,7) ! Dirichlet
     nDirichletBCsides=nDirichletBCsides+1
     DirichletBC(nDirichletBCsides)=SideID
     MaskedSide(SideID)=.TRUE.
   CASE(10,11) !Neumann,
     nNeumannBCsides=nNeumannBCsides+1
     NeumannBC(nNeumannBCsides)=SideID
+  CASE DEFAULT ! unknown BCType
+    CALL abort(__STAMP__,' unknown BC Type in hdg.f90!',IntInfoOpt=BCType)
   END SELECT ! BCType
 END DO
 
@@ -379,7 +381,7 @@ SUBROUTINE InitZeroPotential()
   SELECT CASE(BCType)
   CASE(1) ! periodic
     ! do nothing
-  CASE(2,4,5,6) ! Dirichlet
+  CASE(2,4,5,6,7) ! 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
@@ -628,16 +630,21 @@ 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)
+    CASE(5) ! 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)
       END DO; END DO !p,q
-    CASE(6) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional)
+    CASE(6) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional) for reference state (without zero crossing)
       DO q=0,PP_N; DO p=0,PP_N
         r=q*(PP_N+1) + p+1
         CALL ExactFunc(-2,Face_xGP(:,p,q,SideID),lambda(iVar,r:r,SideID),t=time,iRefState=BCState)
       END DO; END DO !p,q
+    CASE(7) ! exact BC = Dirichlet BC !! ExactFunc via LinState (time is optional)for linear potential function
+      DO q=0,PP_N; DO p=0,PP_N
+        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
     END SELECT ! BCType
   END DO !BCsideID=1,nDirichletBCSides
 #if (PP_nVar!=1)
@@ -898,16 +905,21 @@ SUBROUTINE HDGNewton(time,U_out,td_iter,ForceCGSolverIteration_opt)
       r=q*(PP_N+1) + p+1
       lambda(PP_nVar,r:r,SideID)= 0.
     END DO; END DO !p,q
-  CASE(5) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional)
+  CASE(5) ! 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(PP_nVar,r:r,SideID),t=time,iRefState=BCState)
     END DO; END DO !p,q
-  CASE(6) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional)
+  CASE(6) ! exact BC = Dirichlet BC !! ExactFunc via RefState (time is optional) for reference state (without zero crossing)
     DO q=0,PP_N; DO p=0,PP_N
       r=q*(PP_N+1) + p+1
       CALL ExactFunc(-2,Face_xGP(:,p,q,SideID),lambda(PP_nVar,r:r,SideID),t=time,iRefState=BCState)
     END DO; END DO !p,q
+  CASE(7) ! exact BC = Dirichlet BC !! ExactFunc via LinState (time is optional)for linear potential function
+    DO q=0,PP_N; DO p=0,PP_N
+      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
   END SELECT ! BCType
 END DO !BCsideID=1,nDirichletBCSides
 

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) !Dirichlet
+          CASE(2,4,5,6,7) !Dirichlet
             ! do not consider this side
           CASE(10,11) !Neumann
             HDGSides = HDGSides + 1