1!!***************************************************** 2! 3! Turbo AMG (ramg) interface functions 4! 1. V_cycle 5! 2. plain 6! 3. driver 7! 8!***************************************************** 9 10 !******************************************* 11 ! ramg_V_cycle: 12 ! One V_cycle call. 13 ! Now use SAMG, later will use own code 14 !******************************************* 15 recursive subroutine ramg_V_cycle( 16 & ramg_sol,ramg_rhs,clevel, 17 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper 18 & )!ramg_V_cycle 19 use ramg_data 20 include "common.h" 21 include "mpif.h" 22 23 !Variable Declaration 24 ! arguments 25 integer, intent(in) :: clevel 26 real(kind=8),intent(inout),dimension(amg_nshg(clevel)) 27 & :: ramg_sol,ramg_rhs 28 29 integer,intent(in),dimension(nshg+1) :: colm 30 integer,intent(in),dimension(nnz_tot) :: rowp 31 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 32 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 33 integer,intent(in),dimension(nlwork) :: ilwork 34 integer,intent(in),dimension(nshg) :: iBC,iper 35 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 36 37 38 ! local 39 real(kind=8),dimension(amg_nshg(clevel)) :: myvF,myvE 40 real(kind=8),dimension(:),allocatable :: myvC,myvCS 41 real(kind=8) :: cpusec(10) 42 43 ! In scale 44 45 if (clevel.ne.1) then 46 ramg_rhs = ramg_rhs*amg_ppeDiag(clevel)%p 47 endif 48 49 if (ramg_levelx-clevel .eq. 0) then 50 ! finest level solver 51 call ramg_direct_solve(ramg_rhs,ramg_sol, 52 & colm,rowp,lhsK,lhsP, 53 & ilwork,BC,iBC,iper) 54 !********************** 55 ! higher level smoother 56 !********************** 57 else if (ramg_levelx-clevel .gt. 0) then 58 allocate(myvC(amg_nshg(clevel+1))) 59 allocate(myvCS(amg_nshg(clevel+1))) 60 61 ! smoothing (pre) 62 myvF = 0 ! initial guess is zero 63 call ramg_smoother(clevel,ramg_sol,myvF,ramg_rhs, 64 & colm,rowp,lhsK,lhsP, 65 & ilwork,BC,iBC,iper,2) !dx1 66 67 ! restriction 68 call ramg_calcAv_g(clevel,myvF,ramg_sol,colm,rowp,lhsK,lhsP, 69 & ilwork,BC,iBC,iper,1) 70 myvF = myvF - ramg_rhs! r2 71 call ramg_calcIvFC(clevel,clevel+1,myvF,myvC) 72 ! up one level, solve 73 myvCS = 0 74 call ramg_V_cycle(myvCS,myvC,clevel+1, 75 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 76 ! prolongation 77 call ramg_calcIvCF(clevel,clevel+1,myvCS,myvF, 78 & ilwork,BC,iBC,iper) ! v2hat 79 ramg_sol = ramg_sol - myvF 80 81 ! smoothing (post) 82 call ramg_smoother(clevel,myvF,ramg_sol,ramg_rhs, 83 & colm,rowp,lhsK,lhsP, 84 & ilwork,BC,iBC,iper,1) 85 ramg_sol = myvF 86 87 deallocate(myVC) 88 deallocate(myVCS) 89 end if 90 91 ! Out scale 92 if (clevel.ne.1) then 93 ramg_sol = ramg_sol * amg_ppeDiag(clevel)%p 94 endif 95 96 return 97 98 end subroutine ramg_V_cycle 99 100 101!***************************************************** 102! ramg Cg solver 103! using ramg_V_cycle as preconditioner 104!***************************************************** 105 subroutine ramg_CG( 106 & sol,rhs,level, 107 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper,iterNum) 108 use ramg_data 109 include "common.h" 110 include "mpif.h" 111 include "auxmpi.h" 112 113 integer,intent(in) :: level 114 real(kind=8),intent(in),dimension(amg_nshg(level)) :: rhs 115 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: sol 116 117 integer,intent(in),dimension(nshg+1) :: colm 118 integer,intent(in),dimension(nnz_tot) :: rowp 119 120 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 121 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 122 123 integer,intent(in),dimension(nlwork) :: ilwork 124 integer,intent(in),dimension(nshg) :: iBC,iper 125 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 126 integer,intent(inout) :: iterNum 127 128 ! local 129 real(kind=8),dimension(amg_nshg(level)) :: cgP,cgQ,cgZ,cgR 130 real(kind=8) :: rz,rz_0,pq,alpha,beta 131 real(kind=8) :: tmp,norm_0,norm_1,norm_e,norm_c 132 real(kind=8) :: mres_0,mres_n 133 134 cgR = rhs 135 call ramg_L2_norm_p(level,cgR,1,norm_0) 136 norm_c = norm_0 137 norm_e = norm_0*1e-7 138 139 cgZ = cgR 140 141 call ramg_dot_p(level,cgZ,cgR,1,rz) 142 rz_0 = rz 143 cgP = cgZ 144 sol = 0 145 146 iterNum = 1 147 148 do while ( (norm_c .gt. norm_e).and.(iterNum.lt.500) ) 149 call ramg_calcAv_g(level,cgQ,cgP,colm,rowp,lhsK,lhsP, 150 & ilwork,BC,iBC,iper,1) 151 call ramg_dot_p(level,cgP,cgQ,1,pq) 152 alpha = rz/pq 153 sol = sol + alpha*cgP 154 cgR = cgR - alpha*cgQ 155 call ramg_L2_norm_p(level,cgR,1,tmp) 156 norm_c = tmp 157 cgZ = cgR 158 tmp = rz 159 call ramg_dot_p(level,cgZ,cgR,1,rz) 160 beta = rz/tmp 161 cgP = beta*cgP+cgZ 162 iterNum = iterNum + 1 163 enddo 164 165 end subroutine ! ramg_CG 166 167 subroutine ramg_PCG( 168 & sol,rhs,run_mode, 169 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper 170 & )!ramg_CG 171 use ramg_data 172 include "common.h" 173 include "mpif.h" 174 175 real(kind=8),intent(inout),dimension(nshg) :: sol 176 real(kind=8),intent(in),dimension(nshg) :: rhs 177 integer,intent(in) :: run_mode 178 179 integer,intent(in),dimension(nshg+1) :: colm 180 integer,intent(in),dimension(nnz_tot) :: rowp 181 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 182 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 183 integer,intent(in),dimension(nlwork) :: ilwork 184 integer,intent(in),dimension(nshg) :: iBC,iper 185 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 186 187 188 ! local 189 real(kind=8),dimension(amg_nshg(1)) :: cgP,cgQ,cgZ,cgR 190 real(kind=8) :: rz,rz_0,pq,alpha,beta 191 real(kind=8) :: tmp,norm_0,norm_p,norm_e,norm_c 192 real(kind=8) :: mres_0,mres_n 193 real(kind=8),dimension(nshg,4) :: diag 194 integer :: iterNum 195 196 ! controls 197 198 real(kind=8),dimension(maxiters) :: rconv 199 real(kind=8) :: avgconv 200 logical :: vcycle,gcycle,restart 201 character cflagv,cflagg 202 203 cgR = rhs 204 call ramg_L2_norm_p(1,cgR,1,norm_0) 205 norm_c = norm_0 206 norm_e = norm_0 * epstol(2) 207 write(*,*)'norm_0 ', norm_0 208 209 vcycle = .true. 210 gcycle = .true. 211 restart = .false. 212 213 ! amg preconditioning control 214 215 if (vcycle) then 216 call ramg_V_cycle(cgZ,cgR,1, 217 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 218 if (gcycle) then 219 call ramg_G_cycle(cgZ,cgR,1, 220 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 221 endif 222 else 223 cgZ = cgR 224 endif 225 226 call ramg_dot_p(1,cgZ,cgR,1,rz) 227 rz_0 = rz 228 229 cgP = cgZ 230 sol = 0 231 232 iterNum = 1 233 rconv(1) = 1.0 234 norm_p=norm_c 235 avgconv = 0 236 237 do while ( (norm_c.gt.norm_e).and.(iterNum.lt.maxiters)) 238 !do while ((iterNum.lt.maxiters)) 239 call ramg_calcAv_g(1,cgQ,cgP,colm,rowp,lhsK,lhsP, 240 & ilwork,BC,iBC,iper,1) 241 call ramg_dot_p(1,cgP,cgQ,1,pq) 242 alpha = rz/pq 243 sol = sol + alpha*cgP 244 cgR = cgR - alpha*cgQ 245 call ramg_L2_norm_p(1,cgR,1,tmp) 246 norm_c = tmp 247 if (iterNum.eq.1) rconv(1) = norm_p/norm_c 248 249 ! amg preconditioning control 250 if (vcycle) then 251 call ramg_V_cycle(cgZ,cgR,1, 252 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 253 if (gcycle) then 254 call ramg_G_cycle(cgZ,cgR,1, 255 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 256 endif 257 else 258 cgZ = cgR 259 endif 260 261 tmp = rz 262 call ramg_dot_p(1,cgZ,cgR,1,rz) 263 beta = rz/tmp 264 cgP = beta*cgP + cgZ 265 iterNum = iterNum + 1 266 rconv(iterNum) = norm_p/norm_c 267 norm_p = norm_c 268 cflagv='n' 269 cflagg='n' 270 if (vcycle) cflagv='v' 271 if (gcycle) cflagg='g' 272 if (myrank.eq.master) then 273 write(*,"((A7)(I4)(E14.4)(T30)(F8.4)(T40)(A1)(A1))") 274 & 'AMGCG: ', 275 & iterNum,norm_c/norm_0,rconv(iterNum),cflagv,cflagg 276 end if 277 enddo 278 279 end subroutine !ramg_PCG 280 281!******************************************* 282! ramg Interface 283! Interface for libLES 284!******************************************* 285 subroutine ramg_interface( 286 & colm,rowp,lhsK,lhsP,flowDiag, 287 & mcgR,mcgZ, 288 & ilwork,BC,iBC,iper 289 & ) 290 use ramg_data 291 include "common.h" 292 include "mpif.h" 293 integer,intent(in),dimension(nshg+1) :: colm 294 integer,intent(in),dimension(nnz_tot) :: rowp 295 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 296 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 297 real(kind=8),intent(in),dimension(nshg) :: mcgR 298 real(kind=8),intent(inout),dimension(nshg) :: mcgZ 299 integer, intent(in), dimension(nlwork) :: ilwork 300 integer, intent(in),dimension(nshg) :: iBC,iper 301 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 302 real(kind=8),dimension(nshg) :: myR 303 real(kind=8),intent(in),dimension(nshg,4) :: flowDiag 304 305 integer::i,j,k,amgmode 306 logical :: precflag 307 !character*10 :: fname 308 309 mcgZ = 0 310 amgmode = 11 311 312 ramg_winct = mod(ramg_winct+1,2) 313 314 if (irun_amg_prec.eq.0) then 315 precflag = .false. 316 else if ((irun_amg_prec.ge.2).and.(ramg_flag.eq.1)) then 317 ! first solve attempt with CG if using smart solver 318 precflag=.false. 319 else 320 precflag = .true. 321 endif 322 323 if (precflag) then !.and.(ramg_redo.ne.0)) then 324 myR = mcgR*amg_ppeDiag(1)%p 325 call ramg_driver(colm,rowp,lhsK,lhsP, 326 & nshg,nnz_tot,nflow, 327 & myR, 328 & nlwork,ilwork,ndofBC,BC,iBC,iper,mcgZ,1, 329 & ramg_eps,amgmode) 330 mcgZ = mcgZ*amg_ppeDiag(1)%p 331 else 332 mcgZ = mcgR 333 endif 334 335 end subroutine ramg_interface 336 337!******************************************* 338! ramg Driver 339! 1. extract PPE / system 340! 2. call coarsening 341! 3. solve 342!******************************************** 343 344 !************************************* 345 ! ramg_driver 346 ! Input: global matrix,# of systems 347 ! control params 348 ! Output: solution 349 !************************************* 350 subroutine ramg_driver( 351 & colm,rowp,lhsK,lhsP, 352 & nshg,nnz_tot,nflow, 353 & pperhs, 354 & nlwork,ilwork,ndofBC,BC,iBC,iper, 355 & ramg_sol,n_sol, 356 & amg_eps,amg_mode 357 & ) 358 359 use ramg_data 360 implicit none 361 362 !***********parameters************** 363 !the matrix 364 integer,intent(in) :: nshg 365 integer,intent(in) :: nnz_tot 366 integer,intent(in) :: nflow 367 !the matrix 368 integer,intent(in),dimension(nshg+1) :: colm 369 integer,intent(in),dimension(nnz_tot) :: rowp 370 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 371 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 372 !the forcing term, rhs 373 real(kind=8),intent(in),dimension(nshg) :: pperhs 374 ! the boundary info 375 integer, intent(in) :: nlwork 376 integer, intent(in), dimension(nlwork) :: ilwork 377 integer, intent(in) :: ndofBC 378 integer, intent(in),dimension(nshg) :: iBC,iper 379 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 380 !the output solution 381 integer, intent(in) :: n_sol 382 real(kind=8),intent(inout),dimension(nshg,n_sol) :: ramg_sol 383 384 !the tolerance 385 real(kind=8),intent(in) :: amg_eps 386 !control run mode 387 integer,intent(inout) :: amg_mode 388 ! my AMG parameter 389 !*********parameters end************** 390 391 !*****local variables***************** 392 real(kind=8) :: cpusec(10) 393 394 call cpu_time(cpusec(1)) 395 396 if (amg_mode .eq. 11 ) then ! solve PPE les Precondition ramg 397 call cpu_time(cpusec(2)) 398 call ramg_V_cycle(ramg_sol,pperhs,1, 399 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 400 call ramg_G_cycle(ramg_sol,pperhs,1, 401 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 402 call cpu_time(cpusec(3)) 403 ramg_time(4)=ramg_time(4)+cpusec(3)-cpusec(2) 404 else if (amg_mode .eq. 3) then ! solve PPE CG 405 call cpu_time(cpusec(2)) 406 call ramg_PCG(ramg_sol,pperhs,1, 407 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 408 call cpu_time(cpusec(3)) 409 ramg_time(4)=ramg_time(4)+cpusec(3)-cpusec(2) 410 end if 411 412 return 413 414 end subroutine ramg_driver 415!******************************************* 416! <EOF> ramg Driver 417!******************************************* 418