merge with master. all the changes, the merge has removed...

.directory

@@ -0,0 +1,6 @@
+[Dolphin]
+Timestamp=2016,11,16,13,25,9
+Version=3
+
+[Settings]
+HiddenFilesShown=true

src/equations/maxwell/equation.f90

@@ -279,7 +279,7 @@ SUBROUTINE ExactFunc(ExactFunction,t,tDeriv,x,resu)
 REAL, PARAMETER                 :: B0G=1.0,g=3236.706462           ! aux. Constants for Gyrotron
 REAL, PARAMETER                 :: k0=3562.936537,h=1489.378411    ! aux. Constants for Gyrotron
 REAL, PARAMETER                 :: omegaG=3.562936537e+3           ! aux. Constants for Gyrotron
-REAL                            :: spatialWindow,tShift            !> electromagnetic wave shaping vars
+REAL                            :: spatialWindow,tShift,tShiftBC!> electromagnetic wave shaping vars
 REAL                            :: timeFac,temporalWindow
 INTEGER, PARAMETER              :: mG=34,nG=19                     ! aux. Constants for Gyrotron
 REAL                            :: eta, kx,ky,kz
@@ -518,19 +518,20 @@ SUBROUTINE ExactFunc(ExactFunction,t,tDeriv,x,resu)
     RETURN
   END IF
   ! IC or BC
-  tShift=t-ABS(WaveBasePoint(BeamIdir3))/c
+  tShift=t-ABS(WaveBasePoint(BeamIdir3))/c ! shift to wave base point position
   IF(t.LT.dt)THEN ! initial condiction: IC
     spatialWindow = EXP(    -((x(BeamIdir1)-WaveBasePoint(BeamIdir1))**2+                  &
                               (x(BeamIdir2)-WaveBasePoint(BeamIdir2))**2)*omega_0_2inv     &
                          -0.5*(x(BeamIdir3)-WaveBasePoint(BeamIdir3))**2/((sigma_t*c)**2)  )
     timeFac=COS(BeamWaveNumber*DOT_PRODUCT(WaveVector,x-WaveBasePoint)-BeamOmegaW*tShift)
     resu(1:3)=BeamAmpFac*spatialWindow*E_0*timeFac
   ELSE ! boundary condiction: BC
+    tShiftBC=t+(WaveBasePoint(BeamIdir3)-x(3))/c ! shift to wave base point position
     ! intensity * Gaussian filter in transversal and longitudinal direction
     spatialWindow = EXP(-((x(BeamIdir1)-WaveBasePoint(BeamIdir1))**2+&
                           (x(BeamIdir2)-WaveBasePoint(BeamIdir2))**2)*omega_0_2inv)
     timeFac =COS(BeamWaveNumber*DOT_PRODUCT(WaveVector,x-WaveBasePoint)-BeamOmegaW*tShift)
-    temporalWindow=EXP(-0.5*(tShift/sigma_t)**2)
+    temporalWindow=EXP(-0.5*(tShiftBC/sigma_t)**2)
     resu(1:3)=BeamAmpFac*spatialWindow*E_0*timeFac*temporalWindow
   END IF
   resu(4:6)=c_inv*CROSS(WaveVector,resu(1:3)) 

src/interpolation/eval_xyz.f90

@@ -159,7 +159,7 @@ SUBROUTINE eval_xyz_curved(x_in,NVar,N_in,U_In,U_Out,ElemID,PartID)
 END SUBROUTINE eval_xyz_curved
 
 
-SUBROUTINE eval_xyz_elemcheck(x_in,xi,ElemID,DoReUseMap)
+SUBROUTINE eval_xyz_elemcheck(x_in,xi,ElemID,DoReUseMap,ForceMode)
 !===================================================================================================================================
 ! interpolate a 3D tensor product Lagrange basis defined by (N_in+1) 1D interpolation point positions x
 ! first get xi,eta,zeta from x,y,z...then do tensor product interpolation
@@ -180,6 +180,7 @@ SUBROUTINE eval_xyz_elemcheck(x_in,xi,ElemID,DoReUseMap)
 INTEGER,INTENT(IN)          :: ElemID                                 ! elem index
 REAL,INTENT(IN)             :: x_in(3)                                  ! physical position of particle 
 LOGICAL,INTENT(IN),OPTIONAL :: DoReUseMap
+LOGICAL,INTENT(IN),OPTIONAL :: ForceMode
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! OUTPUT VARIABLES
 REAL,INTENT(INOUT)          :: xi(1:3)
@@ -197,11 +198,11 @@ SUBROUTINE eval_xyz_elemcheck(x_in,xi,ElemID,DoReUseMap)
 END IF
 
 IF(CurvedElem(ElemID))THEN
-  CALL RefElemNewton(Xi,X_In,wBaryCL_NGeo,XiCL_NGeo,XCL_NGeo(:,:,:,:,ElemID),dXCL_NGeo(:,:,:,:,:,ElemID),NGeo,ElemID,Mode=2)
+  CALL RefElemNewton(Xi,X_In,wBaryCL_NGeo,XiCL_NGeo,XCL_NGeo(:,:,:,:,ElemID),dXCL_NGeo(:,:,:,:,:,ElemID),NGeo,ElemID,Mode=iMode)
 ELSE
   ! fill dummy XCL_NGeo1
   IF(NGeo.EQ.1)THEN
-    CALL RefElemNewton(Xi,X_In,wBaryCL_NGeo,XiCL_NGeo,XCL_NGeo(:,:,:,:,ElemID),dXCL_NGeo(:,:,:,:,:,ElemID),NGeo,ElemID,Mode=2)
+    CALL RefElemNewton(Xi,X_In,wBaryCL_NGeo,XiCL_NGeo,XCL_NGeo(:,:,:,:,ElemID),dXCL_NGeo(:,:,:,:,:,ElemID),NGeo,ElemID,Mode=iMode)
   ELSE
     XCL_NGeo1(1:3,0,0,0) = XCL_NGeo(1:3, 0  , 0  , 0  ,ElemID)
     XCL_NGeo1(1:3,1,0,0) = XCL_NGeo(1:3,NGeo, 0  , 0  ,ElemID)
@@ -220,7 +221,7 @@ SUBROUTINE eval_xyz_elemcheck(x_in,xi,ElemID,DoReUseMap)
     dXCL_NGeo1(1:3,1:3,1,0,1) = dXCL_NGeo(1:3,1:3,NGeo, 0  ,NGeo,ElemID)
     dXCL_NGeo1(1:3,1:3,0,1,1) = dXCL_NGeo(1:3,1:3, 0  ,NGeo,NGeo,ElemID)
     dXCL_NGeo1(1:3,1:3,1,1,1) = dXCL_NGeo(1:3,1:3,NGeo,NGeo,NGeo,ElemID)
-    CALL RefElemNewton(Xi,X_In,wBaryCL_NGeo1,XiCL_NGeo1,XCL_NGeo1,dXCL_NGeo1,1,ElemID,Mode=2)
+    CALL RefElemNewton(Xi,X_In,wBaryCL_NGeo1,XiCL_NGeo1,XCL_NGeo1,dXCL_NGeo1,1,ElemID,Mode=iMode)
   END IF
 END IF
 
@@ -434,11 +435,12 @@ SUBROUTINE RefElemNewton(Xi,X_In,wBaryCL_N_In,XiCL_N_In,XCL_N_In,dXCL_N_In,N_In,
 REAL,INTENT(INOUT)               :: Xi(3) ! position in reference element
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
-REAL                             :: Lag(1:3,0:N_In), F(1:3)
+REAL                             :: Lag(1:3,0:N_In), F(1:3),Xi_Old(1:3)
 INTEGER                          :: NewTonIter,i,j,k
 REAL                             :: deltaXi(1:3),deltaXi2
 REAL                             :: Jac(1:3,1:3),sdetJac,sJac(1:3,1:3)
-REAL                             :: buff,buff2
+REAL                             :: buff,buff2, Norm_F, Norm_F_old,lambda
+INTEGER                          :: iArmijo
 !===================================================================================================================================
 
 
@@ -463,6 +465,8 @@ SUBROUTINE RefElemNewton(Xi,X_In,wBaryCL_N_In,XiCL_N_In,XCL_N_In,dXCL_N_In,N_In,
   deltaXi2=1. !HUGE(1.0)
 END IF
 
+Norm_F=DOT_PRODUCT(F,F)
+Norm_F_old=Norm_F
 NewtonIter=0
 !abortCrit=ElemRadiusN_in(ElemID)*ElemRadiusN_in(ElemID)*RefMappingEps
 DO WHILE((deltaXi2.GT.RefMappingEps).AND.(NewtonIter.LT.100))
@@ -504,9 +508,11 @@ SUBROUTINE RefElemNewton(Xi,X_In,wBaryCL_N_In,XiCL_N_In,XCL_N_In,dXCL_N_In,N_In,
   ! Iterate Xi using Newton step
   ! Use FAIL
   !Xi = Xi - MATMUL(sJac,F)
+
+  ! Armijo step size control
   deltaXi=MATMUL(sJac,F)
-  Xi = Xi - deltaXI!MATMUL(sJac,F)
   deltaXi2=DOT_PRODUCT(deltaXi,deltaXi)
+  Xi_Old=Xi
 
   Norm_F_old=Norm_F
   Norm_F=Norm_F*2.
@@ -549,31 +555,16 @@ SUBROUTINE RefElemNewton(Xi,X_In,wBaryCL_N_In,XiCL_N_In,XCL_N_In,dXCL_N_In,N_In,
       IF(PRESENT(PartID)) IPWRITE(UNIT_stdOut,*) ' implicit?', PartisImplicit(PartID)
       IF(PRESENT(PartID)) IPWRITE(UNIT_stdOut,*) ' last?', LastPartPos(PartID,1:3)
 #endif
-      CALL abort(&
-__STAMP__&
-,'Particle Not inSide of Element, ElemID,',ElemID)
+        CALL abort(&
+  __STAMP__&
+  ,'Particle Not inSide of Element, ElemID,',ElemID)
     ELSE
       EXIT
     END IF
   END IF
-
-  ! Compute function value
-  CALL LagrangeInterpolationPolys(Xi(1),N_In,XiCL_N_in,wBaryCL_N_in,Lag(1,:))
-  CALL LagrangeInterpolationPolys(Xi(2),N_In,XiCL_N_in,wBaryCL_N_in,Lag(2,:))
-  CALL LagrangeInterpolationPolys(Xi(3),N_In,XiCL_N_in,wBaryCL_N_in,Lag(3,:))
-  ! F(xi) = x(xi) - x_in
-  F=-x_in ! xRp
-  DO k=0,N_In
-    DO j=0,N_In
-      buff=Lag(2,j)*Lag(3,k)
-      DO i=0,N_In
-        buff2=Lag(1,i)*buff
-        F=F+XCL_N_in(:,i,j,k)*buff2
-      END DO !l=0,N_In
-    END DO !i=0,N_In
-  END DO !j=0,N_In
+
 END DO !newton
-!print*,'newton iter', newtoniter
+
 
 END SUBROUTINE RefElemNewton
 

src/mesh/mapping.f90

@@ -152,21 +152,21 @@ SUBROUTINE InitMappings(N_in,V2S,V2SIJK,V2S2,CV2S,S2V,S2V2,CS2V,FS2M)
 DO j=0,N_in; DO i=0,N_in
   DO f=0,4
     DO s=1,6
-      S2V2(:,i,j,f,s) = SideToVol2(i,j,f,s)
+      S2V2(:,i,j,f,s) = SideToVol2(N_in,i,j,f,s)
     END DO
   END DO
 END DO; END DO
 
 DO j=0,N_in; DO i=0,N_in
   DO s=1,6
-    CS2V(:,i,j,s) = CGNS_SideToVol2(i,j,s)
+    CS2V(:,i,j,s) = CGNS_SideToVol2(N_in,i,j,s)
   END DO
 END DO; END DO
 
 ! Flip_S2M
 DO j=0,N_in; DO i=0,N_in
   DO f=0,4
-    FS2M(:,i,j,f) = Flip_S2M(i,j,f)
+    FS2M(:,i,j,f) = Flip_S2M(N_in,i,j,f)
   END DO
 END DO; END DO
 
@@ -187,7 +187,7 @@ SUBROUTINE InitMappings(N_in,V2S,V2SIJK,V2S2,CV2S,S2V,S2V2,CS2V,FS2M)
 END SUBROUTINE InitMappings
 
 
-FUNCTION Flip_S2M(p, q, flip)
+FUNCTION Flip_S2M(N_in,p, q, flip)
 !===================================================================================================================================
 ! Transforms Coordinates from RHS of Slave to RHS of Master 
 !    input: p,q in Slave-RHS, flip;  
@@ -198,7 +198,7 @@ FUNCTION Flip_S2M(p, q, flip)
 ! IMPLICIT VARIABLE HANDLING
 IMPLICIT NONE
 ! INPUT VARIABLES
-INTEGER,INTENT(IN) :: p,q,flip
+INTEGER,INTENT(IN) :: N_in,p,q,flip
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! INPUT/OUTPUT VARIABLES
 !-----------------------------------------------------------------------------------------------------------------------------------
@@ -213,16 +213,16 @@ FUNCTION Flip_S2M(p, q, flip)
   CASE(1) 
     Flip_S2M = (/     q,     p/)
   CASE(2) 
-    Flip_S2M = (/PP_N-p,     q/)
+    Flip_S2M = (/N_in-p,     q/)
   CASE(3) 
-    Flip_S2M = (/PP_N-q,PP_N-p/)
+    Flip_S2M = (/N_in-q,N_in-p/)
   CASE(4) 
-    Flip_S2M = (/     p,PP_N-q/)
+    Flip_S2M = (/     p,N_in-q/)
 END SELECT
 END FUNCTION Flip_S2M
 
 
-FUNCTION Flip_M2S(p, q, flip)
+FUNCTION Flip_M2S(N_in,p, q, flip)
 !===================================================================================================================================
 ! Transforms Coordinates from RHS of Master to RHS of Slave
 !   actualy this is the same function as Flip_S2M
@@ -232,7 +232,7 @@ FUNCTION Flip_M2S(p, q, flip)
 ! IMPLICIT VARIABLE HANDLING
 IMPLICIT NONE
 ! INPUT VARIABLES
-INTEGER,INTENT(IN) :: p,q,flip
+INTEGER,INTENT(IN) :: N_in,p,q,flip
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! INPUT/OUTPUT VARIABLES
 !-----------------------------------------------------------------------------------------------------------------------------------
@@ -241,7 +241,7 @@ FUNCTION Flip_M2S(p, q, flip)
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
 !===================================================================================================================================
-Flip_M2S=Flip_S2M(p,q,flip)
+Flip_M2S=Flip_S2M(N_in,p,q,flip)
 END FUNCTION Flip_M2S
 
 
@@ -322,7 +322,7 @@ FUNCTION CGNS_SideToVol(l, p, q, locSideID)
 END FUNCTION CGNS_SideToVol
 
 
-FUNCTION CGNS_SideToVol2(p, q, locSideID)
+FUNCTION CGNS_SideToVol2(N_in,p, q, locSideID)
 !===================================================================================================================================
 ! Transforms RHS-Coordinates of Side (CGNS-Notation) into Volume-Coordinates 
 ! input: l, p,q, locSideID
@@ -335,7 +335,7 @@ FUNCTION CGNS_SideToVol2(p, q, locSideID)
 ! IMPLICIT VARIABLE HANDLING
 IMPLICIT NONE
 ! INPUT VARIABLES
-INTEGER,INTENT(IN) :: p,q,locSideID
+INTEGER,INTENT(IN) :: p,q,locSideID,N_in
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! INPUT/OUTPUT VARIABLES
 !-----------------------------------------------------------------------------------------------------------------------------------
@@ -352,7 +352,7 @@ FUNCTION CGNS_SideToVol2(p, q, locSideID)
   CASE(ETA_MINUS)  
     CGNS_SideToVol2 = (/p,q/)
   CASE(ETA_PLUS)   
-    CGNS_SideToVol2 = (/PP_N-p,q/)
+    CGNS_SideToVol2 = (/N_in-p,q/)
   CASE(ZETA_MINUS) 
     CGNS_SideToVol2 = (/q,p/)
   CASE(ZETA_PLUS)  
@@ -384,7 +384,7 @@ FUNCTION VolToSide(i,j,k, flip, locSideID)
 INTEGER,DIMENSION(3) :: pq
 !===================================================================================================================================
 pq = CGNS_VolToSide(i,j,k,locSideID)
-VolToSide(1:2) = Flip_S2M(pq(1),pq(2),flip)
+VolToSide(1:2) = Flip_S2M(PP_N,pq(1),pq(2),flip)
 VolToSide(3) = pq(3)
 END FUNCTION VolToSide
 
@@ -449,7 +449,7 @@ FUNCTION VolToSide2(ijk1,ijk2, flip, locSideID)
 CASE(ZETA_MINUS,ZETA_PLUS)
   pq = CGNS_VolToSide(ijk1,ijk2,0,locSideID)
 END SELECT
-VolToSide2(1:2) = Flip_S2M(pq(1),pq(2),flip)
+VolToSide2(1:2) = Flip_S2M(PP_N,pq(1),pq(2),flip)
 END FUNCTION VolToSide2
 
 FUNCTION SideToVol(l, p, q, flip, locSideID)
@@ -475,7 +475,7 @@ FUNCTION SideToVol(l, p, q, flip, locSideID)
 ! LOCAL VARIABLES
 INTEGER,DIMENSION(2) :: pq
 !===================================================================================================================================
-pq = Flip_M2S(p,q,flip)
+pq = Flip_M2S(PP_N,p,q,flip)
 SideToVol = CGNS_SideToVol(l,pq(1),pq(2),locSideID)
 END FUNCTION SideToVol
 
@@ -505,7 +505,7 @@ FUNCTION SideToAdjointLocSide(p, q, flip, locSideID)
 INTEGER,DIMENSION(2) :: pq
 !===================================================================================================================================
 
-pq = SideToVol2(p, q, flip, locSideID)
+pq = SideToVol2(PP_N,p, q, flip, locSideID)
 SELECT CASE(locSideID)
 CASE(XI_MINUS)
   IF(pq(1).EQ.0)THEN
@@ -538,7 +538,7 @@ FUNCTION SideToAdjointLocSide(p, q, flip, locSideID)
 END FUNCTION SideToAdjointLocSide
 
 
-FUNCTION SideToVol2(p, q, flip, locSideID)
+FUNCTION SideToVol2(N_in,p, q, flip, locSideID)
 !===================================================================================================================================
 ! Transform RHS-Coordinates of Master to Volume-Coordinates. This is: SideToVol = CGNS_SideToVol(Flip_M2S(...))
 ! input: l, p,q, flip, locSideID
@@ -551,7 +551,7 @@ FUNCTION SideToVol2(p, q, flip, locSideID)
 ! IMPLICIT VARIABLE HANDLING
 IMPLICIT NONE
 ! INPUT VARIABLES
-INTEGER,INTENT(IN) :: p,q,flip,locSideID
+INTEGER,INTENT(IN) :: N_in,p,q,flip,locSideID
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! INPUT/OUTPUT VARIABLES
 !-----------------------------------------------------------------------------------------------------------------------------------
@@ -561,8 +561,8 @@ FUNCTION SideToVol2(p, q, flip, locSideID)
 ! LOCAL VARIABLES
 INTEGER,DIMENSION(2) :: pq
 !===================================================================================================================================
-pq = Flip_M2S(p,q,flip)
-SideToVol2 = CGNS_SideToVol2(pq(1),pq(2),locSideID)
+pq = Flip_M2S(N_in,p,q,flip)
+SideToVol2 = CGNS_SideToVol2(N_in,pq(1),pq(2),locSideID)
 END FUNCTION SideToVol2
 
 

src/mesh/mesh_readin.f90

@@ -807,9 +807,15 @@ SUBROUTINE ReadMesh(FileString)
       aSide%Ind=ABS(SideInfo(SIDE_ID,iSide))
       IF(oriented)THEN !oriented side
         aSide%flip=0
+#ifdef PARTICLES
+        aSide%BC_Alpha=99
+#endif /*PARTICLES*/
       ELSE !not oriented
         aSide%flip=MOD(Sideinfo(SIDE_Flip,iSide),10)
         IF((aSide%flip.LT.0).OR.(aSide%flip.GT.4)) STOP 'NodeID doesnt belong to side'
+#ifdef PARTICLES
+        aSide%BC_Alpha=-99
+#endif /*PARTICLES*/
       END IF
     ELSE !mortartype>0
       DO iMortar=1,aSide%nMortars
@@ -818,12 +824,14 @@ SUBROUTINE ReadMesh(FileString)
         IF(SideInfo(SIDE_ID,iSide).LT.0) STOP 'Problem in Mortar readin,should be flip=0'
         aSide%mortarSide(iMortar)%sp%flip=0
         aSide%mortarSide(iMortar)%sp%Ind=ABS(SideInfo(SIDE_ID,iSide))
+#ifdef PARTICLES
+        aSide%BC_Alpha=99
+#endif /*PARTICLES*/
       END DO !iMortar
     END IF
   END DO !i=1,locnSides
 END DO !iElem
 
-
 ! build up side connection
 DO iElem=FirstElemInd,LastElemInd
   aElem=>Elems(iElem)%ep
@@ -848,6 +856,9 @@ SUBROUTINE ReadMesh(FileString)
            (BoundaryType(aSide%BCindex,BC_TYPE).NE.100))THEN ! Reassignement from periodic to non-periodic
           doConnection=.FALSE.
           aSide%flip  =0
+#ifdef PARTICLES
+          aSide%BC_Alpha=0
+#endif /*PARTICLES*/
           IF(iMortar.EQ.0) aSide%mortarType  = 0
           IF(iMortar.EQ.0) aSide%nMortars    = 0
           elemID            = 0

src/mesh/mesh_vars.f90

@@ -206,6 +206,9 @@ MODULE MOD_Mesh_Vars
   INTEGER                      :: NbProc 
   INTEGER                      :: BCindex         !< index in BoundaryType array! 
   INTEGER                      :: flip 
+#ifdef PARTICLES
+  INTEGER                      :: BC_Alpha        !< inital value for periodic displacement before mapping in pos. bc-index range
+#endif /*PARTICLES*/
   INTEGER                      :: nMortars        !< number of slave mortar sides associated with master mortar
   INTEGER                      :: MortarType      !< type of mortar: Type1 : 1-4 , Type 2: 1-2 in eta, Type 2: 1-2 in xi
   TYPE(tSidePtr),POINTER       :: MortarSide(:)   !< array of side pointers to slave mortar sides