|
XBeach
|
00001 !============================================================================== 00002 ! MODULE SOLVER 00003 !============================================================================== 00004 00005 ! DATE AUTHOR CHANGES 00006 ! 00007 ! october 2009 Pieter Bart Smit New module 00008 00009 module solver_module 00010 00011 implicit none 00012 save 00013 00014 ! If mpi is defined, the non-hydrostatic module is NOT included in the compilation 00015 ! to avoid unwanted side effects. 00016 00017 00018 !****************************************************************************** 00019 ! INTERFACE 00020 !****************************************************************************** 00021 00022 private 00023 00024 !----------------------------- PARAMETERS ----------------------------------- 00025 include 'nh_pars.inc' 00026 00027 !----------------------------- VARIABLES ----------------------------------- 00028 00029 !--- PRIVATE VARIABLES --- 00030 logical :: initialized = .false. 00031 00032 integer(kind=iKind) :: itmea = 0 ! mean number of iterations 00033 integer(kind=iKind) :: itmin = 0 ! minimum number of iterations 00034 integer(kind=iKind) :: itmax = 0 ! maximum number of iterations 00035 integer(kind=iKind) :: ittot = 0 ! total number of iterations 00036 integer(kind=iKind) :: itcal = 0 ! total number calls 00037 integer(kind=iKind) :: itnconv = 0 ! total number of matrix calls which didn't converge 00038 00039 real(kind=rKind) :: reps = 0.005_rKind 00040 real(kind=rKind) :: alpha = 0.94_rKind 00041 integer(kind=iKind) :: maxit = 30 00042 00043 00044 real(kind=rKind),dimension(:,:) ,allocatable :: residual ! Residual vector 00045 real(kind=rKind),dimension(:,:,:),allocatable :: work ! work matrix 00046 00047 !--- PUBLIC VARIABLES --- 00048 00049 ! NONE 00050 00051 !--- PUBLIC SUBROUTINES --- 00052 00053 public solver_init !Allocates resources 00054 public solver_free !Free's resources 00055 public solver_solvemat !Solve system 00056 public solver_tridiag 00057 public solver_sip 00058 00059 !--- PRIVATE SUBROUTINES 00060 00061 ! NONE 00062 00063 contains 00064 ! 00065 !****************************************************************************** 00066 ! SUBROUTINES/FUNCTIONS 00067 !****************************************************************************** 00068 00069 ! 00070 !============================================================================== 00071 subroutine solver_init(nx,ny,par) 00072 !============================================================================== 00073 ! 00074 00075 00076 ! DATE AUTHOR CHANGES 00077 ! 00078 ! october 2009 Pieter Bart Smit New module 00079 00080 !------------------------------------------------------------------------------- 00081 ! DECLARATIONS 00082 !------------------------------------------------------------------------------- 00083 ! 00084 !-------------------------- PURPOSE ---------------------------- 00085 ! 00086 ! Initializes Solver 00087 ! 00088 00089 !-------------------------- DEPENDENCIES ---------------------------- 00090 00091 use xmpi_module, only: Halt_Program 00092 use params, only: parameters 00093 use logging_module 00094 use paramsconst 00095 00096 !-------------------------- ARGUMENTS ---------------------------- 00097 ! 00098 type(parameters),intent(in) :: par 00099 integer, intent(in) :: nx !Number of x-meshes 00100 integer, intent(in) :: ny !Number of y-meshes 00101 ! 00102 !-------------------------- LOCAL VARIABLES ---------------------------- 00103 00104 ! - NONE - 00105 00106 !------------------------------------------------------------------------------- 00107 ! IMPLEMENTATION 00108 !------------------------------------------------------------------------------- 00109 reps = par%solver_acc 00110 alpha = par%solver_urelax 00111 maxit = par%solver_maxit 00112 00113 ! If sol. met. ok -> allocate resources 00114 if (par%solver == SOLVER_SIPP) then !Solver is SIP 00115 allocate( work(5,1:nx+1,1:ny+1)); work = 0.0_rKind 00116 allocate(residual( 1:nx+1,1:ny+1)); residual = 0.0_rKind 00117 elseif (par%solver == SOLVER_TRIDIAGG) then !Solver is TRI-DIAG, check if possible 00118 allocate( work(5,1:nx+1,1:ny+1)); work = 0.0_rKind 00119 endif 00120 00121 initialized = .true. 00122 end subroutine solver_init 00123 00124 00125 00126 00127 00128 ! 00129 !============================================================================== 00130 subroutine solver_free 00131 !============================================================================== 00132 ! 00133 if (allocated(residual)) deallocate(residual) 00134 if (allocated(work)) deallocate(work) 00135 00136 end subroutine solver_free 00137 00138 00139 00140 00141 00142 00143 ! 00144 !============================================================================== 00145 subroutine solver_solvemat( amat , rhs , x , nx, ny, par) 00146 !============================================================================== 00147 ! 00148 00149 ! DATE AUTHOR CHANGES 00150 ! 00151 ! october 2009 Pieter Bart Smit New module 00152 00153 !------------------------------------------------------------------------------- 00154 ! DECLARATIONS 00155 !------------------------------------------------------------------------------- 00156 ! 00157 !-------------------------- PURPOSE ---------------------------- 00158 ! 00159 ! solves matrix 00160 ! 00161 00162 !-------------------------- DEPENDENCIES ---------------------------- 00163 00164 use xmpi_module 00165 use params, only: parameters 00166 use paramsconst 00167 00168 !-------------------------- ARGUMENTS ---------------------------- 00169 ! 00170 integer, intent(in) :: nx !Number of x-meshes 00171 integer, intent(in) :: ny !Number of y-meshes 00172 00173 real(kind=rKind),dimension(5,nx+1,ny+1) :: amat !the coefficient matrix used in the linear system 00174 real(kind=rKind),dimension(nx+1,ny+1) :: rhs !the right-hand side vector of the system of equations 00175 real(kind=rKind),dimension(nx+1,ny+1) ,intent(inout) :: x !solution of the linear system 00176 type(parameters),intent(in) :: par 00177 ! 00178 !-------------------------- LOCAL VARIABLES ---------------------------- 00179 00180 integer(kind=iKind) :: it 00181 00182 !------------------------------------------------------------------------------- 00183 ! IMPLEMENTATION 00184 !------------------------------------------------------------------------------- 00185 00186 if (.not. initialized) call solver_init(nx,ny,par) 00187 00188 if (par%solver == SOLVER_SIPP) then 00189 ! 00190 itcal = itcal+1 !Number of times the solver procedure is called 00191 residual = 0. 00192 00193 call solver_sip ( amat , rhs , x , residual , work , it ,nx, ny) !,reps) 00194 00195 ittot = ittot+it !Total number of iterations 00196 itmin = min(it,itmin) !Minimum number of iterations 00197 itmax = max(it,itmax) !Maximum number of iterations 00198 itmea = ittot/itcal !Mean number of iterations 00199 if (it>=par%solver_maxit) then 00200 itnconv = itnconv+1 !Number of times the solver did not converge 00201 endif 00202 ! 00203 elseif (par%solver == SOLVER_TRIDIAGG) then 00204 ! 00205 #ifdef USEMPI 00206 call xmpi_shift_zs(rhs) 00207 do it=1,3 00208 call xmpi_shift_zs(amat(it,:,:)) 00209 enddo 00210 #endif 00211 call solver_tridiag(amat,rhs,x,work,nx,ny) 00212 #ifdef USEMPI 00213 call xmpi_shift_zs(x) 00214 #endif 00215 ! 00216 endif 00217 00218 end subroutine solver_solvemat 00219 00220 ! 00221 !============================================================================== 00222 subroutine solver_tridiag ( amat , rhs , x ,cmat ,nx ,ny,fixshallow ) 00223 !============================================================================== 00224 ! 00225 00226 ! DATE AUTHOR CHANGES 00227 ! 00228 ! october 2009 Pieter Bart Smit New module 00229 00230 !------------------------------------------------------------------------------- 00231 ! DECLARATIONS 00232 !------------------------------------------------------------------------------- 00233 ! 00234 !-------------------------- PURPOSE ---------------------------- 00235 ! 00236 ! Solves matrix use the thomas (or tri-diagonal) algorithm. The solver is only 00237 ! applicable in the 1-DH case but is substantially faster than the SIP method in 00238 ! this case 00239 ! 00240 ! Algorithm is not the most efficient possible as the matrix Amat is still stored 00241 ! with 5 diagonals and 3 y-points. (Amat(5,s%nx+1,s%ny+1) as opposed to (Amat(3,s%nx+1)) 00242 ! This is easier to incorporate in the more general code. 00243 ! 00244 ! NOTE: rhs and matrix are not changed. 00245 ! 00246 !-------------------------- DEPENDENCIES ---------------------------- 00247 ! 00248 ! - NONE - 00249 ! 00250 !-------------------------- ARGUMENTS ---------------------------- 00251 ! 00252 integer, intent(in) :: nx !Number of x-meshes 00253 integer, intent(in) :: ny !Number of y-meshes 00254 logical, intent(in),optional :: fixshallow 00255 00256 real(kind=rKind),dimension(5,nx+1,ny+1),intent(in) :: amat !the coefficient matrix used in the linear system 00257 real(kind=rKind),dimension(nx+1,ny+1) ,intent(in) :: rhs !the right-hand side vector of the system of equations 00258 real(kind=rKind),dimension(nx+1,ny+1) ,intent(inout) :: x !solution of the linear system 00259 real(kind=rKind),dimension(1:nx+1) ,intent(inout) :: cmat !work vector 00260 00261 ! 00262 !-------------------------- LOCAL VARIABLES ---------------------------- 00263 00264 integer(kind=iKind) :: i !Index variable 00265 integer(kind=iKind) :: jindex !Index variable for superfast1D 00266 real (kind=rKind) :: fac !Auxillary variable 00267 logical :: lfixshallow 00268 00269 !------------------------------------------------------------------------------- 00270 ! IMPLEMENTATION 00271 !------------------------------------------------------------------------------- 00272 00273 if (ny>0) then 00274 jindex = 2 00275 else 00276 jindex = 1 00277 endif 00278 00279 if (present(fixshallow)) then 00280 lfixshallow = fixshallow 00281 else 00282 lfixshallow = .false. 00283 endif 00284 00285 fac = amat(1,1,jindex) 00286 if (abs(fac)<tiny(0.d0) .and. lfixshallow) then 00287 fac = sign(tiny(0.d0),fac) 00288 endif 00289 x(1,jindex) = rhs(1,jindex)/fac 00290 00291 !forward elimination 00292 do i=2,nx+1 00293 cmat(i) = amat(3,i-1,jindex)/fac 00294 fac = amat(1,i,jindex)-amat(2,i,jindex)*cmat(i) 00295 if (abs(fac)<tiny(0.d0) .and. lfixshallow) then 00296 fac = sign(tiny(0.d0),fac) 00297 endif 00298 x(i,jindex) = (rhs(i,jindex)-amat(2,i,jindex)*x(i-1,jindex))/fac 00299 enddo 00300 00301 !Backward substitution 00302 do i=nx,1,-1 00303 x(i,jindex) = x(i,jindex)-cmat(i+1)*x(i+1,jindex) 00304 enddo 00305 end subroutine solver_tridiag 00306 00307 ! 00308 !============================================================================== 00309 subroutine solver_sip ( amat , rhs , x , res , cmat , it ,nx, ny) !, acc) 00310 !============================================================================== 00311 ! 00312 ! programmer Marcel Zijlema 00313 ! 00314 ! Version 1.0 Date 01-07-2002 HYDRO01: first release 00315 ! Version 1.1 Date 17-07-2009 Adapted for XBeach (Pieter Smit) 00316 ! 00317 ! ********************************************************************** 00318 ! 00319 ! DESCRIPTION 00320 ! 00321 ! Solves system of equations for the Poisson equation 00322 ! for one layer by means of Stone's SIP solver 00323 ! ********************************************************************** 00324 ! 00325 ! INPUT / OUTPUT ARGUMENTS 00326 ! 00327 00328 use xmpi_module 00329 implicit none 00330 00331 integer(kind=iKind) ,intent(out) :: it !iteration count 00332 ! real(kind=rKind) ,intent(out) :: acc !iteration count 00333 integer, intent(in) :: nx !Number of x-meshes 00334 integer, intent(in) :: ny !Number of y-meshes 00335 00336 real(kind=rKind),dimension(5,nx+1,ny+1),intent(in) :: amat !the coefficient matrix used in the linear system 00337 real(kind=rKind),dimension(nx+1,ny+1) ,intent(in) :: rhs !the right-hand side vector of the system of equations 00338 real(kind=rKind),dimension(nx+1,ny+1) ,intent(inout) :: x !solution of the linear system 00339 real(kind=rKind),dimension(5,1:nx+1,1:ny+1),intent(inout) :: cmat !the matrix containing an ILU factorization 00340 real(kind=rKind),dimension(1:nx+1,1:ny+1) ,intent(inout) :: res !the residual vector 00341 ! 00342 ! LOCAL VARIABLES 00343 ! 00344 logical :: iconv = .false. ! indicator for convergence 00345 integer(kind=iKind) :: i ! X-direction 00346 integer(kind=iKind) :: j ! Y-direction 00347 real(kind=rKind) :: bnorm ! 2-norm of right-hand side vector 00348 real(kind=rKind) :: epslin ! required accuracy in the linear solver 00349 real(kind=rKind) :: p1 ! auxiliary factor 00350 real(kind=rKind) :: p2 ! auxiliary factor 00351 real(kind=rKind) :: p3 ! auxiliary factor 00352 real(kind=rKind) :: rnorm ! 2-norm of residual vector 00353 real(kind=rKind) :: ueps ! minimal accuracy based on machine precision 00354 integer :: imin,imax,jmin,jmax 00355 00356 ! 00357 ! PARAMETERS 00358 ! 00359 00360 real(kind=rKind),parameter :: small=1.e-15_rKind ! a small number 00361 #ifdef USEMPI 00362 logical, parameter :: dompi = .true. ! use mpi or not use mpi 00363 #endif 00364 00365 ! ********************************************************************** 00366 ! 00367 ! I/O 00368 ! 00369 ! none 00370 ! ********************************************************************** 00371 ! 00372 ! SUBROUTINES CALLED 00373 ! 00374 ! 00375 ! ********************************************************************** 00376 ! 00377 ! ERROR MESSAGES 00378 ! 00379 ! none 00380 ! ********************************************************************** 00381 ! 00382 ! PSEUDO CODE 00383 ! 00384 ! The system of equations is solved using an incomplete 00385 ! factorization technique called Strongly Implicit Procedure 00386 ! (SIP) as described in 00387 ! 00388 ! H.L. Stone 00389 ! Iterative solution of implicit approximations of 00390 ! multidimensional partial differential equations 00391 ! SIAM J. of Numer. Anal., vol. 5, 530-558, 1968 00392 ! 00393 ! This method constructs an incomplete lower-upper factorization 00394 ! that has the same sparsity as the original matrix. Hereby, a 00395 ! parameter 0 <= alpha <= 1 is used, which should be around 0.92 00396 ! (when alpha > 0.95, the method may diverge). Furthermore, 00397 ! alpha = 0 means standard ILU decomposition. 00398 ! 00399 ! Afterward, the resulting system is solved in an iterative manner 00400 ! by forward and backward substitutions. 00401 ! ********************************************************************** 00402 ! 00403 00404 imin = 2 00405 imax = nx 00406 jmin = 2 00407 jmax = ny 00408 #ifdef USEMPI 00409 if (dompi) then 00410 imin = 3 00411 jmin = 3 00412 imax = nx-1 00413 jmax = ny-1 00414 00415 if (xmpi_istop) imin = 2 00416 if (xmpi_isbot) imax = nx 00417 if (xmpi_isleft) jmin = 2 00418 if (xmpi_isright) jmax = ny 00419 endif 00420 #endif 00421 00422 it = 0 00423 iconv = .false. 00424 00425 ! --- construct L and U matrices (stored in cmat) 00426 00427 bnorm = 0. 00428 do j = jmin,jmax 00429 do i = imin,imax 00430 p1 = alpha*cmat(5,i-1,j) 00431 p2 = alpha*cmat(3,i,j-1) 00432 cmat(2,i,j) = amat(2,i,j)/(1.+p1) 00433 cmat(4,i,j) = amat(4,i,j)/(1.+p2) 00434 p1 = p1*cmat(2,i,j) 00435 p2 = p2*cmat(4,i,j) 00436 p3 = amat(1,i,j) + p1 + p2 & 00437 - cmat(2,i,j)*cmat(3,i-1,j) & 00438 - cmat(4,i,j)*cmat(5,i,j-1) & 00439 + small 00440 cmat(1,i,j) = 1./p3 00441 cmat(3,i,j) = (amat(3,i,j)-p2)*cmat(1,i,j) 00442 cmat(5,i,j) = (amat(5,i,j)-p1)*cmat(1,i,j) 00443 bnorm = bnorm + rhs(i,j)*rhs(i,j) 00444 enddo 00445 enddo 00446 00447 00448 #ifdef USEMPI 00449 if(dompi) then 00450 call xmpi_allreduce(bnorm,mpi_sum) 00451 endif 00452 #endif 00453 00454 bnorm = sqrt(bnorm) 00455 00456 epslin = reps*bnorm 00457 ueps = 1000.*tiny(0.)*bnorm 00458 if ( epslin < ueps .and. bnorm > 0. ) then 00459 epslin = ueps 00460 end if 00461 00462 ! --- solve the system by forward and backward substitutions 00463 ! in an iterative manner 00464 00465 iconv = .false. 00466 00467 do while ( .not. iconv .and. it < maxit ) 00468 00469 it = it + 1 00470 rnorm = 0. 00471 00472 do j = jmin, jmax 00473 do i = imin, imax 00474 res(i,j) = rhs(i,j)-amat(1,i,j)*x(i,j) & 00475 -amat(2,i,j)*x(i-1,j) & 00476 -amat(3,i,j)*x(i+1,j) & 00477 -amat(4,i,j)*x(i,j-1) & 00478 -amat(5,i,j)*x(i,j+1) 00479 !Calculate norm 00480 rnorm = rnorm + res(i,j)*res(i,j) 00481 res(i,j) = (res(i,j) - cmat(2,i,j)*res(i-1,j) & 00482 - cmat(4,i,j)*res(i,j-1))* & 00483 cmat(1,i,j) 00484 enddo 00485 enddo 00486 00487 #ifdef USEMPI 00488 if (dompi) then 00489 call xmpi_allreduce(rnorm,mpi_sum) 00490 endif 00491 #endif 00492 00493 rnorm=sqrt(rnorm) 00494 00495 do j = ny, 2, -1 00496 do i = nx, 2, -1 00497 res(i,j) = res(i,j) - cmat(3,i,j)*res(i+1,j) - cmat(5,i,j)*res(i,j+1) 00498 x(i,j) = x(i,j) + res(i,j) 00499 enddo 00500 enddo 00501 00502 #ifdef USEMPI 00503 if (dompi) then 00504 call xmpi_shift_zs(x) 00505 endif 00506 #endif 00507 00508 iconv = rnorm .lt. epslin 00509 00510 enddo 00511 00512 end subroutine solver_sip 00513 00514 end module solver_module
1.7.6.1
