Changeset 778 for palm/trunk/SOURCE

Timestamp:

Nov 7, 2011 2:18:25 PM (13 years ago)

Author:

fricke

Message:

Modifications of the multigrid pressure solver

Location:

palm/trunk/SOURCE

Files:

• init_pegrid.f90 (modified) (8 diffs)
• modules.f90 (modified) (3 diffs)
• poismg.f90 (modified) (16 diffs)
• pres.f90 (modified) (2 diffs)

palm/trunk/SOURCE/init_pegrid.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! Calculation of subdomain_size now considers the number of ghost points.
+! Further coarsening on PE0 is now possible for multigrid solver if the 
+! collected field has more grid points than the subdomain of an PE.
 !
 ! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
@@ -114 +116 @@
     IMPLICIT NONE
 
-    INTEGER ::  gathered_size, i, id_inflow_l, id_recycling_l, ind(5), j, k, &
+    INTEGER ::  i, id_inflow_l, id_recycling_l, ind(5), j, k,                &
                 maximum_grid_level_l, mg_switch_to_pe0_level_l, mg_levels_x, &
                 mg_levels_y, mg_levels_z, nnx_y, nnx_z, nny_x, nny_z, nnz_x, &
                 nnz_y, numproc_sqr, nx_total, nxl_l, nxr_l, nyn_l, nys_l,    &
-                nzb_l, nzt_l, omp_get_num_threads, subdomain_size
+                nzb_l, nzt_l, omp_get_num_threads
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::  ind_all, nxlf, nxrf, nynf, nysf
@@ -879 +881 @@
 
           ELSE
-
              mg_switch_to_pe0_level_l = 0
              maximum_grid_level_l = maximum_grid_level
@@ -889 +890 @@
 !--       by user
           IF ( mg_switch_to_pe0_level == 0 )  THEN
-
              IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
                 mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
@@ -922 +922 @@
 
        grid_level_count = 0
+
        nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
 
@@ -952 +953 @@
 !--          The size of this gathered array must not be larger than the
 !--          array tend, which is used in the multigrid scheme as a temporary
-!--          array
-             subdomain_size = ( nxr - nxl + 3 )     * ( nyn - nys + 3 )     * &
-                              ( nzt - nzb + 2 )
+!--          array. Therefore the subdomain size of an PE is calculated and 
+!--          the size of the gathered grid. These values are used in  
+!--          routines pres and poismg
+             subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
+                              ( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
              gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
                               ( nzt_l - nzb + 2 )
 
-             IF ( gathered_size > subdomain_size )  THEN
-                message_string = 'not enough memory for storing ' // &
-                                 'gathered multigrid data on PE0'
-                CALL message( 'init_pegrid', 'PA0236', 1, 2, 0, 6, 0 )
-             ENDIF
 #else
              message_string = 'multigrid gather/scatter impossible ' // &
@@ -981 +979 @@
           nyn_l = nyn_l / 2 
           nzt_l = nzt_l / 2 
+
        ENDDO
 
@@ -987 +986 @@
        maximum_grid_level = 0
 
+    ENDIF
+
+!-- Temporary problem: In the moment the summation to calculate maxerror
+!-- in routine poismg crashes the program if the first coarsement is on PE0.
+!-- Further work required.
+    IF ( maximum_grid_level == mg_switch_to_pe0_level )  THEN
+       message_string = 'At least one coarser grid must be calculated ' // &
+                        'on the subdomain of each PE'
+       CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
     ENDIF
 

palm/trunk/SOURCE/modules.f90

@@ -5 +5 @@
 ! Current revisions:
 ! -----------------
+! +gathered_size, subdomain_size
 !
 ! Former revisions:
@@ -855 +856 @@
                 nxlu, nxr, nxra, nxrg, nx_on_file, nny, ny = 0, ny_a, ny_o,    &
                 nya, nyn, nyna, nyng, nys, nysg, nysv, ny_on_file, nnz, nz = 0,&
-                nza, nzb, nzb_diff, nzt, nzta, nzt_diff
+                nza, nzb, nzb_diff, nzt, nzta, nzt_diff               
 
 
@@ -1250 +1251 @@
 #endif
 
-    INTEGER ::  comm1dx, comm1dy, comm2d, comm_inter, comm_palm, ierr, myidx,  &
-                myidy, ndim = 2, ngp_a, ngp_o, ngp_xy, ngp_y,                  &
+    INTEGER ::  comm1dx, comm1dy, comm2d, comm_inter, comm_palm, gathered_size,&
+                ierr, myidx, myidy, ndim = 2, ngp_a, ngp_o, ngp_xy, ngp_y,     &
                 pleft, pnorth, pright, psouth,                                 &
                 sendrecvcount_xy, sendrecvcount_yz, sendrecvcount_zx,          &
-                sendrecvcount_zyd, sendrecvcount_yxd,                          &
+                sendrecvcount_zyd, sendrecvcount_yxd, subdomain_size,          &
                 type_x, type_x_int, type_xy, type_y, type_y_int
 

palm/trunk/SOURCE/poismg.f90

@@ -7 +7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Allocation of p3 changes when multigrid is used and the collected field on PE0
+! has more grid points than the subdomain of an PE.
 !
 ! Former revisions:
@@ -75 +76 @@
     REAL    ::  maxerror, maximum_mgcycles, residual_norm
 
-    REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r
+    REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r 
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  p3
@@ -81 +82 @@
 
     CALL cpu_log( log_point_s(29), 'poismg', 'start' )
-
 !
 !-- Initialize arrays and variables used in this subroutine
-    ALLOCATE ( p3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
-
+
+!--    If the number of grid points of the gathered grid, which is collected
+!--    on PE0, is larger than the number of grid points of an PE, than array
+!--    p3 will be enlarged. 
+    IF ( gathered_size > subdomain_size )  THEN 
+       ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(         &
+                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,   &
+                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                    &
+                    mg_switch_to_pe0_level)+1) )
+    ELSE
+       ALLOCATE ( p3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+    ENDIF
 !
 !-- Ghost boundaries have to be added to divergence array.
@@ -111 +121 @@
     DO WHILE ( residual_norm > residual_limit  .OR. &
                mgcycles < maximum_mgcycles )
-
+ 
        CALL next_mg_level( d, p_loc, p3, r)
-        
+
 !
 !--    Calculate the residual if the user has not preset the number of
@@ -120 +130 @@
           CALL resid( d, p_loc, r )
           maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
+
 #if defined( __parallel )
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
+             CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
                               comm2d, ierr)
 #else
-          residual_norm = maxerror
+             residual_norm = maxerror
 #endif
           residual_norm = SQRT( residual_norm )
@@ -612 +623 @@
 
           IF ( .NOT. unroll )  THEN
+
              CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
 
@@ -964 +976 @@
     CALL exchange_horiz( p_mg, 1)
 
+
  END SUBROUTINE redblack
 
@@ -1107 +1120 @@
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::  f2, f2_sub, p2, p2_sub
 
+
 !
 !-- Restriction to the coarsest grid
@@ -1115 +1129 @@
 !--    iterations in order to get a more accurate solution.
        ngsrb = 2 * ngsrb
+
        CALL redblack( f_mg, p_mg )
+
        ngsrb = ngsrb / 2
+
 
     ELSEIF ( grid_level /= 1 )  THEN
@@ -1137 +1154 @@
        grid_level = grid_level - 1
        nxl = nxl_mg(grid_level)
+       nys = nys_mg(grid_level)
        nxr = nxr_mg(grid_level)
-       nys = nys_mg(grid_level)
        nyn = nyn_mg(grid_level)
        nzt = nzt_mg(grid_level)
@@ -1150 +1167 @@
 
        IF ( grid_level == mg_switch_to_pe0_level )  THEN
+
 !
 !--       From this level on, calculations are done on PE0 only.
@@ -1156 +1174 @@
 !--       in between (otherwise, the restrict routine would expect
 !--       the gathered array)
+
           nxl_mg_save = nxl_mg(grid_level)
           nxr_mg_save = nxr_mg(grid_level)
@@ -1190 +1209 @@
           nyn = nyn_mg(grid_level)
           nzt = nzt_mg(grid_level)
-
 !
 !--       Gather all arrays from the subdomains on PE0
@@ -1231 +1249 @@
 
        ELSE
-
           CALL restrict( f2, r )
 
@@ -1334 +1351 @@
 
        ELSE
-
           CALL prolong( p2, p3 )
 
@@ -1360 +1376 @@
     ENDIF
 
+
 !
 !-- The following few lines serve the steering of the multigrid scheme 

palm/trunk/SOURCE/pres.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! New allocation of tend when multigrid is used and the collected field on PE0
+! has more grid points than the subdomain of an PE.
 !
 ! Former revisions:
@@ -526 +527 @@
 !--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
 !--    ( nbgp ghost points ).
+
+!--    If the number of grid points of the gathered grid, which is collected
+!--    on PE0, is larger than the number of grid points of an PE, than array
+!--    tend will be enlarged. 
+       IF ( gathered_size > subdomain_size )  THEN
+          DEALLOCATE( tend )
+          ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(          &
+                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
+                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
+                    mg_switch_to_pe0_level)+1) )
+       ENDIF
+
        CALL poismg( tend )
+
+       IF ( gathered_size > subdomain_size )  THEN
+          DEALLOCATE( tend )
+          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       ENDIF
 
 !