Merge branch '8-correct-sfr-geometry-issues-in-module-salt' into 'main'

#8: Correct a domain translation issue with TSPC tracking

See merge request dragon/5.1!26

3 files changed, 65 insertions, 15 deletions

diff --git a/Dragon/src/EVODRV.f b/Dragon/src/EVODRV.f
index afe2f35..cb9b40c 100644
--- a/Dragon/src/EVODRV.f
+++ b/Dragon/src/EVODRV.f
@@ -552,7 +552,11 @@
         ENDIF
         IF(ITIM.EQ.0) THEN
           IF(NTIM.GT.0) THEN
-            IF(XT(IP).LT.TIMES(NTIM)) CALL XABORT('EVODRV: INVALID X1.')
+            IF(XT(IP).LT.TIMES(NTIM)) THEN
+              WRITE(HSMG,'(19HEVODRV: INVALID XT(,I4,2H)=,1P,E12.4,2H <,
+     >        E12.4,1H.)') IP,XT(IP),TIMES(NTIM)
+              CALL XABORT(HSMG)
+            ENDIF
           ENDIF
           NTIM=NTIM+1
           TIMES(NTIM)=XT(IP)
@@ -603,7 +607,11 @@
         ENDIF
         IF(ITIM.EQ.0) THEN
           IF(NTIM.GT.0) THEN
-            IF(XTI.LT.TIMES(NTIM)) CALL XABORT('EVODRV: INVALID X1.')
+            IF(XTI.LT.TIMES(NTIM)) THEN
+              WRITE(HSMG,'(20HEVODRV: INVALID XTI=,1P,E12.4,2H <,E12.4,
+     >        1H.)') XTI,TIMES(NTIM)
+              CALL XABORT(HSMG)
+            ENDIF
           ENDIF
           NTIM=NTIM+1
           TIMES(NTIM)=XTI
diff --git a/Dragon/src/SAL_GEOMETRY_MOD.f90 b/Dragon/src/SAL_GEOMETRY_MOD.f90
index 4bce6e8..a56376f 100644
--- a/Dragon/src/SAL_GEOMETRY_MOD.f90
+++ b/Dragon/src/SAL_GEOMETRY_MOD.f90
@@ -154,6 +154,7 @@ CONTAINS
     ! local variable
     ! ***************
     INTEGER :: ELEM,OK,I,TYPE
+    REAL(PDB) :: X1,Y1,X2,Y2,XMIN,XMAX,YMIN,YMAX
     CHARACTER(LEN=4) :: HSYM
     REAL(PDB),PARAMETER :: CONV=PI/180._PDB
     INTEGER, PARAMETER :: FOUT =6
@@ -188,10 +189,48 @@ CONTAINS
     ENDIF
     ALLOCATE (GG%IPAR(NIPAR,GG%NB_ELEM), GG%RPAR(NRPAR,GG%NB_ELEM), STAT=OK)
     IF(OK/=0) CALL XABORT('SAL110: not enough memory I,R')
-
+    !
     !*    read surfacic file
     CALL SALINP(GG)
     !
+    ! translate the domain for cyclic cases
+    XMIN=1.E10_PDB; YMIN=1.E10_PDB; XMAX=-1.E10_PDB; YMAX=-1.E10_PDB;
+    DO ELEM=1,GG%NB_ELEM
+      TYPE=GG%IPAR(1,ELEM)
+      IF(TYPE==1) THEN
+        X1=GG%RPAR(1,ELEM); Y1=GG%RPAR(2,ELEM);
+        XMIN=MIN(XMIN,X1); YMIN=MIN(YMIN,Y1); XMAX=MAX(XMAX,X1); YMAX=MAX(YMAX,Y1);
+        X2=X1+GG%RPAR(3,ELEM); Y2=Y1+GG%RPAR(4,ELEM);
+        XMIN=MIN(XMIN,X2); YMIN=MIN(YMIN,Y2); XMAX=MAX(XMAX,X2); YMAX=MAX(YMAX,Y2);
+      ENDIF
+    ENDDO
+    IF(ABS(XMIN).LT.1.D-10) XMIN=0.D0
+    IF(ABS(YMIN).LT.1.D-10) YMIN=0.D0
+    SELECT CASE(TYPGEO)
+    CASE(4,9)
+      DO ELEM=1,GG%NB_ELEM
+        GG%RPAR(1,ELEM)=GG%RPAR(1,ELEM)-(XMIN+XMAX)/2.D0
+        GG%RPAR(2,ELEM)=GG%RPAR(2,ELEM)-(XMIN+XMAX)/2.D0
+      ENDDO
+      DO I=1,GG%NBBCDA
+        GG%BCDATAREAD(I)%BCDATA(1)=GG%BCDATAREAD(I)%BCDATA(1)-(XMIN+XMAX)/2.D0
+        GG%BCDATAREAD(I)%BCDATA(2)=GG%BCDATAREAD(I)%BCDATA(2)-(XMIN+XMAX)/2.D0
+        IF(ABS(GG%BCDATAREAD(I)%BCDATA(1)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(1)=0.D0
+        IF(ABS(GG%BCDATAREAD(I)%BCDATA(2)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(2)=0.D0
+      ENDDO
+    CASE(5,6,7,8,10,11,12)
+      DO ELEM=1,GG%NB_ELEM
+        GG%RPAR(1,ELEM)=GG%RPAR(1,ELEM)-XMIN
+        GG%RPAR(2,ELEM)=GG%RPAR(2,ELEM)-YMIN
+      ENDDO
+      DO I=1,GG%NBBCDA
+        GG%BCDATAREAD(I)%BCDATA(1)=GG%BCDATAREAD(I)%BCDATA(1)-XMIN
+        GG%BCDATAREAD(I)%BCDATA(2)=GG%BCDATAREAD(I)%BCDATA(2)-YMIN
+        IF(ABS(GG%BCDATAREAD(I)%BCDATA(1)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(1)=0.D0
+        IF(ABS(GG%BCDATAREAD(I)%BCDATA(2)).LT.1.D-10) GG%BCDATAREAD(I)%BCDATA(2)=0.D0
+      ENDDO
+    END SELECT
+    !
     !*    unite the ends of elements, redefine elements
     CALL SAL128(GG%RPAR,GG%IPAR,GG%NB_ELEM)
     !
@@ -1818,7 +1857,7 @@ CONTAINS
     !
     !---------------------------------------------------------------------
     !
-    USE SAL_GEOMETRY_TYPES,     ONLY : TYPGEO,ANGGEO,EPS,LX=>LENGTHX,LY=>LENGTHY,G_BC_TYPE
+    USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO,ANGGEO,EPS,LX=>LENGTHX,LY=>LENGTHY,G_BC_TYPE
     USE SAL_NUMERIC_MOD, ONLY : SAL141
     !****
     IMPLICIT NONE
@@ -1840,6 +1879,7 @@ CONTAINS
     REAL(PDB), DIMENSION(NPERIM,6) :: AUX_DIST
     INTEGER   :: I,J,K,M,ELEM,TYPBC,IBC,IDATA,NBE(6),OK,NAXES
     REAL(PDB) :: ANGLE,X,Y,D
+    CHARACTER(LEN=131) :: HSMG
     !****
     NAXES=0
     IF((TYPGEO==5).OR.(TYPGEO==6).OR.(TYPGEO==11)) THEN
@@ -1852,13 +1892,6 @@ CONTAINS
     !*    calculate number of elements on axes,
     !     and their distances to the origin of the axis
     NBE=0
-
-    DO I=1,NPERIM
-       ELEM=PERIM(I)
-       IBC=IBC2_ELEM(ELEM)
-       TYPBC=TYPE_BC2(IBC)
-       IDATA=IDATA_BC2(IBC)
-    ENDDO
     DO I=1,NPERIM
        ELEM=PERIM(I)
        IBC=IBC2_ELEM(ELEM)
@@ -2048,7 +2081,10 @@ CONTAINS
           CASE DEFAULT
              CALL XABORT('SAL130_10: boundary condition not implemented.')
        END SELECT
-       IF(M>=NAXES+1) CALL XABORT('SAL130_10: element not on axes')
+       IF(M>=NAXES+1) THEN
+         WRITE(HSMG,'(42HSAL130_10: element not on axes for typgeo=,i3,1h.)') TYPGEO
+         CALL XABORT(HSMG)
+       ENDIF
        NBE(M)=NBE(M)+1
        LIST_ELEMS(NBE(M),M)=ELEM
        !        sort the element list according their distance to
diff --git a/Dragon/src/SAL_TRAJECTORY_MOD.f90 b/Dragon/src/SAL_TRAJECTORY_MOD.f90
index ac1714d..9408d1b 100644
--- a/Dragon/src/SAL_TRAJECTORY_MOD.f90
+++ b/Dragon/src/SAL_TRAJECTORY_MOD.f90
@@ -68,7 +68,7 @@ CONTAINS
     !
     !---------------------------------------------------------------------
     !
-    USE SAL_GEOMETRY_TYPES,  ONLY : ISPEC
+    USE SAL_GEOMETRY_TYPES,  ONLY : ISPEC,TYPGEO
     USE SAL_TRACKING_TYPES, ONLY : NNN,NMAX2,ITRAC2,ANGTAB,ELMTAB,CNT,CNT0,NB_TOT,DNEW,DINIT, &
                                    NNEW,LNEW,IERR,LGMORE,DD0,NTRACK,EPS1,EX0,EY0,LGOK,IPART,DELX, &
                                    N_AXIS
@@ -110,7 +110,7 @@ CONTAINS
       ELSE
         WRITE(*,*) 'PPERIM_MAC2(N_AXIS+1),PPERIM_MAC2(N_AXIS) :',PPERIM_MAC2(N_AXIS+1),PPERIM_MAC2(N_AXIS)
         WRITE(*,*) 'DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1) :',DIST_AXIS(PPERIM_MAC2(N_AXIS+1)-1)
-        WRITE(*,*) 'DELX :',DELX,' RADIA=',RADIA
+        WRITE(*,*) 'DELX :',DELX,' RADIA=',RADIA,' TYPGEO=',TYPGEO
         CALL XABORT('SALTRA: Cant find entry point')
       ENDIF
     ENDIF
@@ -164,6 +164,7 @@ CONTAINS
     !
     !---------------------------------------------------------------------
     !
+    USE SAL_GEOMETRY_TYPES, ONLY : TYPGEO
     USE SAL_TRACKING_TYPES, ONLY : IPART,N_AXIS,DNEW,DELX,NNEW,LNEW,COSINE,AX, &
                                    AY,HX,HY,BX,BY,EX,EY
     INTEGER,   INTENT(IN)                  :: NPERIM
@@ -171,6 +172,7 @@ CONTAINS
     REAL(PDB), INTENT(IN), DIMENSION(:)    :: DIST_AXIS
     INTEGER,   INTENT(IN), DIMENSION(:,:)  :: IPAR
     INTEGER    :: I,J
+    CHARACTER(LEN=131) :: HSMG
     !***
     LNEW=0
     !*    compute crossed element
@@ -180,7 +182,11 @@ CONTAINS
           EXIT
        ENDIF
     ENDDO
-    IF(LNEW==0) CALL XABORT('SAL241_2: Error of distances on the axis')
+    IF(LNEW==0) THEN
+      WRITE(HSMG,'(52HSAL241_2: Error of distances on the axis for typgeo=, &
+      & i3,1h.)') TYPGEO
+      CALL XABORT(HSMG)
+    ENDIF
     !*    get entered node
     NNEW=IPAR(2,LNEW)
     IF(NNEW<0) NNEW=IPAR(3,LNEW)