From 406dee2e4aa83075d076237916204818f8b873ab Mon Sep 17 00:00:00 2001
From: oawile <oawile@7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9>
Date: Fri, 5 Mar 2010 18:08:48 +0000
Subject: [PATCH] - merged Jens' code to fdsolver_* routines - updated some USE
 statements to the new mapping module

git-svn-id: https://ppm.inf.ethz.ch/svn/ppmnumerics/branches/ngtopo/libppmnumerics@534 7c7fe9aa-52eb-4d9e-b0a8-ba7d787348e9
---
 src/ppm_comp_pp_ring.f       |  2 +-
 src/ppm_fdsolver_fft_bd_2d.f | 12 ++++++++----
 src/ppm_fdsolver_fft_bd_3d.f | 15 +++++++++------
 src/ppm_fdsolver_fft_fd_2d.f | 11 +++++++----
 src/ppm_fdsolver_fft_fd_3d.f | 22 +++++++++++++---------
 src/ppm_fdsolver_finalize.f  | 14 +++++++++++---
 src/ppm_fdsolver_map_2d.f    |  2 +-
 src/ppm_fdsolver_map_3d.f    |  2 +-
 src/ppm_fdsolver_solve_2d.f  | 12 +++++++++---
 src/ppm_fdsolver_solve_3d.f  |  1 -
 src/ppm_hamjac_ext_step_3d.f | 12 +++++++++++-
 11 files changed, 71 insertions(+), 34 deletions(-)

diff --git a/src/ppm_comp_pp_ring.f b/src/ppm_comp_pp_ring.f
index fcf8017..be54ec2 100644
--- a/src/ppm_comp_pp_ring.f
+++ b/src/ppm_comp_pp_ring.f
@@ -68,7 +68,7 @@
       USE ppm_module_substop
       USE ppm_module_error
       USE ppm_module_alloc
-      USE ppm_module_map_part
+      USE ppm_module_map
       USE ppm_module_map_part_ring_shift
       USE ppm_module_comp_pp_doring
       IMPLICIT NONE
diff --git a/src/ppm_fdsolver_fft_bd_2d.f b/src/ppm_fdsolver_fft_bd_2d.f
index c4c8680..f7b3507 100644
--- a/src/ppm_fdsolver_fft_bd_2d.f
+++ b/src/ppm_fdsolver_fft_bd_2d.f
@@ -58,7 +58,10 @@
       !-------------------------------------------------------------------------
       USE ppm_module_data
       USE ppm_module_data_fieldsolver
+      USE ppm_module_substart
+      USE ppm_module_substop
       USE ppm_module_error
+      USE ppm_module_alloc
       IMPLICIT NONE
 #if   __KIND == __SINGLE_PRECISION | __KIND ==__SINGLE_PRECISION_COMPLEX 
       INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -103,8 +106,9 @@
 #ifdef __MATHKEISAN
       ! MATHKEISAN variables for MathKeisan FFTs
       INTEGER                                 :: isign_fft,isys
-#if   __KIND == __SINGLE_PRECISION        | __KIND == __DOUBLE_PRECISION
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
+
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       ! parameters for cfft (default=1)
       INTEGER                                 :: incx, incy
 #endif
@@ -169,8 +173,8 @@
       !-------------------------------------------------------------------------
 #if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
       Nx_out = (Nx_in-1)*2
-#elif __KIND == __SINGLE_PRECISION_COMPLEX |__KIND == __DOUBLE_PRECISION_COMPLEX
-      Nx_out = Nx_in -1
+#el#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+   Nx_out = Nx_in -1
 #endif
       Ny_out=Ny_in
       lda(1)=Nx_out+1
diff --git a/src/ppm_fdsolver_fft_bd_3d.f b/src/ppm_fdsolver_fft_bd_3d.f
index 1e6ebcf..ef387ea 100644
--- a/src/ppm_fdsolver_fft_bd_3d.f
+++ b/src/ppm_fdsolver_fft_bd_3d.f
@@ -73,7 +73,10 @@
       !-------------------------------------------------------------------------
       USE ppm_module_data
       USE ppm_module_data_fieldsolver
+      USE ppm_module_substart
+      USE ppm_module_substop
       USE ppm_module_error
+      USE ppm_module_alloc
       IMPLICIT NONE
 #if   __KIND == __SINGLE_PRECISION | __KIND ==__SINGLE_PRECISION_COMPLEX 
       INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -98,9 +101,9 @@
       ! size of data
       INTEGER, DIMENSION(:)         , INTENT(INOUT) :: lda
       ! output data, inverse fast fourier transformed
-#if   __KIND == __SINGLE_PRECISION        | __KIND == __DOUBLE_PRECISION 
+#if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION 
       REAL(MK), DIMENSION(:,:,:)    , POINTER       :: data_out
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       COMPLEX(MK), DIMENSION(:,:,:) , POINTER       :: data_out 
 #elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND ==__DOUBLE_PRECISION_COMPLEX_Z
       COMPLEX(MK), DIMENSION(:,:,:) , POINTER       :: data_out 
@@ -128,7 +131,7 @@
       INTEGER                                 :: isign_fft,isys
 #if __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       INTEGER                                 :: incx, incy
-#elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND ==__DOUBLE_PRECISION_COMPLEX_Z
+#elif __KIND == __SINGLE_PRECISION_COMPLEX_Z | __KIND == __DOUBLE_PRECISION_COMPLEX_Z
       INTEGER                                 :: incx, incy
 #endif
       ! scale of the transformation
@@ -200,9 +203,9 @@
       !-------------------------------------------------------------------------
 #if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
       Nx_out = (Nx_in-1)*2
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       Nx_out = Nx_in-1       
-#elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND ==__DOUBLE_PRECISION_COMPLEX_Z
+#elif __KIND == __SINGLE_PRECISION_COMPLEX_Z | __KIND == __DOUBLE_PRECISION_COMPLEX_Z
       Nx_out = Nx_in-1       
 #endif
       Ny_out=Ny_in
@@ -240,7 +243,7 @@
       lda_work    = 4*Nx_out
 #elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND== __DOUBLE_PRECISION_COMPLEX
       lda_work(1) = 6*Nx_out
-#elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND== __DOUBLE_PRECISION_COMPLEX_Z
+#elif __KIND == __SINGLE_PRECISION_COMPLEX_Z | __KIND== __DOUBLE_PRECISION_COMPLEX_Z
       lda_work(1) = 6*Nx_out
 #endif
       CALL ppm_alloc(work,lda_work,ppm_param_alloc_fit,info)
diff --git a/src/ppm_fdsolver_fft_fd_2d.f b/src/ppm_fdsolver_fft_fd_2d.f
index 4969725..29e21c4 100644
--- a/src/ppm_fdsolver_fft_fd_2d.f
+++ b/src/ppm_fdsolver_fft_fd_2d.f
@@ -58,7 +58,10 @@
       !-------------------------------------------------------------------------
       USE ppm_module_data
       USE ppm_module_data_fieldsolver
+      USE ppm_module_substart
+      USE ppm_module_substop
       USE ppm_module_error
+      USE ppm_module_alloc
       IMPLICIT NONE
 #if   __KIND == __SINGLE_PRECISION | __KIND ==__SINGLE_PRECISION_COMPLEX 
       INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -77,8 +80,8 @@
       ! input data
 #if   __KIND == __SINGLE_PRECISION | __KIND == __DOUBLE_PRECISION
       REAL(MK), DIMENSION(:,:)      , INTENT(IN   ) :: data_in
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
-      COMPLEX(MK), DIMENSION(:,:)   , INTENT(IN   ) :: data_in
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
+      COMPLEX(MK), DIMENSION(:,:)       , INTENT(IN   ) :: data_in
 #endif
       ! size of the array
       INTEGER, DIMENSION(:)         , INTENT(INOUT) :: lda
@@ -168,13 +171,13 @@
       !-------------------------------------------------------------------------
 #if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
       Nx_out = Nx_in/2 + 1
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       Nx_out = Nx_in
 #endif
       Ny_out=Ny_in
 #if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
       lda(1) = Nx_out
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       lda(1) = Nx_out+1      ! to fit ppm-convention
 #endif
       lda(2)=Ny_out
diff --git a/src/ppm_fdsolver_fft_fd_3d.f b/src/ppm_fdsolver_fft_fd_3d.f
index 25d7ad6..847e68b 100644
--- a/src/ppm_fdsolver_fft_fd_3d.f
+++ b/src/ppm_fdsolver_fft_fd_3d.f
@@ -76,7 +76,10 @@
       !-------------------------------------------------------------------------
       USE ppm_module_data
       USE ppm_module_data_fieldsolver
+      USE ppm_module_substart
+      USE ppm_module_substop
       USE ppm_module_error
+      USE ppm_module_alloc
       IMPLICIT NONE
 #if   __KIND == __SINGLE_PRECISION | __KIND ==__SINGLE_PRECISION_COMPLEX 
       INTEGER, PARAMETER :: MK = ppm_kind_single
@@ -99,9 +102,9 @@
       ! input data
 #if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
       REAL(MK), DIMENSION(:,:,:)    , INTENT(IN   ) :: data_in
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       COMPLEX(MK), DIMENSION(:,:,:) , INTENT(IN   ) :: data_in
-#elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND ==__DOUBLE_PRECISION_COMPLEX_Z
+#elif __KIND == __SINGLE_PRECISION_COMPLEX_Z|__KIND == __DOUBLE_PRECISION_COMPLEX_Z
       COMPLEX(MK), DIMENSION(:,:,:) , INTENT(IN   ) :: data_in
 #endif
       ! size of array
@@ -131,7 +134,7 @@
       INTEGER                                 :: isign_fft,isys
 #if __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       INTEGER                                 :: incx, incy
-#elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND==__DOUBLE_PRECISION_COMPLEX_Z
+#elif __KIND == __SINGLE_PRECISION_COMPLEX_Z | __KIND == __DOUBLE_PRECISION_COMPLEX_Z
       INTEGER                                 :: incx, incy
 #endif
       !scale of the transformation
@@ -202,20 +205,20 @@
       !-------------------------------------------------------------------------
       !  Allocate result array
       !-------------------------------------------------------------------------
-#if   __KIND == __SINGLE_PRECISION        | __KIND == __DOUBLE_PRECISION
+#if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
       Nx_out = Nx_in/2 + 1
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       Nx_out = Nx_in     
-#elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND ==__DOUBLE_PRECISION_COMPLEX_Z
+#elif __KIND == __SINGLE_PRECISION_COMPLEX_Z | __KIND == __DOUBLE_PRECISION_COMPLEX_Z
       Nx_out = Nx_in     
 #endif
       Ny_out = Ny_in
       Nz_out = Nz_in
 #if   __KIND == __SINGLE_PRECISION         | __KIND == __DOUBLE_PRECISION
       lda(1) = Nx_out
-#elif __KIND == __SINGLE_PRECISION_COMPLEX| __KIND == __DOUBLE_PRECISION_COMPLEX
+#elif __KIND == __SINGLE_PRECISION_COMPLEX | __KIND == __DOUBLE_PRECISION_COMPLEX
       lda(1) = Nx_out+1      ! to fit ppm-convention
-#elif __KIND==__SINGLE_PRECISION_COMPLEX_Z|__KIND ==__DOUBLE_PRECISION_COMPLEX_Z
+#elif __KIND == __SINGLE_PRECISION_COMPLEX_Z| __KIND == __DOUBLE_PRECISION_COMPLEX_Z
       lda(1) = Nx_out+1      ! to fit ppm-convention
 #endif
       lda(2) = Ny_out
@@ -338,7 +341,8 @@
          ENDDO
       ENDDO     
 #endif
-#if __KIND ==__SINGLE_PRECISION_COMPLEX_Z| __KIND ==__DOUBLE_PRECISION_COMPLEX_Z
+
+#if __KIND == __SINGLE_PRECISION_COMPLEX_Z | __KIND == __DOUBLE_PRECISION_COMPLEX_Z
       !-------------------------------------------------------------------------
       !  Copy margin to conform with PPM convention
       !-------------------------------------------------------------------------
diff --git a/src/ppm_fdsolver_finalize.f b/src/ppm_fdsolver_finalize.f
index 15ea205..9e6e3eb 100644
--- a/src/ppm_fdsolver_finalize.f
+++ b/src/ppm_fdsolver_finalize.f
@@ -54,8 +54,16 @@
       !-------------------------------------------------------------------------
       !  Modules
       !-------------------------------------------------------------------------
-	  USE ppm_module_data_fieldsolver
-      USE ppm_module_data
+
+       USE ppm_module_substart
+       USE ppm_module_substop
+       USE ppm_module_write
+       USE ppm_module_error
+       USE ppm_module_alloc
+       USE ppm_module_data_fieldsolver
+
+
+
       IMPLICIT NONE
       !-------------------------------------------------------------------------
       !  FFTW include
@@ -70,7 +78,7 @@
       !-------------------------------------------------------------------------
       !  Local Variables
       !-------------------------------------------------------------------------
-      REAL(ppm_kind_double)                        :: t0
+      REAL(ppm_kind_double)                                      :: t0
       INTEGER, DIMENSION(3)                        :: lda
       !-------------------------------------------------------------------------
       !  Initialise
diff --git a/src/ppm_fdsolver_map_2d.f b/src/ppm_fdsolver_map_2d.f
index fc856fc..7634052 100644
--- a/src/ppm_fdsolver_map_2d.f
+++ b/src/ppm_fdsolver_map_2d.f
@@ -83,7 +83,7 @@
       !  Modules
       !-------------------------------------------------------------------------
       USE ppm_module_data
-      USE ppm_module_map_field
+      USE ppm_module_map
       USE ppm_module_check_topoid
       USE ppm_module_check_meshid
       USE ppm_module_error
diff --git a/src/ppm_fdsolver_map_3d.f b/src/ppm_fdsolver_map_3d.f
index 8510835..593d7fd 100644
--- a/src/ppm_fdsolver_map_3d.f
+++ b/src/ppm_fdsolver_map_3d.f
@@ -94,7 +94,7 @@
       !  Modules
       !-------------------------------------------------------------------------
       USE ppm_module_data
-      USE ppm_module_map_field
+      USE ppm_module_map
       USE ppm_module_check_topoid
       USE ppm_module_check_meshid
       USE ppm_module_error
diff --git a/src/ppm_fdsolver_solve_2d.f b/src/ppm_fdsolver_solve_2d.f
index be04257..3c1fde2 100644
--- a/src/ppm_fdsolver_solve_2d.f
+++ b/src/ppm_fdsolver_solve_2d.f
@@ -183,7 +183,11 @@
       !-------------------------------------------------------------------------
       USE ppm_module_data
       USE ppm_module_data_mesh
-      USE ppm_module_typedef
+      USE ppm_module_substart
+      USE ppm_module_substop
+      USE ppm_module_write
+      USE ppm_module_error
+      USE ppm_module_alloc
       USE ppm_module_mktopo
       USE ppm_module_mesh_define
       USE ppm_module_fdsolver_map
@@ -274,8 +278,7 @@
       !  Initialise
       !-------------------------------------------------------------------------
       CALL substart('ppm_fdsolver_solve_2d',t0,info)
-      f_topo => ppm_topo(field_topoid)%t
-      f_mesh => f_topo%mesh(mesh_id_data)
+
       !-------------------------------------------------------------------------
       ! Check arguments
       !-------------------------------------------------------------------------
@@ -293,6 +296,9 @@
               GOTO 9999
           ENDIF
       ENDIF
+      
+      f_topo => ppm_topo(field_topoid)%t
+      f_mesh => f_topo%mesh(mesh_id_data)
       !-------------------------------------------------------------------------
       !  Check if FFTW-Library is available of if NEC Library
       !-------------------------------------------------------------------------
diff --git a/src/ppm_fdsolver_solve_3d.f b/src/ppm_fdsolver_solve_3d.f
index 9345558..49bd31c 100644
--- a/src/ppm_fdsolver_solve_3d.f
+++ b/src/ppm_fdsolver_solve_3d.f
@@ -232,7 +232,6 @@
       !-------------------------------------------------------------------------
       USE ppm_module_data
       USE ppm_module_data_mesh
-      USE ppm_module_typedef
       USE ppm_module_mktopo
       USE ppm_module_mesh_define
       USE ppm_module_fdsolver_map
diff --git a/src/ppm_hamjac_ext_step_3d.f b/src/ppm_hamjac_ext_step_3d.f
index 077550d..d83fc5f 100644
--- a/src/ppm_hamjac_ext_step_3d.f
+++ b/src/ppm_hamjac_ext_step_3d.f
@@ -147,6 +147,15 @@
            DO k=1,ndata(3,isubl)
               DO j=1,ndata(2,isubl)
                  DO i=1,ndata(1,isubl)
+                    IF(ABS(phi(i,j,k,isub)).GT.7.0_mk*dx(1)) CYCLE
+                    ! TODO replace the hardcoded 7 by an argument
+                    IF(ABS(phi(i,j,k,isub)).LT.dx(1)) THEN
+#if   __MODE == __SCA
+                       tpsi(i,j,k,isub) = psi(i,j,k,isub)
+#elif __MODE == __VEC
+                       tpsi(1:lda,i,j,k,isub) = psi(1:lda,i,j,k,isub)
+#endif
+                    ELSE
                     phimid(1) = phi(i+1,j,k,isub)-phi(i-1,j,k,isub)
                     phimid(2) = phi(i,j+1,k,isub)-phi(i,j-1,k,isub)
                     phimid(3) = phi(i,j,k+1,isub)-phi(i,j,k-1,isub)
@@ -186,7 +195,8 @@
                     tpsi(1:lda,i,j,k,isub) = psi(1:lda,i,j,k,isub) &
                          & - wenotau * dphi_dt(1:lda)
                     rms = MAX(rms,SUM(ABS(dphi_dt)))
-#endif                    
+#endif
+                 ENDIF
                  END DO
               END DO
            END DO
-- 
GitLab