Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLO_ModElliptSmoothing.F90
Go to the documentation of this file.
1 ! *********************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4 ! * *
5 ! * Illinois Rocstar LLC *
6 ! * Champaign, IL *
7 ! * www.illinoisrocstar.com *
8 ! * sales@illinoisrocstar.com *
9 ! * *
10 ! * License: See LICENSE file in top level of distribution package or *
11 ! * http://opensource.org/licenses/NCSA *
12 ! *********************************************************************
13 ! *********************************************************************
14 ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15 ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16 ! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17 ! * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18 ! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19 ! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20 ! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21 ! * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22 ! *********************************************************************
23 ! ******************************************************************************
24 !
25 ! Purpose: Suite for grid smoothing routines based on elliptic PDE.
26 !
27 ! Description: None.
28 !
29 ! ******************************************************************************
30 !
31 ! $Id: RFLO_ModElliptSmoothing.F90,v 1.15 2008/12/06 08:44:15 mtcampbe Exp $
32 !
33 ! Copyright: (c) 2004 by the University of Illinois
34 !
35 ! ******************************************************************************
36 
38 
39  USE modglobal, ONLY : t_global
40  USE moddatastruct, ONLY: t_region
41  USE modgrid, ONLY : t_grid
42  USE modbndpatch, ONLY : t_patch
43  USE modparameters
44  USE moddatatypes
45  USE moderror
46  USE modmpi
47 
48  IMPLICIT NONE
49 
50  PRIVATE
51  PUBLIC :: rflo_elliptgridsmoo, &
62 
63 ! private :
64 
65 ! ******************************************************************************
66 ! Declarations and definitions
67 ! ******************************************************************************
68 
69  CHARACTER(CHRLEN) :: RCSIdentString = &
70  '$RCSfile: RFLO_ModElliptSmoothing.F90,v $ $Revision: 1.15 $'
71 
72 ! ******************************************************************************
73 ! Routines
74 ! ******************************************************************************
75 
76  CONTAINS
77 
78 !******************************************************************************
79 !
80 ! Purpose: smooth the distribution of grid points by solving 3D elliptic PDE
81 ! in physical space.
82 !
83 ! Description: none.
84 !
85 ! Input: regions = data of all grid regions.
86 !
87 ! Output: regions%levels%grid%xyz = new grid coordinates
88 ! resid = convergence of the Jacobi iteration.
89 !
90 ! Notes: none.
91 !
92 !******************************************************************************
93 
94 SUBROUTINE rflo_elliptgridsmoo( regions,resid )
95 
103  IMPLICIT NONE
104 
105 #include "Indexing.h"
106 
107 ! ... parameters
108  REAL(RFREAL) :: resid
109  TYPE(t_region), POINTER :: regions(:)
110 
111 ! ... loop variables
112  INTEGER :: ireg, ipatch, ijk, i, j, k
113 
114 ! ... local variables
115  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend, inoff, ijnoff
116  INTEGER :: bctype, iregsrc, ipatchsrc
117 
118  REAL(RFREAL) :: dx, dy, dz
119  REAL(RFREAL), POINTER :: xyz(:,:), xyzold(:,:), xyzorig(:,:)
120 
121  TYPE(t_global), POINTER :: global
122  TYPE(t_grid), POINTER :: grid, gridold, gridsrc
123  TYPE(t_patch), POINTER :: patch, patchsrc
124 
125 !******************************************************************************
126 
127  global => regions(1)%global
128 
129  CALL registerfunction( global,'RFLO_ElliptGridSmoo',&
130  'RFLO_ModElliptSmoothing.F90' )
131 
132 ! smooth grid region-wise -----------------------------------------------------
133 
134  DO ireg=1,global%nRegions
135  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
136  regions(ireg)%active==active .AND. & ! on my processor
137  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
138 
139  CALL rflo_gridphysgrad3d( regions(ireg) )
140 
141 ! --- copy coordinates
142 
143  xyz => regions(ireg)%levels(1)%grid%xyz
144  xyzold => regions(ireg)%levels(1)%grid%xyzOld
145  xyzold = xyz
146 
147 ! --- compute grid control map
148 
149  CALL rflo_elliptmatrixcoeffs3d( regions(ireg) )
150 
151  CALL rflo_elliptgridjac3d( regions(ireg) )
152 ! CALL RFLO_ElliptGridGauss3D( regions(iReg) )
153 ! CALL RFLO_ElliptGridSOR3D( regions(iReg) )
154 
155 ! --- transform grid coordinates to deformations
156 
157  xyz => regions(ireg)%levels(1)%grid%xyz
158  xyzorig => regions(ireg)%levels(1)%gridOld%xyz
159 
160  DO ijk=lbound(xyz,2),ubound(xyz,2)
161  xyz(xcoord,ijk) = xyz(xcoord,ijk) - xyzorig(xcoord,ijk)
162  xyz(ycoord,ijk) = xyz(ycoord,ijk) - xyzorig(ycoord,ijk)
163  xyz(zcoord,ijk) = xyz(zcoord,ijk) - xyzorig(zcoord,ijk)
164  ENDDO
165 
166  DO ijk=lbound(xyz,2),ubound(xyz,2)
167  xyzold(xcoord,ijk) = xyzold(xcoord,ijk) - xyzorig(xcoord,ijk)
168  xyzold(ycoord,ijk) = xyzold(ycoord,ijk) - xyzorig(ycoord,ijk)
169  xyzold(zcoord,ijk) = xyzold(zcoord,ijk) - xyzorig(zcoord,ijk)
170  ENDDO
171 
172 ! --- zero out movements along certain boundaries, comment it out
173 ! in case deformations extended to boundary points
174 
175 ! DO iPatch=1,regions(iReg)%nPatches
176 ! patch => regions(iReg)%levels(1)%patches(iPatch)
177 ! bcType = patch%bcType
178 ! IF ((bcType>=BC_INFLOW .AND. bcType<=BC_INFLOW +BC_RANGE) .OR. &
179 ! (bcType>=BC_OUTFLOW .AND. bcType<=BC_OUTFLOW +BC_RANGE) .OR. &
180 ! (bcType>=BC_SLIPWALL .AND. bcType<=BC_SLIPWALL +BC_RANGE) .OR. &
181 ! (bcType>=BC_NOSLIPWALL .AND. bcType<=BC_NOSLIPWALL+BC_RANGE) .OR. &
182 ! (bcType>=BC_FARFIELD .AND. bcType<=BC_FARFIELD +BC_RANGE) .OR. &
183 ! (bcType>=BC_INJECTION .AND. bcType<=BC_INJECTION +BC_RANGE) .OR. &
184 ! (bcType>=BC_TRA_PERI .AND. bcType<=BC_TRA_PERI +BC_RANGE) .OR. &
185 ! (bcType>=BC_ROT_PERI .AND. bcType<=BC_ROT_PERI +BC_RANGE)) THEN
186 ! CALL RFLO_ElliptGridPatch( regions(iReg),patch )
187 ! ENDIF ! bcType
188 ! ENDDO ! iPatch
189 
190  CALL rflo_getdimensphysnodes( regions(ireg),1,ipnbeg,ipnend, &
191  jpnbeg,jpnend,kpnbeg,kpnend )
192  CALL rflo_getnodeoffset( regions(ireg),1,inoff,ijnoff )
193 
194  resid = 0._rfreal
195  DO k=kpnbeg,kpnend
196  DO j=jpnbeg,jpnend
197  DO i=ipnbeg,ipnend
198  ijk = indijk(i,j,k,inoff,ijnoff)
199  dx = xyz(xcoord,ijk) - xyzold(xcoord,ijk)
200  dy = xyz(ycoord,ijk) - xyzold(ycoord,ijk)
201  dz = xyz(zcoord,ijk) - xyzold(zcoord,ijk)
202 
203  resid = resid + dx*dx + dy*dy +dz*dz
204  ENDDO
205  ENDDO
206  ENDDO
207 
208  ENDIF ! region on this processor and active, grid moving
209  ENDDO ! iReg
210 
211 ! fix interfaces between regions ----------------------------------------------
212 ! copy / send deformations
213 
214  DO ireg=1,global%nRegions
215  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
216  regions(ireg)%active==active .AND. & ! on my processor
217  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
218 
219  grid => regions(ireg)%levels(1)%grid
220  gridold => regions(ireg)%levels(1)%gridOld
221 
222  DO ipatch=1,regions(ireg)%nPatches
223  patch => regions(ireg)%levels(1)%patches(ipatch)
224  bctype = patch%bcType
225  IF ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
226  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range) .OR. &
227  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range)) THEN
228  iregsrc = patch%srcRegion
229  ipatchsrc = patch%srcPatch
230  patchsrc => regions(iregsrc)%levels(1)%patches(ipatchsrc)
231  gridsrc => regions(iregsrc)%levels(1)%grid
232 
233  IF (regions(iregsrc)%procid == global%myProcid) THEN
234  CALL rflo_exchangednodecopy( regions(ireg),regions(iregsrc), &
235  patch,patchsrc,.true., &
236  grid%xyz,gridsrc%xyz )
237  ELSE
238  CALL rflo_exchangednodesend( regions(ireg),regions(iregsrc), &
239  patch,grid%xyz )
240  ENDIF
241  ENDIF ! bcType
242  ENDDO ! iPatch
243 
244  ENDIF ! region on this processor and active, grid moving
245  ENDDO ! iReg
246 
247 ! receive deformations
248 
249  DO ireg=1,global%nRegions
250  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
251  regions(ireg)%active==active .AND. & ! on my processor
252  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
253 
254  grid => regions(ireg)%levels(1)%grid
255  gridold => regions(ireg)%levels(1)%gridOld
256 
257  DO ipatch=1,regions(ireg)%nPatches
258  patch => regions(ireg)%levels(1)%patches(ipatch)
259  bctype = patch%bcType
260  IF ((bctype>=bc_regionconf .AND. bctype<=bc_regionconf+bc_range) .OR. &
261  (bctype>=bc_tra_peri .AND. bctype<=bc_tra_peri +bc_range) .OR. &
262  (bctype>=bc_rot_peri .AND. bctype<=bc_rot_peri +bc_range)) THEN
263  iregsrc = patch%srcRegion
264  ipatchsrc = patch%srcPatch
265  patchsrc => regions(iregsrc)%levels(1)%patches(ipatchsrc)
266  gridsrc => regions(iregsrc)%levels(1)%grid
267 
268  IF (regions(iregsrc)%procid /= global%myProcid) THEN
269  CALL rflo_exchangednoderecv( regions(ireg),regions(iregsrc), &
270  patch,patchsrc,.true.,grid%xyz )
271  ENDIF
272  ENDIF ! bcType
273  ENDDO ! iPatch
274 
275  ENDIF ! region on this processor and active, grid moving
276  ENDDO ! iReg
277 
278 ! clear send requests
279 
280  DO ireg=1,global%nRegions
281  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
282  regions(ireg)%active==active .AND. & ! on my processor
283  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
284  CALL rflo_clearsendrequests( regions,ireg,.true. )
285  ENDIF
286  ENDDO
287 
288 ! update grid, dummy, corner and edge cells -----------------------------------
289 
290  DO ireg=1,global%nRegions
291  IF (regions(ireg)%procid==global%myProcid .AND. & ! region active and
292  regions(ireg)%active==active .AND. & ! on my processor
293  regions(ireg)%mixtInput%moveGrid) THEN ! and moving
294 
295 ! --- change xyz from deformations to coordinates
296 
297  xyz => regions(ireg)%levels(1)%grid%xyz
298  xyzorig => regions(ireg)%levels(1)%gridOld%xyz
299 
300  DO ijk=lbound(xyz,2),ubound(xyz,2)
301  xyz(xcoord,ijk) = xyz(xcoord,ijk) + xyzorig(xcoord,ijk)
302  xyz(ycoord,ijk) = xyz(ycoord,ijk) + xyzorig(ycoord,ijk)
303  xyz(zcoord,ijk) = xyz(zcoord,ijk) + xyzorig(zcoord,ijk)
304  ENDDO
305 
306 ! --- update coarse grids and dummy cells
307 
308  CALL rflo_generatecoarsegrids( regions(ireg) ) ! coarsen finest grid
309  CALL rflo_copygeometrydummy( regions(ireg) ) ! copy to dummy nodes
310  CALL rflo_extrapolategeometry( regions(ireg) ) ! extrapolate
311  ENDIF ! region on this processor and active, grid moving
312  ENDDO ! iReg
313 
314  CALL rflo_exchangegeometry( regions ) ! exchange geometry
315 
316 ! finalize --------------------------------------------------------------------
317 
318  CALL deregisterfunction( global )
319 
320 END SUBROUTINE rflo_elliptgridsmoo
321 
322 
323 !******************************************************************************
324 !
325 ! Purpose: smooth the distribution of grid points regionwise by solving
326 ! 3D elliptic PDE in physical space.
327 !
328 ! Description: none.
329 !
330 ! Input: region = data of current region.
331 !
332 ! Output: regions%levels%grid%xyz = new grid coordinates
333 ! resid = convergence of the Jacobi iteration.
334 !
335 ! Notes: none.
336 !
337 !******************************************************************************
338 
339 SUBROUTINE rflo_elliptgridsmooregion( region,resid )
340 
343 
344  IMPLICIT NONE
345 
346 #include "Indexing.h"
347 
348 ! ... parameters
349  REAL(RFREAL) :: resid
350  TYPE(t_region) :: region
351 
352 ! ... loop variables
353  INTEGER :: ijk, i, j, k
354 
355 ! ... local variables
356  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend, inoff, ijnoff
357 
358  REAL(RFREAL) :: dx, dy, dz
359  REAL(RFREAL), POINTER :: xyz(:,:), xyzold(:,:)
360 
361  TYPE(t_global), POINTER :: global
362 
363 !******************************************************************************
364 
365  global => region%global
366 
367  CALL registerfunction( global,'RFLO_ElliptGridSmooRegion',&
368  'RFLO_ModElliptSmoothing.F90' )
369 
370 ! compute gradients of grid coordinates ---------------------------------------
371 
372  CALL rflo_gridphysgrad3d( region )
373 
374 ! copy coordinates
375 
376  xyz => region%levels(1)%grid%xyz
377  xyzold => region%levels(1)%grid%xyzOld
378  xyzold = xyz
379 
380 ! compute grid control map
381 
382  CALL rflo_elliptmatrixcoeffs3d( region )
383 
384 ! solve for xyz
385 
386  CALL rflo_elliptgridjac3d( region )
387 ! CALL RFLO_ElliptGridGauss3D( region )
388 ! CALL RFLO_ElliptGridSOR3D( region )
389 
390 ! compute residual
391 
392  CALL rflo_getdimensphysnodes( region,1,ipnbeg,ipnend, &
393  jpnbeg,jpnend,kpnbeg,kpnend )
394  CALL rflo_getnodeoffset( region,1,inoff,ijnoff )
395 
396  resid = 0._rfreal
397  DO k=kpnbeg,kpnend
398  DO j=jpnbeg,jpnend
399  DO i=ipnbeg,ipnend
400  ijk = indijk(i,j,k,inoff,ijnoff)
401  dx = xyz(xcoord,ijk) - xyzold(xcoord,ijk)
402  dy = xyz(ycoord,ijk) - xyzold(ycoord,ijk)
403  dz = xyz(zcoord,ijk) - xyzold(zcoord,ijk)
404 
405  resid = resid + dx*dx + dy*dy +dz*dz
406  ENDDO
407  ENDDO
408  ENDDO
409 
410 ! finalize --------------------------------------------------------------------
411 
412  CALL deregisterfunction( global )
413 
414 END SUBROUTINE rflo_elliptgridsmooregion
415 
416 
417 
418 !******************************************************************************
419 !
420 ! Purpose: zero out movements obtained by previous treatments along a boundary.
421 !
422 ! Description: none.
423 !
424 ! Input: region = data of current region, grid movements
425 ! patch = current patch.
426 !
427 ! Output: region%levels%grid%xyz = new grid movements.
428 !
429 ! Notes: none.
430 !
431 !******************************************************************************
432 
433 SUBROUTINE rflo_elliptgridpatch( region,patch )
434 
436 
437  IMPLICIT NONE
438 
439 #include "Indexing.h"
440 
441 ! ... parameters
442  TYPE(t_region) :: region
443  TYPE(t_patch) :: patch
444 
445 ! ... loop variables
446  INTEGER :: i, j, k
447 
448 ! ... local variables
449  INTEGER :: ilev
450  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, inoff, ijnoff, ijknb
451 
452  REAL(RFREAL), POINTER :: xyz(:,:), xyzold(:,:)
453 
454 !******************************************************************************
455 
456  CALL registerfunction( region%global,'RFLO_ElliptGridPatch',&
457  'RFLO_ModElliptSmoothing.F90' )
458 
459 ! get dimensions and pointers
460 
461  ilev = 1
462 
463  CALL rflo_getpatchindicesnodes( region,patch,ilev,ibeg,iend, &
464  jbeg,jend,kbeg,kend )
465  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
466 
467  xyz => region%levels(ilev)%grid%xyz
468  xyzold => region%levels(ilev)%grid%xyzOld
469 
470 ! for grid coordinates: new = old
471 
472 ! DO k=kbeg,kend
473 ! DO j=jbeg,jend
474 ! DO i=ibeg,iend
475 ! ijkNB = IndIJK(i,j,k,iNOff,ijNOff)
476 ! xyz(XCOORD,ijkNB) = xyzOld(XCOORD,ijkNB)
477 ! xyz(YCOORD,ijkNB) = xyzOld(YCOORD,ijkNB)
478 ! xyz(ZCOORD,ijkNB) = xyzOld(ZCOORD,ijkNB)
479 ! ENDDO
480 ! ENDDO
481 ! ENDDO
482 
483 ! for deformations: zero out deformations
484 
485  DO k=kbeg,kend
486  DO j=jbeg,jend
487  DO i=ibeg,iend
488  ijknb = indijk(i,j,k,inoff,ijnoff)
489  xyz(xcoord,ijknb) = 0._rfreal
490  xyz(ycoord,ijknb) = 0._rfreal
491  xyz(zcoord,ijknb) = 0._rfreal
492  ENDDO
493  ENDDO
494  ENDDO
495 
496 ! finalize
497 
498  CALL deregisterfunction( region%global )
499 
500 END SUBROUTINE rflo_elliptgridpatch
501 
502 
503 !******************************************************************************
504 !
505 ! Purpose: compute matrix coefficients for elliptic PDE grid smooting.
506 !
507 ! Description: Equation 4.68 in Handbook of Grid Generation, by
508 ! Joe F. Thompson, Bharat K. Soni, and Nigel P. Weatherill,
509 ! CRC Press 1999, ISBN 0-8493-2687-7.
510 !
511 ! Input: region = grid data of current region
512 !
513 ! Output: grid%aijk, aimjk, ... = matrix coefficients
514 !
515 ! Notes: none.
516 !
517 !******************************************************************************
518 
519 SUBROUTINE rflo_elliptmatrixcoeffs3d( region )
520 
523 
524  IMPLICIT NONE
525 
526 #include "Indexing.h"
527 
528 ! ... parameters
529  TYPE(t_region) :: region
530 
531 ! ... loop variables
532  INTEGER :: i, j, k
533 
534 ! ... local variables
535  INTEGER :: ilev, ijkn, inoff, ijnoff
536  INTEGER :: idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend
537  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend
538 
539  REAL(RFREAL) :: di, dj, dk, rdi2, rdj2, rdk2, r2di, r2dj, r2dk, q1, q2, q3
540  REAL(RFREAL) :: al11, al12, al13, al22, al23, al33
541  REAL(RFREAL) :: au11, au12, au13, au22, au23, au33
542  REAL(RFREAL), POINTER :: stui(:,:), stuj(:,:), stuk(:,:), pmat(:,:,:)
543  REAL(RFREAL), POINTER :: aijk(:), aimjk(:), aipjk(:), aijmk(:), aijpk(:)
544  REAL(RFREAL), POINTER :: aijkm(:), aijkp(:), aimjmk(:), aipjmk(:)
545  REAL(RFREAL), POINTER :: aimjpk(:), aipjpk(:), aimjkm(:), aipjkm(:)
546  REAL(RFREAL), POINTER :: aimjkp(:), aipjkp(:), aijmkm(:), aijpkm(:)
547  REAL(RFREAL), POINTER :: aijmkp(:), aijpkp(:)
548 
549  TYPE(t_global), POINTER :: global
550 
551 !******************************************************************************
552 
553  global => region%global
554 
555  CALL registerfunction( global,'RFLO_ElliptMatrixCoeffs3D',&
556  'RFLO_ModElliptSmoothing.F90' )
557 
558 ! get dimensions and pointers -------------------------------------------------
559 
560  ilev = 1
561  CALL rflo_getdimensdummynodes( region,ilev,idnbeg,idnend, &
562  jdnbeg,jdnend,kdnbeg,kdnend )
563  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend, &
564  jpnbeg,jpnend,kpnbeg,kpnend )
565  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
566 
567  stui => region%levels(ilev)%grid%stui
568  stuj => region%levels(ilev)%grid%stuj
569  stuk => region%levels(ilev)%grid%stuk
570  pmat => region%levels(ilev)%grid%pmat
571 
572  aijk => region%levels(ilev)%grid%aijk
573  aimjk => region%levels(ilev)%grid%aimjk
574  aipjk => region%levels(ilev)%grid%aipjk
575  aijmk => region%levels(ilev)%grid%aijmk
576  aijpk => region%levels(ilev)%grid%aijpk
577  aijkm => region%levels(ilev)%grid%aijkm
578  aijkp => region%levels(ilev)%grid%aijkp
579  aimjmk => region%levels(ilev)%grid%aimjmk
580  aipjmk => region%levels(ilev)%grid%aipjmk
581  aimjpk => region%levels(ilev)%grid%aimjpk
582  aipjpk => region%levels(ilev)%grid%aipjpk
583  aimjkm => region%levels(ilev)%grid%aimjkm
584  aipjkm => region%levels(ilev)%grid%aipjkm
585  aimjkp => region%levels(ilev)%grid%aimjkp
586  aipjkp => region%levels(ilev)%grid%aipjkp
587  aijmkm => region%levels(ilev)%grid%aijmkm
588  aijpkm => region%levels(ilev)%grid%aijpkm
589  aijmkp => region%levels(ilev)%grid%aijmkp
590  aijpkp => region%levels(ilev)%grid%aijpkp
591 
592  di = 1._rfreal/(ipnend-ipnbeg)
593  dj = 1._rfreal/(jpnend-jpnbeg)
594  dk = 1._rfreal/(kpnend-kpnbeg)
595 
596  rdi2 = 1._rfreal/(di*di)
597  rdj2 = 1._rfreal/(dj*dj)
598  rdk2 = 1._rfreal/(dk*dk)
599  r2di = 1._rfreal/(2*di)
600  r2dj = 1._rfreal/(2*dj)
601  r2dk = 1._rfreal/(2*dk)
602 
603 ! start -----------------------------------------------------------------------
604 
605  DO k=kdnbeg,kdnend
606  DO j=jdnbeg,jdnend
607  DO i=idnbeg,idnend
608  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
609 
610  al11 = dot_product( stui(:,ijkn),stui(:,ijkn) )
611  al12 = dot_product( stui(:,ijkn),stuj(:,ijkn) )
612  al13 = dot_product( stui(:,ijkn),stuk(:,ijkn) )
613  al22 = dot_product( stuj(:,ijkn),stuj(:,ijkn) )
614  al23 = dot_product( stuj(:,ijkn),stuk(:,ijkn) )
615  al33 = dot_product( stuk(:,ijkn),stuk(:,ijkn) )
616  au11 = al22*al33 - al23*al23
617  au12 = al13*al23 - al12*al33
618  au13 = al12*al23 - al13*al22
619  au22 = al11*al33 - al13*al13
620  au23 = al13*al12 - al11*al23
621  au33 = al11*al22 - al12*al12
622 
623  q1 = au11*pmat(xcoord,1,ijkn) + 2._rfreal*au12*pmat(xcoord,2,ijkn) + &
624  2._rfreal*au13*pmat(xcoord,3,ijkn) + au22*pmat(xcoord,4,ijkn) + &
625  2._rfreal*au23*pmat(xcoord,5,ijkn) + au33*pmat(xcoord,6,ijkn)
626  q2 = au11*pmat(ycoord,1,ijkn) + 2._rfreal*au12*pmat(ycoord,2,ijkn) + &
627  2._rfreal*au13*pmat(ycoord,3,ijkn) + au22*pmat(ycoord,4,ijkn) + &
628  2._rfreal*au23*pmat(ycoord,5,ijkn) + au33*pmat(ycoord,6,ijkn)
629  q3 = au11*pmat(zcoord,1,ijkn) + 2._rfreal*au12*pmat(zcoord,2,ijkn) + &
630  2._rfreal*au13*pmat(zcoord,3,ijkn) + au22*pmat(zcoord,4,ijkn) + &
631  2._rfreal*au23*pmat(zcoord,5,ijkn) + au33*pmat(zcoord,6,ijkn)
632 
633  aijk( ijkn) = 2._rfreal*(au11*rdi2 + au22*rdj2 + au33*rdk2)
634  aimjk( ijkn) = au11*rdi2 - q1*r2di
635  aipjk( ijkn) = au11*rdi2 + q1*r2di
636  aijmk( ijkn) = au22*rdj2 - q2*r2dj
637  aijpk( ijkn) = au22*rdj2 + q2*r2dj
638  aijkm( ijkn) = au33*rdk2 - q3*r2dk
639  aijkp( ijkn) = au33*rdk2 + q3*r2dk
640  aimjmk(ijkn) = 2._rfreal*au12*r2di*r2dj
641  aipjmk(ijkn) = -2._rfreal*au12*r2di*r2dj
642  aimjpk(ijkn) = -2._rfreal*au12*r2di*r2dj
643  aipjpk(ijkn) = 2._rfreal*au12*r2di*r2dj
644  aimjkm(ijkn) = 2._rfreal*au13*r2di*r2dk
645  aipjkm(ijkn) = -2._rfreal*au13*r2di*r2dk
646  aimjkp(ijkn) = -2._rfreal*au13*r2di*r2dk
647  aipjkp(ijkn) = 2._rfreal*au13*r2di*r2dk
648  aijmkm(ijkn) = 2._rfreal*au23*r2dj*r2dk
649  aijpkm(ijkn) = -2._rfreal*au23*r2dj*r2dk
650  aijmkp(ijkn) = -2._rfreal*au23*r2dj*r2dk
651  aijpkp(ijkn) = 2._rfreal*au23*r2dj*r2dk
652 !$BTEST
653 ! write(*,*) 'pmat1',region%iRegionGlobal,i,j,k,pmat(1:3,1,ijkN)
654 ! write(*,*) 'pmat2',region%iRegionGlobal,i,j,k,pmat(1:3,2,ijkN)
655 ! write(*,*) 'pmat3',region%iRegionGlobal,i,j,k,pmat(1:3,3,ijkN)
656 ! write(*,*) 'pmat4',region%iRegionGlobal,i,j,k,pmat(1:3,4,ijkN)
657 ! write(*,*) 'pmat5',region%iRegionGlobal,i,j,k,pmat(1:3,5,ijkN)
658 ! write(*,*) 'pmat6',region%iRegionGlobal,i,j,k,pmat(1:3,6,ijkN)
659 !$ETEST
660  ENDDO
661  ENDDO
662  ENDDO
663 
664 ! finalize --------------------------------------------------------------------
665 
666  CALL deregisterfunction( region%global )
667 
668 END SUBROUTINE rflo_elliptmatrixcoeffs3d
669 
670 
671 !******************************************************************************
672 !
673 ! Purpose: compute matrix coefficients for elliptic PDE grid smooting.
674 !
675 ! Description: Equation 4.20 in Handbook of Grid Generation, by
676 ! Joe F. Thompson, Bharat K. Soni, and Nigel P. Weatherill,
677 ! CRC Press 1999, ISBN 0-8493-2687-7.
678 !
679 ! Input: region = grid data of current region
680 !
681 ! Output: grid%aijk, aimjk, ... = matrix coefficients
682 !
683 ! Notes: none.
684 !
685 !******************************************************************************
686 
687 SUBROUTINE rflo_elliptmatrixcoeffs2d( region,patch,iPatch )
688 
691 
692  IMPLICIT NONE
693 
694 #include "Indexing.h"
695 
696 ! ... parameters
697  TYPE(t_region) :: region
698  TYPE(t_patch) :: patch
699  INTEGER :: ipatch
700 
701 ! ... loop variables
702  INTEGER :: i, j
703 
704 ! ... local variables
705  INTEGER :: h1, h2
706 
707  REAL(RFREAL) :: di, dj, rdi2, rdj2, r2di, r2dj, pco, qco, rco, sco, tco
708  REAL(RFREAL), POINTER :: sti(:,:,:), stj(:,:,:), pfun(:,:,:,:)
709  REAL(RFREAL), POINTER :: aij(:,:), aimj(:,:), aipj(:,:), aijm(:,:), aijp(:,:)
710  REAL(RFREAL), POINTER :: aimjm(:,:), aipjm(:,:), aimjp(:,:), aipjp(:,:)
711 
712  TYPE(t_global), POINTER :: global
713 
714 !******************************************************************************
715 
716  global => region%global
717 
718  CALL registerfunction( global,'RFLO_ElliptMatrixCoeffs2D',&
719  'RFLO_ModElliptSmoothing.F90' )
720 
721 ! get dimensions and pointers -------------------------------------------------
722 
723  sti => patch%sti
724  stj => patch%stj
725  pfun => patch%pfun
726 
727  aij => patch%aij
728  aimj => patch%aimj
729  aipj => patch%aipj
730  aijm => patch%aijm
731  aijp => patch%aijp
732  aimjm => patch%aimjm
733  aipjm => patch%aipjm
734  aimjp => patch%aimjp
735  aipjp => patch%aipjp
736 
737  h1 = patch%l1end-patch%l1beg+2
738  h2 = patch%l2end-patch%l2beg+2
739 
740  di = 1._rfreal/(h1-1)
741  dj = 1._rfreal/(h2-1)
742 
743  rdi2 = 1._rfreal/(di*di)
744  rdj2 = 1._rfreal/(dj*dj)
745  r2di = 1._rfreal/(2*di)
746  r2dj = 1._rfreal/(2*dj)
747 
748 ! start -----------------------------------------------------------------------
749 
750  DO j=1,h2
751  DO i=1,h1
752  pco = dot_product( stj(:,i,j),stj(:,i,j) )
753  qco = dot_product( sti(:,i,j),stj(:,i,j) )
754  rco = dot_product( sti(:,i,j),sti(:,i,j) )
755  sco = pco*pfun(1,1,i,j) - 2._rfreal*qco*pfun(1,2,i,j) + rco*pfun(1,3,i,j)
756  tco = pco*pfun(2,1,i,j) - 2._rfreal*qco*pfun(2,2,i,j) + rco*pfun(2,3,i,j)
757 
758  aij( i,j) = 2._rfreal*(pco*rdi2 + rco*rdj2)
759  aimj( i,j) = pco*rdi2 - sco*r2di
760  aipj( i,j) = pco*rdi2 + sco*r2di
761  aijm( i,j) = rco*rdj2 - tco*r2dj
762  aijp( i,j) = rco*rdj2 + tco*r2dj
763  aimjm(i,j) = -2._rfreal*qco*r2di*r2dj
764  aipjm(i,j) = 2._rfreal*qco*r2di*r2dj
765  aimjp(i,j) = 2._rfreal*qco*r2di*r2dj
766  aipjp(i,j) = -2._rfreal*qco*r2di*r2dj
767 !$BTEST
768 ! IF ((region%iRegionGlobal == 24 .OR. &
769 ! region%iRegionGlobal == 27 .OR. &
770 ! region%iRegionGlobal == 30 .OR. &
771 ! region%iRegionGlobal == 33) .AND. iPatch==2) &
772 ! write(*,*)region%iRegionGlobal,iPatch,i,j, &
773 ! pfun(1:2,1,i,j), &
774 ! pfun(1:2,2,i,j), &
775 ! pfun(1:2,3,i,j)
776 !$ETEST
777  ENDDO ! i
778  ENDDO ! j
779 
780 ! finalize --------------------------------------------------------------------
781 
782  CALL deregisterfunction( region%global )
783 
784 END SUBROUTINE rflo_elliptmatrixcoeffs2d
785 
786 
787 !******************************************************************************
788 !
789 ! Purpose: perform one Jacobi iteration to solve for physical space grid
790 ! coordinate by elliptic PDE.
791 !
792 ! Description: the method used is Gauss-Jacobi method.
793 !
794 ! Input: region = grid data of current region
795 !
796 ! Output: xyz = new grid in physical space
797 !
798 ! Notes: none.
799 !
800 !******************************************************************************
801 
802 SUBROUTINE rflo_elliptgridjac3d( region )
803 
805 
806  IMPLICIT NONE
807 
808 #include "Indexing.h"
809 
810 ! ... parameters
811  TYPE(t_region) :: region
812 
813 ! ... loop variables
814  INTEGER :: i, j, k, l
815 
816 ! ... local variables
817  INTEGER :: ilev, inoff, ijnoff
818  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend
819  INTEGER :: ijkn, imjkn, ipjkn, ijmkn, ijpkn, ijkmn, ijkpn
820  INTEGER :: imjmkn, ipjmkn, imjpkn, ipjpkn, imjkmn, ipjkmn, imjkpn, ipjkpn
821  INTEGER :: ijmkmn, ijpkmn, ijmkpn, ijpkpn
822 
823  REAL(RFREAL) :: raijk
824  REAL(RFREAL), POINTER :: aijk(:), aimjk(:), aipjk(:), aijmk(:), aijpk(:)
825  REAL(RFREAL), POINTER :: aijkm(:), aijkp(:), aimjmk(:), aipjmk(:), aimjpk(:)
826  REAL(RFREAL), POINTER :: aipjpk(:), aimjkm(:), aipjkm(:), aimjkp(:)
827  REAL(RFREAL), POINTER :: aipjkp(:), aijmkm(:), aijpkm(:), aijmkp(:)
828  REAL(RFREAL), POINTER :: aijpkp(:), xyz(:,:), xyzold(:,:)
829 
830  TYPE(t_global), POINTER :: global
831 
832 !******************************************************************************
833 
834  global => region%global
835 
836  CALL registerfunction( global,'RFLO_ElliptGridJac3D',&
837  'RFLO_ModElliptSmoothing.F90' )
838 
839 ! get dimensions and pointers -------------------------------------------------
840 
841  ilev = 1
842  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend, &
843  jpnbeg,jpnend,kpnbeg,kpnend )
844  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
845 
846  aijk => region%levels(ilev)%grid%aijk
847  aimjk => region%levels(ilev)%grid%aimjk
848  aipjk => region%levels(ilev)%grid%aipjk
849  aijmk => region%levels(ilev)%grid%aijmk
850  aijpk => region%levels(ilev)%grid%aijpk
851  aijkm => region%levels(ilev)%grid%aijkm
852  aijkp => region%levels(ilev)%grid%aijkp
853  aimjmk => region%levels(ilev)%grid%aimjmk
854  aipjmk => region%levels(ilev)%grid%aipjmk
855  aimjpk => region%levels(ilev)%grid%aimjpk
856  aipjpk => region%levels(ilev)%grid%aipjpk
857  aimjkm => region%levels(ilev)%grid%aimjkm
858  aipjkm => region%levels(ilev)%grid%aipjkm
859  aimjkp => region%levels(ilev)%grid%aimjkp
860  aipjkp => region%levels(ilev)%grid%aipjkp
861  aijmkm => region%levels(ilev)%grid%aijmkm
862  aijpkm => region%levels(ilev)%grid%aijpkm
863  aijmkp => region%levels(ilev)%grid%aijmkp
864  aijpkp => region%levels(ilev)%grid%aijpkp
865 
866  xyz => region%levels(ilev)%grid%xyz
867  xyzold => region%levels(ilev)%grid%xyzOld
868 
869 ! perform Gauss-Jacobi --------------------------------------------------------
870 
871  DO k=kpnbeg+1,kpnend-1
872  DO j=jpnbeg+1,jpnend-1
873  DO i=ipnbeg+1,ipnend-1
874  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
875  imjkn = indijk(i-1,j ,k ,inoff,ijnoff)
876  ipjkn = indijk(i+1,j ,k ,inoff,ijnoff)
877  ijmkn = indijk(i ,j-1,k ,inoff,ijnoff)
878  ijpkn = indijk(i ,j+1,k ,inoff,ijnoff)
879  ijkmn = indijk(i ,j ,k-1,inoff,ijnoff)
880  ijkpn = indijk(i ,j ,k+1,inoff,ijnoff)
881  imjmkn = indijk(i-1,j-1,k ,inoff,ijnoff)
882  ipjmkn = indijk(i+1,j-1,k ,inoff,ijnoff)
883  imjpkn = indijk(i-1,j+1,k ,inoff,ijnoff)
884  ipjpkn = indijk(i+1,j+1,k ,inoff,ijnoff)
885  imjkmn = indijk(i-1,j ,k-1,inoff,ijnoff)
886  ipjkmn = indijk(i+1,j ,k-1,inoff,ijnoff)
887  imjkpn = indijk(i-1,j ,k+1,inoff,ijnoff)
888  ipjkpn = indijk(i+1,j ,k+1,inoff,ijnoff)
889  ijmkmn = indijk(i ,j-1,k-1,inoff,ijnoff)
890  ijpkmn = indijk(i ,j+1,k-1,inoff,ijnoff)
891  ijmkpn = indijk(i ,j-1,k+1,inoff,ijnoff)
892  ijpkpn = indijk(i ,j+1,k+1,inoff,ijnoff)
893 
894  raijk = 1._rfreal/aijk(ijkn)
895 
896  DO l = xcoord,zcoord
897  xyz(l,ijkn) = aimjk( ijkn)*xyzold(l, imjkn) + &
898  aipjk( ijkn)*xyzold(l, ipjkn) + &
899  aijmk( ijkn)*xyzold(l, ijmkn) + &
900  aijpk( ijkn)*xyzold(l, ijpkn) + &
901  aijkm( ijkn)*xyzold(l, ijkmn) + &
902  aijkp( ijkn)*xyzold(l, ijkpn) + &
903  aimjmk( ijkn)*xyzold(l,imjmkn) + &
904  aipjmk( ijkn)*xyzold(l,ipjmkn) + &
905  aimjpk( ijkn)*xyzold(l,imjpkn) + &
906  aipjpk( ijkn)*xyzold(l,ipjpkn) + &
907  aimjkm( ijkn)*xyzold(l,imjkmn) + &
908  aipjkm( ijkn)*xyzold(l,ipjkmn) + &
909  aimjkp( ijkn)*xyzold(l,imjkpn) + &
910  aipjkp( ijkn)*xyzold(l,ipjkpn) + &
911  aijmkm( ijkn)*xyzold(l,ijmkmn) + &
912  aijpkm( ijkn)*xyzold(l,ijpkmn) + &
913  aijmkp( ijkn)*xyzold(l,ijmkpn) + &
914  aijpkp( ijkn)*xyzold(l,ijpkpn)
915  xyz(l,ijkn) = xyz(l,ijkn)*raijk
916  ENDDO ! l
917  ENDDO ! i
918  ENDDO ! j
919  ENDDO ! k
920 
921 ! finalize --------------------------------------------------------------------
922 
923  CALL deregisterfunction( region%global )
924 
925 END SUBROUTINE rflo_elliptgridjac3d
926 
927 
928 !******************************************************************************
929 !
930 ! Purpose: perform one Gauss-Seidel iteration to solve for physical space grid
931 ! coordinate by elliptic PDE.
932 !
933 ! Description: the method used is Gauss-Seidel method.
934 !
935 ! Input: region = grid data of current region
936 !
937 ! Output: xyz = new grid in physical space
938 !
939 ! Notes: none.
940 !
941 !******************************************************************************
942 
943 SUBROUTINE rflo_elliptgridgauss3d( region )
944 
946 
947  IMPLICIT NONE
948 
949 #include "Indexing.h"
950 
951 ! ... parameters
952  TYPE(t_region) :: region
953 
954 ! ... loop variables
955  INTEGER :: i, j, k, l
956 
957 ! ... local variables
958  INTEGER :: ilev, inoff, ijnoff
959  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend
960  INTEGER :: ijkn, imjkn, ipjkn, ijmkn, ijpkn, ijkmn, ijkpn
961  INTEGER :: imjmkn, ipjmkn, imjpkn, ipjpkn, imjkmn, ipjkmn, imjkpn, ipjkpn
962  INTEGER :: ijmkmn, ijpkmn, ijmkpn, ijpkpn
963 
964  REAL(RFREAL) :: raijk
965  REAL(RFREAL), POINTER :: aijk(:), aimjk(:), aipjk(:), aijmk(:), aijpk(:)
966  REAL(RFREAL), POINTER :: aijkm(:), aijkp(:), aimjmk(:), aipjmk(:), aimjpk(:)
967  REAL(RFREAL), POINTER :: aipjpk(:), aimjkm(:), aipjkm(:), aimjkp(:)
968  REAL(RFREAL), POINTER :: aipjkp(:), aijmkm(:), aijpkm(:), aijmkp(:)
969  REAL(RFREAL), POINTER :: aijpkp(:), xyz(:,:)
970 
971  TYPE(t_global), POINTER :: global
972 
973 !******************************************************************************
974 
975  global => region%global
976 
977  CALL registerfunction( global,'RFLO_ElliptGridGauss3D',&
978  'RFLO_ModElliptSmoothing.F90' )
979 
980 ! get dimensions and pointers -------------------------------------------------
981 
982  ilev = 1
983  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend, &
984  jpnbeg,jpnend,kpnbeg,kpnend )
985  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
986 
987  aijk => region%levels(ilev)%grid%aijk
988  aimjk => region%levels(ilev)%grid%aimjk
989  aipjk => region%levels(ilev)%grid%aipjk
990  aijmk => region%levels(ilev)%grid%aijmk
991  aijpk => region%levels(ilev)%grid%aijpk
992  aijkm => region%levels(ilev)%grid%aijkm
993  aijkp => region%levels(ilev)%grid%aijkp
994  aimjmk => region%levels(ilev)%grid%aimjmk
995  aipjmk => region%levels(ilev)%grid%aipjmk
996  aimjpk => region%levels(ilev)%grid%aimjpk
997  aipjpk => region%levels(ilev)%grid%aipjpk
998  aimjkm => region%levels(ilev)%grid%aimjkm
999  aipjkm => region%levels(ilev)%grid%aipjkm
1000  aimjkp => region%levels(ilev)%grid%aimjkp
1001  aipjkp => region%levels(ilev)%grid%aipjkp
1002  aijmkm => region%levels(ilev)%grid%aijmkm
1003  aijpkm => region%levels(ilev)%grid%aijpkm
1004  aijmkp => region%levels(ilev)%grid%aijmkp
1005  aijpkp => region%levels(ilev)%grid%aijpkp
1006 
1007  xyz => region%levels(ilev)%grid%xyz
1008 
1009 ! perform Gauss iteration ----------------------------------------------------
1010 
1011  DO k=kpnbeg+1,kpnend-1
1012  DO j=jpnbeg+1,jpnend-1
1013  DO i=ipnbeg+1,ipnend-1
1014  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1015  imjkn = indijk(i-1,j ,k ,inoff,ijnoff)
1016  ipjkn = indijk(i+1,j ,k ,inoff,ijnoff)
1017  ijmkn = indijk(i ,j-1,k ,inoff,ijnoff)
1018  ijpkn = indijk(i ,j+1,k ,inoff,ijnoff)
1019  ijkmn = indijk(i ,j ,k-1,inoff,ijnoff)
1020  ijkpn = indijk(i ,j ,k+1,inoff,ijnoff)
1021  imjmkn = indijk(i-1,j-1,k ,inoff,ijnoff)
1022  ipjmkn = indijk(i+1,j-1,k ,inoff,ijnoff)
1023  imjpkn = indijk(i-1,j+1,k ,inoff,ijnoff)
1024  ipjpkn = indijk(i+1,j+1,k ,inoff,ijnoff)
1025  imjkmn = indijk(i-1,j ,k-1,inoff,ijnoff)
1026  ipjkmn = indijk(i+1,j ,k-1,inoff,ijnoff)
1027  imjkpn = indijk(i-1,j ,k+1,inoff,ijnoff)
1028  ipjkpn = indijk(i+1,j ,k+1,inoff,ijnoff)
1029  ijmkmn = indijk(i ,j-1,k-1,inoff,ijnoff)
1030  ijpkmn = indijk(i ,j+1,k-1,inoff,ijnoff)
1031  ijmkpn = indijk(i ,j-1,k+1,inoff,ijnoff)
1032  ijpkpn = indijk(i ,j+1,k+1,inoff,ijnoff)
1033 
1034  raijk = 1._rfreal/aijk(ijkn)
1035 
1036  DO l = xcoord,zcoord
1037  xyz(l,ijkn) = aimjk( ijkn)*xyz(l, imjkn) + &
1038  aipjk( ijkn)*xyz(l, ipjkn) + &
1039  aijmk( ijkn)*xyz(l, ijmkn) + &
1040  aijpk( ijkn)*xyz(l, ijpkn) + &
1041  aijkm( ijkn)*xyz(l, ijkmn) + &
1042  aijkp( ijkn)*xyz(l, ijkpn) + &
1043  aimjmk( ijkn)*xyz(l,imjmkn) + &
1044  aipjmk( ijkn)*xyz(l,ipjmkn) + &
1045  aimjpk( ijkn)*xyz(l,imjpkn) + &
1046  aipjpk( ijkn)*xyz(l,ipjpkn) + &
1047  aimjkm( ijkn)*xyz(l,imjkmn) + &
1048  aipjkm( ijkn)*xyz(l,ipjkmn) + &
1049  aimjkp( ijkn)*xyz(l,imjkpn) + &
1050  aipjkp( ijkn)*xyz(l,ipjkpn) + &
1051  aijmkm( ijkn)*xyz(l,ijmkmn) + &
1052  aijpkm( ijkn)*xyz(l,ijpkmn) + &
1053  aijmkp( ijkn)*xyz(l,ijmkpn) + &
1054  aijpkp( ijkn)*xyz(l,ijpkpn)
1055  xyz(l,ijkn) = xyz(l,ijkn)*raijk
1056  ENDDO ! l
1057  ENDDO ! i
1058  ENDDO ! j
1059  ENDDO ! k
1060 
1061 ! finalize --------------------------------------------------------------------
1062 
1063  CALL deregisterfunction( region%global )
1064 
1065 END SUBROUTINE rflo_elliptgridgauss3d
1066 
1067 
1068 !******************************************************************************
1069 !
1070 ! Purpose: perform one SOR iteration to solve for physical space grid
1071 ! coordinate by elliptic PDE.
1072 !
1073 ! Description: the method used is SOR method.
1074 !
1075 ! Input: region = grid data of current region
1076 !
1077 ! Output: xyz = new grid in physical space
1078 !
1079 ! Notes: none.
1080 !
1081 !******************************************************************************
1082 
1083 SUBROUTINE rflo_elliptgridsor3d( region )
1084 
1086 
1087  IMPLICIT NONE
1088 
1089 #include "Indexing.h"
1090 
1091 ! ... parameters
1092  TYPE(t_region) :: region
1093 
1094 ! ... loop variables
1095  INTEGER :: i, j, k, l
1096 
1097 ! ... local variables
1098  INTEGER :: ilev, inoff, ijnoff
1099  INTEGER :: ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend
1100  INTEGER :: ijkn, imjkn, ipjkn, ijmkn, ijpkn, ijkmn, ijkpn
1101  INTEGER :: imjmkn, ipjmkn, imjpkn, ipjpkn, imjkmn, ipjkmn, imjkpn, ipjkpn
1102  INTEGER :: ijmkmn, ijpkmn, ijmkpn, ijpkpn
1103 
1104  REAL(RFREAL) :: raijk, xyztemp, omega, omomg
1105  REAL(RFREAL), POINTER :: aijk(:), aimjk(:), aipjk(:), aijmk(:), aijpk(:)
1106  REAL(RFREAL), POINTER :: aijkm(:), aijkp(:), aimjmk(:), aipjmk(:), aimjpk(:)
1107  REAL(RFREAL), POINTER :: aipjpk(:), aimjkm(:), aipjkm(:), aimjkp(:)
1108  REAL(RFREAL), POINTER :: aipjkp(:), aijmkm(:), aijpkm(:), aijmkp(:)
1109  REAL(RFREAL), POINTER :: aijpkp(:), xyz(:,:)
1110 
1111  TYPE(t_global), POINTER :: global
1112 
1113 !******************************************************************************
1114 
1115  global => region%global
1116 
1117  CALL registerfunction( global,'RFLO_ElliptGridSOR3D',&
1118  'RFLO_ModElliptSmoothing.F90' )
1119 
1120 ! get dimensions and pointers -------------------------------------------------
1121 
1122  ilev = 1
1123  CALL rflo_getdimensphysnodes( region,ilev,ipnbeg,ipnend, &
1124  jpnbeg,jpnend,kpnbeg,kpnend )
1125  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
1126 
1127  aijk => region%levels(ilev)%grid%aijk
1128  aimjk => region%levels(ilev)%grid%aimjk
1129  aipjk => region%levels(ilev)%grid%aipjk
1130  aijmk => region%levels(ilev)%grid%aijmk
1131  aijpk => region%levels(ilev)%grid%aijpk
1132  aijkm => region%levels(ilev)%grid%aijkm
1133  aijkp => region%levels(ilev)%grid%aijkp
1134  aimjmk => region%levels(ilev)%grid%aimjmk
1135  aipjmk => region%levels(ilev)%grid%aipjmk
1136  aimjpk => region%levels(ilev)%grid%aimjpk
1137  aipjpk => region%levels(ilev)%grid%aipjpk
1138  aimjkm => region%levels(ilev)%grid%aimjkm
1139  aipjkm => region%levels(ilev)%grid%aipjkm
1140  aimjkp => region%levels(ilev)%grid%aimjkp
1141  aipjkp => region%levels(ilev)%grid%aipjkp
1142  aijmkm => region%levels(ilev)%grid%aijmkm
1143  aijpkm => region%levels(ilev)%grid%aijpkm
1144  aijmkp => region%levels(ilev)%grid%aijmkp
1145  aijpkp => region%levels(ilev)%grid%aijpkp
1146 
1147  xyz => region%levels(ilev)%grid%xyz
1148 
1149 ! perform point SOR -----------------------------------------------------------
1150 
1151  omega = 1.3_rfreal
1152  omomg = 1._rfreal - omega
1153 
1154  DO k=kpnbeg+1,kpnend-1
1155  DO j=jpnbeg+1,jpnend-1
1156  DO i=ipnbeg+1,ipnend-1
1157  ijkn = indijk(i ,j ,k ,inoff,ijnoff)
1158  imjkn = indijk(i-1,j ,k ,inoff,ijnoff)
1159  ipjkn = indijk(i+1,j ,k ,inoff,ijnoff)
1160  ijmkn = indijk(i ,j-1,k ,inoff,ijnoff)
1161  ijpkn = indijk(i ,j+1,k ,inoff,ijnoff)
1162  ijkmn = indijk(i ,j ,k-1,inoff,ijnoff)
1163  ijkpn = indijk(i ,j ,k+1,inoff,ijnoff)
1164  imjmkn = indijk(i-1,j-1,k ,inoff,ijnoff)
1165  ipjmkn = indijk(i+1,j-1,k ,inoff,ijnoff)
1166  imjpkn = indijk(i-1,j+1,k ,inoff,ijnoff)
1167  ipjpkn = indijk(i+1,j+1,k ,inoff,ijnoff)
1168  imjkmn = indijk(i-1,j ,k-1,inoff,ijnoff)
1169  ipjkmn = indijk(i+1,j ,k-1,inoff,ijnoff)
1170  imjkpn = indijk(i-1,j ,k+1,inoff,ijnoff)
1171  ipjkpn = indijk(i+1,j ,k+1,inoff,ijnoff)
1172  ijmkmn = indijk(i ,j-1,k-1,inoff,ijnoff)
1173  ijpkmn = indijk(i ,j+1,k-1,inoff,ijnoff)
1174  ijmkpn = indijk(i ,j-1,k+1,inoff,ijnoff)
1175  ijpkpn = indijk(i ,j+1,k+1,inoff,ijnoff)
1176 
1177  raijk = 1._rfreal/aijk(ijkn)
1178 
1179  DO l = xcoord,zcoord
1180  xyztemp = aimjk( ijkn)*xyz(l, imjkn) + &
1181  aipjk( ijkn)*xyz(l, ipjkn) + &
1182  aijmk( ijkn)*xyz(l, ijmkn) + &
1183  aijpk( ijkn)*xyz(l, ijpkn) + &
1184  aijkm( ijkn)*xyz(l, ijkmn) + &
1185  aijkp( ijkn)*xyz(l, ijkpn) + &
1186  aimjmk( ijkn)*xyz(l,imjmkn) + &
1187  aipjmk( ijkn)*xyz(l,ipjmkn) + &
1188  aimjpk( ijkn)*xyz(l,imjpkn) + &
1189  aipjpk( ijkn)*xyz(l,ipjpkn) + &
1190  aimjkm( ijkn)*xyz(l,imjkmn) + &
1191  aipjkm( ijkn)*xyz(l,ipjkmn) + &
1192  aimjkp( ijkn)*xyz(l,imjkpn) + &
1193  aipjkp( ijkn)*xyz(l,ipjkpn) + &
1194  aijmkm( ijkn)*xyz(l,ijmkmn) + &
1195  aijpkm( ijkn)*xyz(l,ijpkmn) + &
1196  aijmkp( ijkn)*xyz(l,ijmkpn) + &
1197  aijpkp( ijkn)*xyz(l,ijpkpn)
1198  xyztemp = xyztemp*raijk
1199  xyz(l,ijkn) = omega*xyztemp+omomg*xyz(l,ijkn)
1200  ENDDO ! l
1201  ENDDO ! i
1202  ENDDO ! j
1203  ENDDO ! k
1204 
1205 ! finalize --------------------------------------------------------------------
1206 
1207  CALL deregisterfunction( region%global )
1208 
1209 END SUBROUTINE rflo_elliptgridsor3d
1210 
1211 
1212 !******************************************************************************
1213 !
1214 ! Purpose: perform Jacobi iteration to solve for physical space grid
1215 ! coordinate on patches by elliptic PDE.
1216 !
1217 ! Description: the method used is Gauss-Jacobi method.
1218 !
1219 ! Input: region = grid data of current region
1220 !
1221 ! Output: xyz = new grid in physical space
1222 !
1223 ! Notes: none.
1224 !
1225 !******************************************************************************
1226 
1227 SUBROUTINE rflo_elliptgridjac2d( region,patch,iPatch )
1228 
1231 
1232  IMPLICIT NONE
1233 
1234 #include "Indexing.h"
1235 
1236 ! ... parameters
1237  TYPE(t_region) :: region
1238  TYPE(t_patch) :: patch
1239  INTEGER :: ipatch
1240 
1241 ! ... loop variables
1242  INTEGER :: i, j, k, l, iter
1243 
1244 ! ... local variables
1245  INTEGER :: ilev, lbound, h1, h2, ng1, ng2, ijkn
1246  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, inoff, ijnoff
1247 
1248  REAL(RFREAL) :: raij, dx, dy, dz, resid
1249  REAL(RFREAL), POINTER :: aij(:,:), aimj(:,:), aipj(:,:), aijm(:,:), aijp(:,:)
1250  REAL(RFREAL), POINTER :: aimjm(:,:), aipjm(:,:), aimjp(:,:), aipjp(:,:)
1251  REAL(RFREAL), POINTER :: xyz(:,:), st(:,:,:), stold(:,:,:)
1252 
1253  TYPE(t_global), POINTER :: global
1254 
1255 !******************************************************************************
1256 
1257  global => region%global
1258 
1259  CALL registerfunction( global,'RFLO_ElliptGridJac2D',&
1260  'RFLO_ModElliptSmoothing.F90' )
1261 
1262 ! get dimensions --------------------------------------------------------------
1263 
1264  ilev = 1
1265  lbound = patch%lbound
1266 
1267  CALL rflo_getpatchindicesnodes( region,patch,ilev, &
1268  ibeg,iend,jbeg,jend,kbeg,kend )
1269  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
1270 
1271  h1 = patch%l1end-patch%l1beg+2
1272  h2 = patch%l2end-patch%l2beg+2
1273 
1274 ! get pointers ----------------------------------------------------------------
1275 
1276  xyz => region%levels(ilev)%grid%xyz
1277  st => patch%st
1278  stold => patch%stOld
1279 
1280  aij => patch%aij
1281  aimj => patch%aimj
1282  aipj => patch%aipj
1283  aijm => patch%aijm
1284  aijp => patch%aijp
1285  aimjm => patch%aimjm
1286  aipjm => patch%aipjm
1287  aimjp => patch%aimjp
1288  aipjp => patch%aipjp
1289 
1290 ! copy xyz to stOld -----------------------------------------------------------
1291 
1292  DO k=kbeg,kend
1293  DO j=jbeg,jend
1294  DO i=ibeg,iend
1295  ijkn = indijk(i,j,k,inoff,ijnoff)
1296  IF (lbound==1 .OR. lbound==2) THEN
1297  ng1 = j - jbeg + 1
1298  ng2 = k - kbeg + 1
1299  ELSE IF (lbound==3 .OR. lbound==4) THEN
1300  ng1 = k - kbeg + 1
1301  ng2 = i - ibeg + 1
1302  ELSE IF (lbound==5 .OR. lbound==6) THEN
1303  ng1 = i - ibeg + 1
1304  ng2 = j - jbeg + 1
1305  ENDIF
1306  stold(1,ng1,ng2) = xyz(xcoord,ijkn)
1307  stold(2,ng1,ng2) = xyz(ycoord,ijkn)
1308  stold(3,ng1,ng2) = xyz(zcoord,ijkn)
1309  ENDDO
1310  ENDDO
1311  ENDDO
1312 
1313 ! perform Gauss-Jacobi --------------------------------------------------------
1314 
1315  DO iter = 1,global%moveGridSiter
1316 
1317  CALL rflo_gridphysgrad2d( region,patch )
1318  CALL rflo_elliptmatrixcoeffs2d( region,patch,ipatch )
1319 
1320  resid = 0._rfreal
1321  DO j=2,h2-1
1322  DO i=2,h1-1
1323 
1324  raij = 1._rfreal/aij(i,j)
1325 
1326  DO l = 1,3
1327  st(l,i,j) = aimj( i,j)*stold(l,i-1,j ) + &
1328  aipj( i,j)*stold(l,i+1,j ) + &
1329  aijm( i,j)*stold(l,i ,j-1) + &
1330  aijp( i,j)*stold(l,i ,j+1) + &
1331  aimjm(i,j)*stold(l,i-1,j-1) + &
1332  aipjm(i,j)*stold(l,i+1,j-1) + &
1333  aimjp(i,j)*stold(l,i-1,j+1) + &
1334  aipjp(i,j)*stold(l,i+1,j+1)
1335  st(l,i,j) = st(l,i,j)*raij
1336  ENDDO ! l
1337 
1338  dx = st(1,i,j) - stold(1,i,j)
1339  dy = st(2,i,j) - stold(2,i,j)
1340  dz = st(3,i,j) - stold(3,i,j)
1341  resid = resid + dx*dx + dy*dy + dz*dz
1342  ENDDO ! i
1343  ENDDO ! j
1344 
1345 !$BTEST
1346 ! write(*,*)region%iRegionGlobal,iPatch,iter,SQRT(resid)
1347 !$ETEST
1348 
1349  stold = st
1350 
1351 ! - copy st back to xyz -------------------------------------------------------
1352 
1353  DO k=kbeg,kend
1354  DO j=jbeg,jend
1355  DO i=ibeg,iend
1356  ijkn = indijk(i,j,k,inoff,ijnoff)
1357  IF (lbound==1 .OR. lbound==2) THEN
1358  ng1 = j - jbeg + 1
1359  ng2 = k - kbeg + 1
1360  ELSE IF (lbound==3 .OR. lbound==4) THEN
1361  ng1 = k - kbeg + 1
1362  ng2 = i - ibeg + 1
1363  ELSE IF (lbound==5 .OR. lbound==6) THEN
1364  ng1 = i - ibeg + 1
1365  ng2 = j - jbeg + 1
1366  ENDIF
1367  xyz(xcoord,ijkn) = st(1,ng1,ng2)
1368  xyz(ycoord,ijkn) = st(2,ng1,ng2)
1369  xyz(zcoord,ijkn) = st(3,ng1,ng2)
1370  ENDDO ! i
1371  ENDDO ! j
1372  ENDDO ! k
1373 
1374  ENDDO ! iter
1375 
1376 ! finalize --------------------------------------------------------------------
1377 
1378  CALL deregisterfunction( region%global )
1379 
1380 END SUBROUTINE rflo_elliptgridjac2d
1381 
1382 
1383 !******************************************************************************
1384 !
1385 ! Purpose: perform Gauss-Seidel iteration to solve for physical space grid
1386 ! coordinate on patches by elliptic PDE.
1387 !
1388 ! Description: the method used is Gauss-Seidel method.
1389 !
1390 ! Input: region = grid data of current region
1391 !
1392 ! Output: xyz = new grid in physical space
1393 !
1394 ! Notes: none.
1395 !
1396 !******************************************************************************
1397 
1398 SUBROUTINE rflo_elliptgridgauss2d( region,patch,iPatch )
1399 
1402 
1403  IMPLICIT NONE
1404 
1405 #include "Indexing.h"
1406 
1407 ! ... parameters
1408  TYPE(t_region) :: region
1409  TYPE(t_patch) :: patch
1410  INTEGER :: ipatch
1411 
1412 ! ... loop variables
1413  INTEGER :: i, j, k, l, iter
1414 
1415 ! ... local variables
1416  INTEGER :: ilev, lbound, h1, h2, ng1, ng2, ijkn
1417  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, inoff, ijnoff
1418 
1419  REAL(RFREAL) :: raij, dx, dy, dz, resid
1420  REAL(RFREAL), POINTER :: aij(:,:), aimj(:,:), aipj(:,:), aijm(:,:), aijp(:,:)
1421  REAL(RFREAL), POINTER :: aimjm(:,:), aipjm(:,:), aimjp(:,:), aipjp(:,:)
1422  REAL(RFREAL), POINTER :: xyz(:,:), st(:,:,:), stold(:,:,:)
1423 
1424  TYPE(t_global), POINTER :: global
1425 
1426 !******************************************************************************
1427 
1428  global => region%global
1429 
1430  CALL registerfunction( global,'RFLO_ElliptGridGauss2D',&
1431  'RFLO_ModElliptSmoothing.F90' )
1432 
1433 ! get dimensions --------------------------------------------------------------
1434 
1435  ilev = 1
1436  lbound = patch%lbound
1437 
1438  CALL rflo_getpatchindicesnodes( region,patch,ilev, &
1439  ibeg,iend,jbeg,jend,kbeg,kend )
1440  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
1441 
1442  h1 = patch%l1end-patch%l1beg+2
1443  h2 = patch%l2end-patch%l2beg+2
1444 
1445 ! get pointers ----------------------------------------------------------------
1446 
1447  xyz => region%levels(ilev)%grid%xyz
1448  st => patch%st
1449  stold => patch%stOld
1450 
1451  aij => patch%aij
1452  aimj => patch%aimj
1453  aipj => patch%aipj
1454  aijm => patch%aijm
1455  aijp => patch%aijp
1456  aimjm => patch%aimjm
1457  aipjm => patch%aipjm
1458  aimjp => patch%aimjp
1459  aipjp => patch%aipjp
1460 
1461 ! copy xyz to st --------------------------------------------------------------
1462 
1463  DO k=kbeg,kend
1464  DO j=jbeg,jend
1465  DO i=ibeg,iend
1466  ijkn = indijk(i,j,k,inoff,ijnoff)
1467  IF (lbound==1 .OR. lbound==2) THEN
1468  ng1 = j - jbeg + 1
1469  ng2 = k - kbeg + 1
1470  ELSE IF (lbound==3 .OR. lbound==4) THEN
1471  ng1 = k - kbeg + 1
1472  ng2 = i - ibeg + 1
1473  ELSE IF (lbound==5 .OR. lbound==6) THEN
1474  ng1 = i - ibeg + 1
1475  ng2 = j - jbeg + 1
1476  ENDIF
1477  st(1,ng1,ng2) = xyz(xcoord,ijkn)
1478  st(2,ng1,ng2) = xyz(ycoord,ijkn)
1479  st(3,ng1,ng2) = xyz(zcoord,ijkn)
1480  ENDDO
1481  ENDDO
1482  ENDDO
1483 
1484 ! perform Gauss-Seidel --------------------------------------------------------
1485 
1486  DO iter = 1,global%moveGridSiter
1487 
1488  CALL rflo_gridphysgrad2d( region,patch )
1489  CALL rflo_elliptmatrixcoeffs2d( region,patch,ipatch )
1490 
1491  stold = st
1492 
1493  resid = 0._rfreal
1494  DO j=2,h2-1
1495  DO i=2,h1-1
1496 
1497  raij = 1._rfreal/aij(i,j)
1498 
1499  DO l = 1,3
1500  st(l,i,j) = aimj( i,j)*st(l,i-1,j ) + &
1501  aipj( i,j)*st(l,i+1,j ) + &
1502  aijm( i,j)*st(l,i ,j-1) + &
1503  aijp( i,j)*st(l,i ,j+1) + &
1504  aimjm(i,j)*st(l,i-1,j-1) + &
1505  aipjm(i,j)*st(l,i+1,j-1) + &
1506  aimjp(i,j)*st(l,i-1,j+1) + &
1507  aipjp(i,j)*st(l,i+1,j+1)
1508  st(l,i,j) = st(l,i,j)*raij
1509  ENDDO ! l
1510 
1511  dx = st(1,i,j) - stold(1,i,j)
1512  dy = st(2,i,j) - stold(2,i,j)
1513  dz = st(3,i,j) - stold(3,i,j)
1514  resid = resid + dx*dx + dy*dy + dz*dz
1515  ENDDO ! i
1516  ENDDO ! j
1517 
1518 !$BTEST
1519 ! write(*,*)region%iRegionGlobal,iPatch,iter,SQRT(resid)
1520 !$ETEST
1521 
1522 ! - copy st back to xyz -------------------------------------------------------
1523 
1524  DO k=kbeg,kend
1525  DO j=jbeg,jend
1526  DO i=ibeg,iend
1527  ijkn = indijk(i,j,k,inoff,ijnoff)
1528  IF (lbound==1 .OR. lbound==2) THEN
1529  ng1 = j - jbeg + 1
1530  ng2 = k - kbeg + 1
1531  ELSE IF (lbound==3 .OR. lbound==4) THEN
1532  ng1 = k - kbeg + 1
1533  ng2 = i - ibeg + 1
1534  ELSE IF (lbound==5 .OR. lbound==6) THEN
1535  ng1 = i - ibeg + 1
1536  ng2 = j - jbeg + 1
1537  ENDIF
1538  xyz(xcoord,ijkn) = st(1,ng1,ng2)
1539  xyz(ycoord,ijkn) = st(2,ng1,ng2)
1540  xyz(zcoord,ijkn) = st(3,ng1,ng2)
1541  ENDDO ! i
1542  ENDDO ! j
1543  ENDDO ! k
1544 
1545  ENDDO ! iter
1546 
1547 ! finalize --------------------------------------------------------------------
1548 
1549  CALL deregisterfunction( region%global )
1550 
1551 END SUBROUTINE rflo_elliptgridgauss2d
1552 
1553 
1554 !******************************************************************************
1555 !
1556 ! Purpose: perform SOR iteration to solve for physical space grid
1557 ! coordinate on patches by elliptic PDE.
1558 !
1559 ! Description: the method used is SOR method.
1560 !
1561 ! Input: region = grid data of current region
1562 !
1563 ! Output: xyz = new grid in physical space
1564 !
1565 ! Notes: none.
1566 !
1567 !******************************************************************************
1568 
1569 SUBROUTINE rflo_elliptgridsor2d( region,patch,iPatch )
1570 
1573 
1574  IMPLICIT NONE
1575 
1576 #include "Indexing.h"
1577 
1578 ! ... parameters
1579  TYPE(t_region) :: region
1580  TYPE(t_patch) :: patch
1581  INTEGER :: ipatch
1582 
1583 ! ... loop variables
1584  INTEGER :: i, j, k, l, iter
1585 
1586 ! ... local variables
1587  INTEGER :: ilev, lbound, h1, h2, ng1, ng2, ijkn
1588  INTEGER :: ibeg, iend, jbeg, jend, kbeg, kend, inoff, ijnoff
1589 
1590  REAL(RFREAL) :: raij, stemp, omega, omomg, dx, dy, dz, resid
1591  REAL(RFREAL), POINTER :: aij(:,:), aimj(:,:), aipj(:,:), aijm(:,:), aijp(:,:)
1592  REAL(RFREAL), POINTER :: aimjm(:,:), aipjm(:,:), aimjp(:,:), aipjp(:,:)
1593  REAL(RFREAL), POINTER :: xyz(:,:), st(:,:,:), stold(:,:,:)
1594 
1595  TYPE(t_global), POINTER :: global
1596 
1597 !******************************************************************************
1598 
1599  global => region%global
1600 
1601  CALL registerfunction( global,'RFLO_ElliptGridSOR2D',&
1602  'RFLO_ModElliptSmoothing.F90' )
1603 
1604 ! get dimensions and parameters -----------------------------------------------
1605 
1606  omega = 1.5_rfreal
1607  omomg = 1._rfreal - omega
1608 
1609  ilev = 1
1610  lbound = patch%lbound
1611 
1612  CALL rflo_getpatchindicesnodes( region,patch,ilev, &
1613  ibeg,iend,jbeg,jend,kbeg,kend )
1614  CALL rflo_getnodeoffset( region,ilev,inoff,ijnoff )
1615 
1616  h1 = patch%l1end-patch%l1beg+2
1617  h2 = patch%l2end-patch%l2beg+2
1618 
1619 ! get pointers ----------------------------------------------------------------
1620 
1621  xyz => region%levels(ilev)%grid%xyz
1622  st => patch%st
1623  stold => patch%stOld
1624 
1625  aij => patch%aij
1626  aimj => patch%aimj
1627  aipj => patch%aipj
1628  aijm => patch%aijm
1629  aijp => patch%aijp
1630  aimjm => patch%aimjm
1631  aipjm => patch%aipjm
1632  aimjp => patch%aimjp
1633  aipjp => patch%aipjp
1634 
1635 ! copy xyz to st --------------------------------------------------------------
1636 
1637  DO k=kbeg,kend
1638  DO j=jbeg,jend
1639  DO i=ibeg,iend
1640  ijkn = indijk(i,j,k,inoff,ijnoff)
1641  IF (lbound==1 .OR. lbound==2) THEN
1642  ng1 = j - jbeg + 1
1643  ng2 = k - kbeg + 1
1644  ELSE IF (lbound==3 .OR. lbound==4) THEN
1645  ng1 = k - kbeg + 1
1646  ng2 = i - ibeg + 1
1647  ELSE IF (lbound==5 .OR. lbound==6) THEN
1648  ng1 = i - ibeg + 1
1649  ng2 = j - jbeg + 1
1650  ENDIF
1651  st(1,ng1,ng2) = xyz(xcoord,ijkn)
1652  st(2,ng1,ng2) = xyz(ycoord,ijkn)
1653  st(3,ng1,ng2) = xyz(zcoord,ijkn)
1654  ENDDO
1655  ENDDO
1656  ENDDO
1657 
1658 ! perform point SOR ----------------------------------------------------------
1659 
1660  DO iter = 1,global%moveGridSiter
1661 
1662  CALL rflo_gridphysgrad2d( region,patch )
1663  CALL rflo_elliptmatrixcoeffs2d( region,patch,ipatch )
1664 
1665  stold = st
1666 
1667  resid = 0._rfreal
1668  DO j=2,h2-1
1669  DO i=2,h1-1
1670 
1671  raij = 1._rfreal/aij(i,j)
1672 
1673  DO l = 1,3
1674  stemp = aimj( i,j)*st(l,i-1,j ) + &
1675  aipj( i,j)*st(l,i+1,j ) + &
1676  aijm( i,j)*st(l,i ,j-1) + &
1677  aijp( i,j)*st(l,i ,j+1) + &
1678  aimjm(i,j)*st(l,i-1,j-1) + &
1679  aipjm(i,j)*st(l,i+1,j-1) + &
1680  aimjp(i,j)*st(l,i-1,j+1) + &
1681  aipjp(i,j)*st(l,i+1,j+1)
1682  stemp = stemp*raij
1683  st(l,i,j) = omega*stemp+omomg*st(l,i,j)
1684  ENDDO ! l
1685 
1686  dx = st(1,i,j) - stold(1,i,j)
1687  dy = st(2,i,j) - stold(2,i,j)
1688  dz = st(3,i,j) - stold(3,i,j)
1689  resid = resid + dx*dx + dy*dy + dz*dz
1690  ENDDO ! i
1691  ENDDO ! j
1692 
1693 !$BTEST
1694 ! write(*,*)region%iRegionGlobal,iPatch,iter,SQRT(resid)
1695 !$ETEST
1696 
1697 ! - copy st back to xyz -------------------------------------------------------
1698 
1699  DO k=kbeg,kend
1700  DO j=jbeg,jend
1701  DO i=ibeg,iend
1702  ijkn = indijk(i,j,k,inoff,ijnoff)
1703  IF (lbound==1 .OR. lbound==2) THEN
1704  ng1 = j - jbeg + 1
1705  ng2 = k - kbeg + 1
1706  ELSE IF (lbound==3 .OR. lbound==4) THEN
1707  ng1 = k - kbeg + 1
1708  ng2 = i - ibeg + 1
1709  ELSE IF (lbound==5 .OR. lbound==6) THEN
1710  ng1 = i - ibeg + 1
1711  ng2 = j - jbeg + 1
1712  ENDIF
1713  xyz(xcoord,ijkn) = st(1,ng1,ng2)
1714  xyz(ycoord,ijkn) = st(2,ng1,ng2)
1715  xyz(zcoord,ijkn) = st(3,ng1,ng2)
1716  ENDDO ! i
1717  ENDDO ! j
1718  ENDDO ! k
1719 
1720  ENDDO ! iter
1721 
1722 ! finalize --------------------------------------------------------------------
1723 
1724  CALL deregisterfunction( region%global )
1725 
1726 END SUBROUTINE rflo_elliptgridsor2d
1727 
1728 
1729 ! ******************************************************************************
1730 ! End
1731 ! ******************************************************************************
1732 
1733 END MODULE rflo_modelliptsmoothing
1734 
1735 
1736 ! ******************************************************************************
1737 !
1738 ! RCS Revision history:
1739 !
1740 ! $Log: RFLO_ModElliptSmoothing.F90,v $
1741 ! Revision 1.15 2008/12/06 08:44:15 mtcampbe
1742 ! Updated license.
1743 !
1744 ! Revision 1.14 2008/11/19 22:17:27 mtcampbe
1745 ! Added Illinois Open Source License/Copyright
1746 !
1747 ! Revision 1.13 2006/03/14 04:37:13 wasistho
1748 ! removed obsolete RFLO_ElliptFlatPatch
1749 !
1750 ! Revision 1.12 2006/03/12 10:28:33 wasistho
1751 ! defined boundFlat in ElliptFlatPatch
1752 !
1753 ! Revision 1.11 2006/03/08 06:39:35 wasistho
1754 ! added moveGridSiter
1755 !
1756 ! Revision 1.10 2006/03/05 19:06:22 wasistho
1757 ! set computational space coordinates from initial grid
1758 !
1759 ! Revision 1.9 2006/03/03 04:17:00 wasistho
1760 ! elliptic smoothings on coords instead of deformations
1761 !
1762 ! Revision 1.8 2006/03/03 00:47:41 wasistho
1763 ! added elliptGridSmooRegion
1764 !
1765 ! Revision 1.7 2006/03/02 00:23:54 wasistho
1766 ! prepared elliptic pde grid motion
1767 !
1768 ! Revision 1.6 2006/02/11 03:42:14 wasistho
1769 ! added ModMoveGridElliptGlo/Fra
1770 !
1771 ! Revision 1.5 2006/02/09 00:25:11 wasistho
1772 ! removed unused pointer xyzTemp
1773 !
1774 ! Revision 1.4 2006/02/08 07:51:56 wasistho
1775 ! increased eps in flatPatch
1776 !
1777 ! Revision 1.3 2005/12/07 08:46:15 wasistho
1778 ! added stuff for surface mesh motion EPDE
1779 !
1780 ! Revision 1.2 2005/12/05 10:49:44 wasistho
1781 ! added RFLO_ElliptFlatPatch
1782 !
1783 ! Revision 1.1 2005/12/03 09:39:47 wasistho
1784 ! initial import
1785 !
1786 !
1787 !
1788 ! ******************************************************************************
1789 
1790 
1791 
1792 
1793 
1794 
1795 
1796 
1797 
1798 
1799 
1800 
1801 
1802 
1803 
1804 
1805 
subroutine, public rflo_elliptgridsmoo(regions, resid)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ibeg
subroutine rflo_copygeometrydummy(region)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE jpnbeg
subroutine, public rflo_elliptgridsmooregion(region, resid)
subroutine rflo_calccellcentroids(region)
j indices k indices k
Definition: Indexing.h:6
NT dx
subroutine, public rflo_elliptgridsor2d(region, patch, iPatch)
NT q1
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE kpnbeg
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
subroutine, public rflo_elliptgridjac2d(region, patch, iPatch)
subroutine rflo_extrapolategeometry(region)
subroutine rflo_exchangednoderecv(region, regionSrc, patch, patchSrc, average, dNode)
subroutine rflo_changeinteriorgrid(region, boundMoved, edgeMoved, arcLen12, arcLen34, arcLen56, xyzOld, xyz)
subroutine, public rflo_gridphysgrad2d(region, patch)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE jdnbeg
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE jpnend
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE idnend
subroutine rflo_getnodeoffset(region, iLev, iNodeOffset, ijNodeOffset)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE jdnend
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE idnbeg
subroutine, public rflo_gridphysgrad3d(region)
Definition: patch.h:74
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE knode iend
subroutine rflo_calcfacevectors(region)
subroutine, public rflo_elliptgridgauss3d(region)
subroutine, public rflo_elliptgridpatch(region, patch)
subroutine, public rflo_elliptmatrixcoeffs3d(region)
subroutine rflo_exchangegeometry(regions)
subroutine, public rflo_elliptgridsor3d(region)
subroutine rflo_generatecoarsegrids(region)
subroutine rflo_getpatchindicesnodes(region, patch, iLev, ibeg, iend, jbeg, jend, kbeg, kend)
blockLoc i
Definition: read.cpp:79
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ipnbeg
subroutine, public rflo_elliptgridjac3d(region)
subroutine rflo_clearsendrequests(regions, iReg, geometry)
RT dz() const
Definition: Direction_3.h:133
subroutine rflo_getdimensphysnodes(region, iLev, ipnbeg, ipnend, jpnbeg, jpnend, kpnbeg, kpnend)
j indices j
Definition: Indexing.h:6
NT dy
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE knode jend
subroutine rflo_exchangednodecopy(region, regionSrc, patch, patchSrc, average, dNode, dNodeSrc)
subroutine rflo_getdimensdummynodes(region, iLev, idnbeg, idnend, jdnbeg, jdnend, kdnbeg, kdnend)
subroutine, public rflo_elliptmatrixcoeffs2d(region, patch, iPatch)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE knode jbeg
long double dot_product(pnt vec1, pnt vec2)
subroutine grid(bp)
Definition: setup_py.f90:257
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE knode kbeg
subroutine deregisterfunction(global)
Definition: ModError.F90:469
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ipnend
subroutine rflo_exchangednodesend(region, regionSrc, patch, dNode)
subroutine, public rflo_elliptgridgauss2d(region, patch, iPatch)
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE kdnbeg
subroutine rflo_calcfacecentroids(region)