Changeset 1056 for palm

Timestamp:

Nov 16, 2012 3:28:04 PM (12 years ago)

Author:

raasch

Message:

bugfix for mg-solver

Location:

palm/trunk/SOURCE

Files:

• init_pegrid.f90 (modified) (2 diffs)
• poismg.f90 (modified) (11 diffs)

palm/trunk/SOURCE/init_pegrid.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Indices for arrays n.._mg start from zero due to definition of arrays f2 and
+! p2 as automatic arrays in recursive subroutine next_mg_level
 !
 ! Former revisions:
@@ -860 +861 @@
        ENDIF
 
-       ALLOCATE( grid_level_count(maximum_grid_level),                   &
-                 nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
-                 nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
-                 nzt_mg(maximum_grid_level) )
+       ALLOCATE( grid_level_count(maximum_grid_level),                       &
+                 nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
+                 nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
+                 nzt_mg(0:maximum_grid_level) )
 
        grid_level_count = 0
+!
+!--    Index zero required as dummy due to definition of arrays f2 and p2 in
+!--    recursive subroutine next_mg_level
+       nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
 
        nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt

palm/trunk/SOURCE/poismg.f90

@@ -23 +23 @@
 ! Current revisions:
 ! -----------------
-!
+! Bugfix: all ghost points have to be used for allocating p3
+! arrays p2, f2, and f2_l changed from allocatable to automatic
 !
 ! Former revisions:
@@ -125 +126 @@
                     mg_switch_to_pe0_level)+1) )
     ELSE
-       ALLOCATE ( p3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
-    ENDIF
+       ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+    ENDIF
+
+    p3 = 0.0
  
 !
@@ -618 +621 @@
                     nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                &
                     nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg, p_mg
-
 
     l = grid_level
@@ -1002 +1004 @@
        ENDDO
     ENDDO
-!$OMP END PARALLEL 
+!$OMP END PARALLEL
 
 !
@@ -1027 +1029 @@
     REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
                     nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
-                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2
+                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2, f2_l
 
     REAL, DIMENSION(nzb:mg_loc_ind(5,myid)+1,                            &
@@ -1033 +1035 @@
                     mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub
 
-    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  f2_l
-
-    ALLOCATE( f2_l(nzb:nzt_mg(grid_level)+1,                            &
-                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
-                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
 
 #if defined( __parallel )
@@ -1073 +1070 @@
                         nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
 
-    DEALLOCATE( f2_l )
-
     CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
 #endif
@@ -1150 +1145 @@
                  nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg, p_mg, p3, r
 
-    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  f2, f2_sub, p2, p2_sub
-
+    REAL, DIMENSION(nzb:nzt_mg(grid_level-1)+1,                                &
+                    nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,             &
+                    nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2, p2
+
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub, p2_sub
 
 !
@@ -1190 +1188 @@
        nyn = nyn_mg(grid_level)
        nzt = nzt_mg(grid_level)
-
-       ALLOCATE( f2(nzb:nzt_mg(grid_level)+1,                    &
-                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
-                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1),  &
-                 p2(nzb:nzt_mg(grid_level)+1,                    &
-                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
-                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
 
        IF ( grid_level == mg_switch_to_pe0_level )  THEN
@@ -1281 +1272 @@
 
        ELSE
+
           CALL restrict( f2, r )
 
@@ -1392 +1384 @@
 
 !
-!--    Temporary arrays for the actual grid are not needed any more
-       DEALLOCATE( p2, f2 )
-
-!
 !--    Computation of the new pressure correction. Therefore, 
 !--    values from prior grids are added up automatically stage by stage.