Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Mar 29, 2011 11:39:40 AM (13 years ago)
Author:
raasch
Message:

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:

Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:

localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/pres.f90

    r695 r707  
    44! Current revisions:
    55! -----------------
    6 !
     6! Calculation of weighted average of p is now handled in the same way
     7! regardless of the number of ghost layers (advection scheme),
     8! multigrid and sor method are using p_loc for iterative advancements of
     9! pressure,
     10! localsum calculation modified for proper OpenMP reduction,
     11! bc_lr/ns replaced by bc_lr/ns_cyc
    712!
    813! Former revisions:
    … …  
    99104
    100105!
    101 !-- Multigrid method expects 1 additional grid point for the arrays
    102 !-- d, tend and p
     106!-- Multigrid method expects array d to have one ghost layer.
     107!--
    103108    IF ( psolver == 'multigrid' )  THEN
    104109     
    105110       DEALLOCATE( d )
    106111       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
     112
     113!
     114!--    Since p is later used to hold the weighted average of the substeps, it
     115!--    cannot be used in the iterative solver. Therefore, its initial value is
     116!--    stored on p_loc, which is then iteratively advanced in every substep.
     117       IF ( intermediate_timestep_count == 1 )  THEN
     118          DO  i = nxl-1, nxr+1
     119             DO  j = nys-1, nyn+1
     120                DO  k = nzb, nzt+1
     121                   p_loc(k,j,i) = p(k,j,i)
     122                ENDDO
     123             ENDDO
     124          ENDDO
     125       ENDIF
    107126       
    108        IF ( ws_scheme_mom  .OR. ws_scheme_sca )  THEN
    109        
    110           DEALLOCATE( tend )
    111           ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    112           DEALLOCATE( p )
    113           ALLOCATE( p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    114          
    115        ENDIF
    116        
     127    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count == 1 )  THEN
     128
     129!
     130!--    Since p is later used to hold the weighted average of the substeps, it
     131!--    cannot be used in the iterative solver. Therefore, its initial value is
     132!--    stored on p_loc, which is then iteratively advanced in every substep.
     133       p_loc = p
     134
    117135    ENDIF
    118136
    … …  
    297315    ENDDO
    298316
    299     localsum = ( localsum + threadsum ) * dt_3d * &
    300                weight_pres(intermediate_timestep_count)
     317    localsum = localsum + threadsum * dt_3d * &
     318                          weight_pres(intermediate_timestep_count)
    301319
    302320    !$OMP END PARALLEL
    … …  
    362380       ENDDO
    363381    ENDDO
    364     localsum = ( localsum + threadsum ) * dt_3d                              &
    365                * weight_pres(intermediate_timestep_count)
     382    localsum = localsum + threadsum * dt_3d * &
     383                         weight_pres(intermediate_timestep_count)
    366384    !$OMP END PARALLEL
    367385#endif
    … …  
    371389!-- of the total domain
    372390    sums_divold_l(0:statistic_regions) = localsum
    373 
    374 !
    375 !-- Determine absolute minimum/maximum (only for test cases, therefore as
    376 !-- comment line)
    377 !    CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &
    378 !                        divmax_ijk )
    379391
    380392    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
    … …  
    496508!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
    497509!--    scheme
    498        CALL sor( d, ddzu_pres, ddzw, p )
    499        tend = p
     510       CALL sor( d, ddzu_pres, ddzw, p_loc )
     511       tend = p_loc
    500512
    501513    ELSEIF ( psolver == 'multigrid' )  THEN
    … …  
    506518!--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
    507519!--    ( nbgp ghost points ).
    508        exchange_mg = .TRUE.
    509520       CALL poismg( tend )
    510        exchange_mg = .FALSE.
     521
    511522!
    512523!--    Restore perturbation pressure on tend because this array is used
    513524!--    further below to correct the velocity fields
    514 
    515        tend = p
    516        IF( ws_scheme_mom .OR. ws_scheme_sca )  THEN
    517 !       
    518 !--       Allocate p to its normal size and restore pressure.         
    519           DEALLOCATE( p )
    520           ALLOCATE( p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    521          
    522        ENDIF
    523 
    524     ENDIF
    525 
    526 !
    527 !-- Store perturbation pressure on array p, used in the momentum equations
    528     IF ( psolver(1:7) == 'poisfft' )  THEN
    529 
    530        IF ( intermediate_timestep_count == 1 )  THEN
    531           !$OMP PARALLEL PRIVATE (i,j,k)
    532           !$OMP DO
    533           DO  i = nxlg, nxrg
    534              DO  j = nysg, nyng
    535                 DO  k = nzb, nzt+1
    536                    p(k,j,i) = tend(k,j,i)                                     &
    537                             * weight_substep(intermediate_timestep_count)
    538                 ENDDO
    539              ENDDO
    540           ENDDO
    541           !$OMP END PARALLEL
    542          
    543        ELSE
    544           !$OMP PARALLEL PRIVATE (i,j,k)
    545           !$OMP DO
    546           DO  i = nxlg, nxrg
    547              DO  j = nysg, nyng
    548                 DO  k = nzb, nzt+1
    549                    p(k,j,i) = p(k,j,i) + tend(k,j,i)                          &
    550                             * weight_substep(intermediate_timestep_count)
    551                 ENDDO
    552              ENDDO
    553           ENDDO
    554           !$OMP END PARALLEL
    555          
    556        ENDIF
     525       DO  i = nxl-1, nxr+1
     526          DO  j = nys-1, nyn+1
     527             DO  k = nzb, nzt+1
     528                tend(k,j,i) = p_loc(k,j,i)
     529             ENDDO
     530          ENDDO
     531       ENDDO
     532
     533    ENDIF
     534
     535!
     536!-- Store perturbation pressure on array p, used for pressure data output.
     537!-- Ghost layers are added in the output routines (except sor-method: see below)
     538    IF ( intermediate_timestep_count == 1 )  THEN
     539       !$OMP PARALLEL PRIVATE (i,j,k)
     540       !$OMP DO
     541       DO  i = nxl-1, nxr+1
     542          DO  j = nys-1, nyn+1
     543             DO  k = nzb, nzt+1
     544                p(k,j,i) = tend(k,j,i) * &
     545                           weight_substep(intermediate_timestep_count)
     546             ENDDO
     547          ENDDO
     548       ENDDO
     549       !$OMP END PARALLEL
     550
     551    ELSE
     552       !$OMP PARALLEL PRIVATE (i,j,k)
     553       !$OMP DO
     554       DO  i = nxl-1, nxr+1
     555          DO  j = nys-1, nyn+1
     556             DO  k = nzb, nzt+1
     557                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
     558                           weight_substep(intermediate_timestep_count)
     559             ENDDO
     560          ENDDO
     561       ENDDO
     562       !$OMP END PARALLEL
     563
     564    ENDIF
    557565       
    558     ENDIF
     566!
     567!-- SOR-method needs ghost layers for the next timestep
     568    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
    559569
    560570!
    561571!-- Correction of the provisional velocities with the current perturbation
    562572!-- pressure just computed
    563     IF ( conserve_volume_flow  .AND. &
    564          ( bc_lr == 'cyclic'  .OR.  bc_ns == 'cyclic' ) )  THEN
     573    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .OR.  bc_ns_cyc ) )  THEN
    565574       volume_flow_l(1) = 0.0
    566575       volume_flow_l(2) = 0.0
    567576    ENDIF
     577
    568578    !$OMP PARALLEL PRIVATE (i,j,k)
    569579    !$OMP DO
    … …  
    571581       DO  j = nys, nyn
    572582          DO  k = nzb_w_inner(j,i)+1, nzt
    573              w(k,j,i) = w(k,j,i) - dt_3d *                                    &
    574                         ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)           &
    575                         * weight_pres(intermediate_timestep_count)
     583             w(k,j,i) = w(k,j,i) - dt_3d *                                 &
     584                        ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
     585                        weight_pres(intermediate_timestep_count)
    576586          ENDDO
    577587          DO  k = nzb_u_inner(j,i)+1, nzt
    578588             u(k,j,i) = u(k,j,i) - dt_3d *                                 &
    579                         ( tend(k,j,i) - tend(k,j,i-1) ) * ddx              &
    580                         * weight_pres(intermediate_timestep_count)
     589                        ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
     590                        weight_pres(intermediate_timestep_count)
    581591          ENDDO
    582592          DO  k = nzb_v_inner(j,i)+1, nzt
    583593             v(k,j,i) = v(k,j,i) - dt_3d *                                 &
    584                         ( tend(k,j,i) - tend(k,j-1,i) ) * ddy              &
    585                         * weight_pres(intermediate_timestep_count)
     594                        ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
     595                        weight_pres(intermediate_timestep_count)
    586596          ENDDO                                                         
    587597!
    588598!--       Sum up the volume flow through the right and north boundary
    589           IF ( conserve_volume_flow  .AND.  bc_lr == 'cyclic'  .AND. &
    590                bc_ns == 'cyclic' .AND. i == nx )  THEN
     599          IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND. &
     600               i == nx )  THEN
    591601             !$OMP CRITICAL
    592602             DO  k = nzb_2d(j,i) + 1, nzt
    … …  
    595605             !$OMP END CRITICAL
    596606          ENDIF
    597           IF ( conserve_volume_flow  .AND.  bc_ns == 'cyclic'  .AND. &
    598                bc_lr == 'cyclic' .AND. j == ny )  THEN
     607          IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND. &
     608               j == ny )  THEN
    599609             !$OMP CRITICAL
    600610             DO  k = nzb_2d(j,i) + 1, nzt
    … …  
    608618    !$OMP END PARALLEL
    609619   
    610     IF ( psolver == 'multigrid' .OR. psolver == 'sor' )  THEN       
    611        IF ( intermediate_timestep_count == 1 .OR. simulated_time == 0)  THEN
    612           !$OMP PARALLEL PRIVATE (i,j,k)
    613           !$OMP DO
    614           DO  i = nxl, nxr
    615              DO  j = nys, nyn
    616                 DO  k = nzb, nzt+1
    617                    p_sub(k,j,i) = tend(k,j,i)                                 &
    618                                 * weight_substep(intermediate_timestep_count)
    619                 ENDDO
    620              ENDDO
    621           ENDDO
    622           !$OMP END PARALLEL
    623        ELSE
    624           !$OMP PARALLEL PRIVATE (i,j,k)
    625           !$OMP DO
    626           DO  i = nxl, nxr
    627              DO  j = nys, nyn
    628                 DO  k = nzb, nzt+1
    629                    p_sub(k,j,i) = p_sub(k,j,i) + tend(k,j,i)                  &
    630                                 * weight_substep(intermediate_timestep_count)
    631                 ENDDO
    632              ENDDO
    633           ENDDO
    634           !$OMP END PARALLEL
    635        ENDIF
    636          
    637        IF ( intermediate_timestep_count == intermediate_timestep_count_max )  &
    638           THEN
    639           !$OMP PARALLEL PRIVATE (i,j,k)
    640           !$OMP DO
    641           DO  i = nxl, nxr
    642              DO  j = nys, nyn
    643                 DO  k = nzb, nzt+1
    644                    p(k,j,i) = p_sub(k,j,i)
    645                 ENDDO
    646              ENDDO
    647           ENDDO
    648           !$OMP END PARALLEL
    649        ENDIF
    650     ENDIF
    651  
    652 !
    653 !-- Resize tend to its normal size in case of multigrid and ws-scheme.
    654     IF ( psolver == 'multigrid' .AND. ( ws_scheme_mom        &
    655                                    .OR. ws_scheme_sca ) )  THEN
    656        DEALLOCATE( tend )
    657        ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    658     ENDIF
    659 
    660 
    661620!
    662621!-- Conserve the volume flow
    663     IF ( conserve_volume_flow  .AND. &
    664          ( bc_lr == 'cyclic'  .AND.  bc_ns == 'cyclic' ) )  THEN
     622    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
    665623
    666624#if defined( __parallel )   
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |