Debugging 1D linear simulations

src/auxval.F90

@@ -4,6 +4,7 @@ subroutine auxval
   USE basic
   USE grid
   USE array
+  USE model
   USE fourier, ONLY: init_grid_distr_and_plans, alloc_local_1, alloc_local_2
   use prec_const
   IMPLICIT NONE
@@ -11,14 +12,18 @@ subroutine auxval
   INTEGER :: irows,irowe, irow, icol, i_
   IF (my_id .EQ. 0) WRITE(*,*) '=== Set auxiliary values ==='
 
-  CALL init_grid_distr_and_plans(Nr,Nz)
+  IF (NON_LIN) THEN
+    CALL init_grid_distr_and_plans(Nr,Nz)
+  ELSE
+    CALL init_1Dgrid_distr
+  ENDIF
 
-  CALL set_krgrid ! MPI Distributed dimension
+  CALL set_pgrid
+  CALL set_jgrid
 
+  CALL set_krgrid ! MPI Distributed dimension
   CALL set_kzgrid
 
-  CALL set_pj
-
   CALL memory ! Allocate memory for global arrays
 
   DO i_ = 0,num_procs-1

src/grid_mod.F90

@@ -55,16 +55,28 @@ MODULE grid
   INTEGER, PUBLIC, PROTECTED ::  ips_i,ipe_i, ijs_i,ije_i
 
   ! Public Functions
-  PUBLIC :: set_pj
+  PUBLIC :: init_1Dgrid_distr
+  PUBLIC :: set_pgrid, set_jgrid
   PUBLIC :: set_krgrid, set_kzgrid
   PUBLIC :: grid_readinputs, grid_outputinputs
   PUBLIC :: bare, bari
 CONTAINS
 
-  SUBROUTINE set_pj
+  SUBROUTINE init_1Dgrid_distr
+
+    write(*,*) Nr
+    local_nkr        = (Nr/2+1)/num_procs
+    write(*,*) local_nkr
+    local_nkr_offset = my_id*local_nkr
+
+    if (my_id .EQ. num_procs-1) local_nkr = (Nr/2+1)-local_nkr_offset
+
+  END SUBROUTINE init_1Dgrid_distr
+
+  SUBROUTINE set_pgrid
     USE prec_const
     IMPLICIT NONE
-    INTEGER :: ip, ij
+    INTEGER :: ip
 
     ips_e = 1; ipe_e = pmaxe/(1+pskip) + 1
     ips_i = 1; ipe_i = pmaxi/(1+pskip) + 1
@@ -73,6 +85,13 @@ CONTAINS
     DO ip = ips_e,ipe_e; parray_e(ip) = (ip-1); END DO
     DO ip = ips_i,ipe_i; parray_i(ip) = (ip-1); END DO
 
+  END SUBROUTINE set_pgrid
+
+  SUBROUTINE set_jgrid
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: ij
+
     ijs_e = 1; ije_e = jmaxe + 1
     ijs_i = 1; ije_i = jmaxi + 1
     ALLOCATE(jarray_e(ijs_e:ije_e))
@@ -81,80 +100,81 @@ CONTAINS
     DO ij = ijs_i,ije_i; jarray_i(ij) = ij-1; END DO
     maxj  = MAX(jmaxi, jmaxe)
 
-  END SUBROUTINE set_pj
-
+  END SUBROUTINE set_jgrid
 
-    SUBROUTINE set_krgrid
-      USE prec_const
-      IMPLICIT NONE
-      INTEGER :: i_
 
-      Nkr = Nr/2+1 ! Defined only on positive kr since fields are real
-      ! Start and END indices of grid
-      ikrs = local_nkr_offset + 1
-      ikre = ikrs + local_nkr - 1
-
-      ! Grid spacings
-      IF (Lr .EQ. 0) THEN
-        deltakr = 1._dp
+  SUBROUTINE set_krgrid
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: i_
+
+    Nkr = Nr/2+1 ! Defined only on positive kr since fields are real
+    ! Start and END indices of grid
+    ikrs = local_nkr_offset + 1
+    ikre = ikrs + local_nkr - 1
+    ALLOCATE(krarray(ikrs:ikre))
+
+    ! Grid spacings
+    IF (Lr .EQ. 0) THEN
+      deltakr = 1._dp
+    ELSE
+      deltakr = 2._dp*PI/Lr
+    ENDIF
+
+    ! Discretized kr positions ordered as dk*(0 1 2 3)
+    DO ikr = ikrs,ikre
+      krarray(ikr) = REAL(ikr-1,dp) * deltakr
+      IF (krarray(ikr) .EQ. 0) ikr_0 = ikr
+    END DO
+
+    ! Orszag 2/3 filter
+    two_third_krmax = 2._dp/3._dp*deltakr*Nkr
+    ALLOCATE(AA_r(ikrs:ikre))
+    DO ikr = ikrs,ikre
+      IF ( (krarray(ikr) .LT. two_third_krmax) ) THEN
+        AA_r(ikr) = 1._dp;
       ELSE
-        deltakr = 2._dp*PI/Lr
+        AA_r(ikr) = 0._dp;
       ENDIF
+    END DO
+  END SUBROUTINE set_krgrid
 
-      ! Discretized kr positions ordered as dk*(0 1 2 3)
-      ALLOCATE(krarray(ikrs:ikre))
-      DO ikr = ikrs,ikre
-        krarray(ikr) = REAL(ikr-1,dp) * deltakr
-        IF (krarray(ikr) .EQ. 0) ikr_0 = ikr
-      END DO
-
-      ! Orszag 2/3 filter
-      two_third_krmax = 2._dp/3._dp*deltakr*Nkr
-      ALLOCATE(AA_r(ikrs:ikre))
-      DO ikr = ikrs,ikre
-        IF ( (krarray(ikr) .LT. two_third_krmax) ) THEN
-          AA_r(ikr) = 1._dp;
-        ELSE
-          AA_r(ikr) = 0._dp;
-        ENDIF
-      END DO
-    END SUBROUTINE set_krgrid
-
-    SUBROUTINE set_kzgrid
-      USE prec_const
-      IMPLICIT NONE
-      INTEGER :: i_
-
-      Nkz = Nz;
-      ! Start and END indices of grid
-      ikzs = 1
-      ikze = Nkz
-
-      ! Grid spacings
-      IF (Lz .EQ. 0) THEN
-        deltakz = 1._dp
-      ELSE
-        deltakz = 2._dp*PI/Lz
-      ENDIF
-      ! Discretized kz positions ordered as dk*(0 1 2 3 -2 -1)
-      ALLOCATE(kzarray(ikzs:ikze))
+  SUBROUTINE set_kzgrid
+    USE prec_const
+    IMPLICIT NONE
+    INTEGER :: i_
+
+    Nkz = Nz;
+    ! Start and END indices of grid
+    ikzs = 1
+    ikze = Nkz
+    ALLOCATE(kzarray(ikzs:ikze))
+
+    ! Grid spacings and discretized kz positions ordered as dk*(0 1 2 3 -2 -1)
+    IF (Lz .EQ. 0) THEN ! 1D linear case
+      deltakz    = 1._dp
+      kzarray(1) = 0
+      ikz_0      = 1
+    ELSE
+      deltakz = 2._dp*PI/Lz
       DO ikz = ikzs,ikze
         kzarray(ikz) = deltakz*(MODULO(ikz-1,Nkz/2)-Nkz/2*FLOOR(2.*real(ikz-1)/real(Nkz)))
         if (ikz .EQ. Nz/2+1)     kzarray(ikz) = -kzarray(ikz)
         IF (kzarray(ikz) .EQ. 0) ikz_0 = ikz
       END DO
-
-      ! Orszag 2/3 filter
-      two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2);
-      ALLOCATE(AA_z(ikzs:ikze))
-      DO ikz = ikzs,ikze
-        IF ( (kzarray(ikz) .GT. -two_third_kzmax) .AND. (kzarray(ikz) .LT. two_third_kzmax) ) THEN
-          AA_z(ikz) = 1._dp;
-        ELSE
-          AA_z(ikz) = 0._dp;
-        ENDIF
-      END DO
-    END SUBROUTINE set_kzgrid
+    ENDIF
+
+    ! Orszag 2/3 filter
+    two_third_kzmax = 2._dp/3._dp*deltakz*(Nkz/2);
+    ALLOCATE(AA_z(ikzs:ikze))
+    DO ikz = ikzs,ikze
+      IF ( (kzarray(ikz) .GT. -two_third_kzmax) .AND. (kzarray(ikz) .LT. two_third_kzmax) ) THEN
+        AA_z(ikz) = 1._dp;
+      ELSE
+        AA_z(ikz) = 0._dp;
+      ENDIF
+    END DO
+  END SUBROUTINE set_kzgrid
 
   SUBROUTINE grid_readinputs
     ! Read the input parameters

src/stepon.F90

@@ -61,7 +61,7 @@ SUBROUTINE stepon
         ! Execution time start
         CALL cpu_time(t0_checkfield)
 
-        IF ( ikrs .EQ. 1 ) CALL enforce_symetry() ! Enforcing symmetry on kr = 0
+        IF ( (ikrs .EQ. 1) .AND. (NON_LIN) ) CALL enforce_symetry() ! Enforcing symmetry on kr = 0
 
         mlend=.FALSE.
         IF(.NOT.nlend) THEN