1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23!=============================================================================== 24! Function: 25! --------- 26 27!> \file richards.f90 28!> 29!> \brief This routine solves the Richards equation, then compute the new 30!> velocities deducted from the gradients of hydraulic head and 31!> from the permeability. 32!> These velocities are used for post-processing, calculation of dispersion 33!> coefficients, convergence criterion of Newton scheme... but not 34!> for transport. 35!> In order to ensure the exact conservation of mass, the mass fluxes are 36!> computed following the procedure of the standard subroutine resopv 37!> (See the <a href="../../theory.pdf#resopv"><b>resopv</b></a> 38!> section of the theory guide for more informations). 39!> 40!> The hydraulic head is contained in the pressure array, while the velocities 41!> are contained in the velocity arrays. 42!> 43!> The Richards equation for underground flows is : 44!> d theta(h) / dt - div( K(h) grad(h) ) = 0 (or source(h) in very 45!> particular cases). 46!> This equation can also be written, denoting C(h) = d theta(h)/dh : 47!> C(h) dh/dt - div( K(h) grad(h) ) = C(h) dh/dt - d theta(h) / dt (1) 48!> The right hand term is close to zero and is not considered in ESTEL. 49!> We consider it here for exact mass conservation. 50!> 51!> (1) is solved with a 'Newton scheme' handled by tridim. The structure used 52!> for this Newton scheme is the same as the one used in the loop over 53!> nterup for standard flows. 54!> For going from time step n to time step n+1, we use sub-iterations 55!> indexed by m : 56!> C(h^n) (h^(n+1,m+1)-h^n)/detla t - div( K(h^(n+1,m)) grad(h^(n+1,m+1)) ) 57!> = C(h^n) (h^(n+1,m)-h^n)/detla t - (theta(h^(n+1,m))-theta(h(n)))/delta t. 58!> These sub-iterations, if they converge, converge to the solution i 59!> of the problem. 60!> 61!> The Darcy velocity q is then computed thanks to the relation : 62!> q = -K(h) grad(h). 63!> 64!> This routine is essentially inspired from navstv, resopv and 65!> cs_equation_iterative_solve_scalar. 66!> 67!> Please refer to the <a href="../../theory.pdf#groundwater"><b>groundwater flows</b></a> 68!> section of the theory guide for more theoretical informations. 69!------------------------------------------------------------------------------- 70 71!------------------------------------------------------------------------------- 72! Arguments 73!______________________________________________________________________________. 74! mode name role ! 75!______________________________________________________________________________! 76!> \param[in] icvrge indicator of convergence 77!> \param[in] dt time step (per cell) 78!_______________________________________________________________________________ 79 80 81subroutine richards & 82 (icvrge, dt) 83 84!=============================================================================== 85! Module files 86!=============================================================================== 87 88use paramx 89use dimens, only: ndimfb 90use numvar 91use entsor 92use cstphy 93use cstnum 94use optcal 95use pointe 96use albase 97use parall 98use period 99use mesh 100use field 101use field_operator 102use cs_c_bindings 103use darcy_module 104use cs_c_bindings 105 106!=============================================================================== 107 108implicit none 109 110! Arguments 111 112integer icvrge 113double precision, pointer, dimension(:) :: dt 114 115! Local variables 116 117integer iccocg, inc , iel, isou, init 118integer imrgrp, nswrgp, imligp, iwarnp 119integer imucpp, ircflp, isweep, isym, lchain 120integer ndircp, niterf, nswmpr 121integer iflmas, iflmab, iesize, idiffp, iconvp, ibsize, imvisp 122integer fid 123integer iflid , iflwgr, f_dim, f_id0, iwgrp, iprev, iitsm 124 125double precision epsrgp, climgp, extrap 126double precision thetap, xdu, xdv, xdw, xnrmul 127double precision relaxp, residu, ressol, epsilp, rnormp 128 129character(len=80) :: chaine 130 131type(solving_info) sinfo 132type(var_cal_opt) :: vcopt_p 133 134double precision rvoid(1) 135 136double precision, dimension(:,:), allocatable :: gradp 137double precision, allocatable, dimension(:,:) :: xam, weighf, uvwk 138double precision, allocatable, dimension(:) :: dpvar, dam, presa, rhs, w1, weighb 139double precision, allocatable, dimension(:) :: rovsdt, rhs0, viscf, viscb 140 141double precision, dimension(:), pointer :: imasfl, bmasfl 142double precision, dimension(:), pointer :: coefa_p, coefb_p 143double precision, dimension(:), pointer :: coefaf_p, coefbf_p 144double precision, dimension(:), pointer :: cpro_permeability, cproa_capacity 145double precision, dimension(:), pointer :: cproa_sat, cpro_sat 146double precision, dimension(:,:), pointer :: cpro_permeability_6 147double precision, dimension(:,:), pointer :: cvar_vel 148double precision, dimension(:), pointer :: cvar_pr, cvara_pr 149double precision, dimension(:), pointer :: cpro_prtot 150double precision, dimension(:,:), pointer :: cpro_wgrec_v 151double precision, dimension(:), pointer :: cpro_wgrec_s 152 153!=============================================================================== 154! 0. Initialization 155!=============================================================================== 156 157if (darcy_convergence_criterion.eq.0) then 158 allocate(uvwk(1,ncelet)) 159else 160 allocate(uvwk(3,ncelet)) 161endif 162allocate(dpvar(ncelet)) 163allocate(presa(ncelet)) 164allocate(rhs(ncelet), rovsdt(ncelet), rhs0(ncelet)) 165allocate(viscf(nfac), viscb(ndimfb)) 166 167! Map field arrays 168call field_get_val_v(ivarfl(iu), cvar_vel) 169 170if (darcy_anisotropic_permeability.eq.0) then 171 call field_get_val_s_by_name('permeability', cpro_permeability) 172else 173 call field_get_id('permeability', fid) 174 call field_get_val_v(fid, cpro_permeability_6) 175endif 176 177call field_get_val_prev_s_by_name('capacity', cproa_capacity) 178call field_get_val_prev_s_by_name('saturation', cproa_sat) 179call field_get_val_s_by_name('saturation', cpro_sat) 180call field_get_val_s(ivarfl(ipr), cvar_pr) 181call field_get_val_prev_s(ivarfl(ipr), cvara_pr) 182 183sinfo%nbivar = 0 184 185! Preparation of convergence criterion for Newton scheme 186if (nterup.gt.1) then 187 !$omp parallel do 188 do iel = 1,ncel 189 if (darcy_convergence_criterion.eq.0) then 190 uvwk(1,iel) = cvar_pr(iel) 191 else 192 uvwk(1,iel) = cvar_vel(1,iel) 193 uvwk(2,iel) = cvar_vel(2,iel) 194 uvwk(3,iel) = cvar_vel(3,iel) 195 endif 196 enddo 197 xnrmul = 0.d0 198 !$omp parallel do reduction(+:xnrmul) 199 do iel = 1, ncel 200 if (darcy_convergence_criterion.eq.0) then 201 xnrmul = xnrmul +cvar_pr(iel)**2*volume(iel) 202 else 203 xnrmul = xnrmul +(cvar_vel(1,iel)**2 & 204 + cvar_vel(2,iel)**2 & 205 + cvar_vel(2,iel)**2)*volume(iel) 206 endif 207 enddo 208 if (irangp.ge.0) then 209 call parsom (xnrmul) 210 endif 211 xnrmu0 = sqrt(xnrmul) 212endif 213 214! Index of the field 215iflid = ivarfl(ipr) 216 217call field_get_key_struct_var_cal_opt(ivarfl(ipr), vcopt_p) 218 219if (vcopt_p%iwgrec.eq.1) then 220 ! Id weighting field for gradient 221 call field_get_key_int(iflid, kwgrec, iflwgr) 222 call field_get_dim(iflwgr, f_dim) 223 if (f_dim.gt.1) then 224 call field_get_val_v(iflwgr, cpro_wgrec_v) 225 else 226 call field_get_val_s(iflwgr, cpro_wgrec_s) 227 endif 228endif 229 230! Preparation of writing 231call field_get_name(ivarfl(ipr), chaine) 232lchain = 16 233 234f_id0 = -1 235iwgrp = 0 236if (vcopt_p%iwgrec.eq.1) iwgrp = 1 237 238! --- Boundary conditions 239call field_get_coefa_s(ivarfl(ipr), coefa_p) 240call field_get_coefb_s(ivarfl(ipr), coefb_p) 241call field_get_coefaf_s(ivarfl(ipr), coefaf_p) 242call field_get_coefbf_s(ivarfl(ipr), coefbf_p) 243 244! --- Physical fluxes 245call field_get_key_int(ivarfl(ipr), kimasf, iflmas) 246call field_get_key_int(ivarfl(ipr), kbmasf, iflmab) 247call field_get_val_s(iflmas, imasfl) 248call field_get_val_s(iflmab, bmasfl) 249 250 251! --- Solving options. Copied from resopv. 252isym = 1 253if (darcy_unsteady.eq.1) isym = 2 254 255! Matrix block size. Copied from reopv. 256ibsize = 1 257iesize = 1 258 259allocate(dam(ncelet), xam(isym,nfac)) 260 261!=============================================================================== 262! 1. Building of the linear system to solve 263!=============================================================================== 264 265! Unsteady term 266if (darcy_unsteady.eq.1) then 267 do iel = 1, ncel 268 rovsdt(iel) = volume(iel)*cproa_capacity(iel)/dt(iel) 269 enddo 270else 271 do iel = 1, ncel 272 rovsdt(iel) = 0.d0 273 enddo 274endif 275 276! We keep the capacity in explicit because the Newton scheme converges 277! much better like this. 278 279! Computation of diffusion coefficients at the centers of faces 280if (darcy_anisotropic_permeability.eq.0) then 281 imvisp = vcopt_p%imvisf 282 call viscfa & 283 ( imvisp , & 284 cpro_permeability , & 285 viscf , viscb ) 286 if (vcopt_p%iwgrec.eq.1) then 287 ! Weighting for gradient 288 do iel = 1, ncel 289 cpro_wgrec_s(iel) = cpro_permeability(iel) 290 enddo 291 call synsca(cpro_wgrec_s) 292 endif 293else if (darcy_anisotropic_permeability.eq.1) then 294 allocate(weighf(2,nfac)) 295 allocate(weighb(ndimfb)) 296 iwarnp = vcopt_p%iwarni 297 call vitens & 298 ( cpro_permeability_6 , iwarnp , & 299 weighf , weighb , & 300 viscf , viscb ) 301 if (vcopt_p%iwgrec.eq.1) then 302 ! Weighting for gradient 303 do iel = 1, ncel 304 do isou = 1, 6 305 cpro_wgrec_v(isou,iel) = cpro_permeability_6(isou,iel) 306 enddo 307 enddo 308 call syntis(cpro_wgrec_v) 309 endif 310endif 311 312iconvp = 0 ! no convection 313idiffp = 1 ! diffusion 314ndircp = vcopt_p%ndircl ! no diagonal stepped aside 315thetap = 1.d0 ! implicit scheme 316imucpp = 0 ! do not multiply the convectiv term by anything 317 318! Note that the computed matrix is always the same in the loop over iterns 319! (loop of the Newton scheme, handled by tridim) 320! Another option would be to store the arrays dam and xam in a global array 321! rather than computing the matrix at each call to 'richards'. 322call matrix & 323( iconvp , idiffp , ndircp , isym , & 324 thetap , imucpp , & 325 coefb_p , coefbf_p , rovsdt , & 326 imasfl , bmasfl , viscf , viscb , & 327 rvoid , dam , xam ) 328 329!=============================================================================== 330! 2. Solving (Loop over the non-orthogonalities) 331!=============================================================================== 332 333! For this part, see documentation on resopv. 334! The way used to handle non-orthogonalities is the same. 335! The difference is the instationnary term. 336 337nswmpr = vcopt_p%nswrsm 338 339! Initialization of dpvar for avoiding warnings 340do iel = 1, ncel 341 dpvar(iel) = 0.d0 342enddo 343 344! We compute the first residue (rhs) 345 346relaxp = vcopt_p%relaxv ! relaxation for loop over isweep 347isweep = 1 ! counter for non-orthogonalities 348iccocg = 1 ! no calculation of cocg. What does it mean? 349init = 1 ! it is an initialization of the residue 350inc = 1 ! 0 increment, 1 otherwise 351imrgrp = vcopt_p%imrgra 352nswrgp = vcopt_p%nswrgr ! number of sweeps for gradients reconstruction 353imligp = vcopt_p%imligr 354iwgrp = vcopt_p%iwgrec 355iwarnp = vcopt_p%iwarni 356epsrgp = vcopt_p%epsrgr 357climgp = vcopt_p%climgr 358ircflp = vcopt_p%ircflu 359extrap = 0 360 361if (darcy_anisotropic_permeability.eq.0) then 362 363 call itrgrp & 364( ivarfl(ipr), init , inc , imrgrp , iccocg , nswrgp , imligp , iphydr , & 365 iwgrp , iwarnp , & 366 epsrgp , climgp , extrap , & 367 rvoid , & 368 cvar_pr , & 369 coefa_p , coefb_p , & 370 coefaf_p , coefbf_p , & 371 viscf , viscb , & 372 cpro_permeability, & 373 rhs ) 374 375else if (darcy_anisotropic_permeability.eq.1) then 376 377 call itrgrv & 378( ivarfl(ipr), init , inc , imrgrp , iccocg , nswrgp , imligp , ircflp , & 379 iphydr , iwgrp , iwarnp , & 380 epsrgp , climgp , extrap , & 381 rvoid , & 382 cvar_pr , & 383 coefa_p , coefb_p , & 384 coefaf_p , coefbf_p , & 385 viscf , viscb , & 386 cpro_permeability_6 , & 387 weighf , weighb , & 388 rhs ) 389 390endif 391 392do iel = 1, ncel 393 if (darcy_unsteady.eq.1) then 394 ! We take care of exact mass conservation (as precise as the convergence 395 ! of the Newton loop is good) 396 rhs0(iel) = cproa_capacity(iel)*volume(iel)/dt(iel) & 397 *cvar_pr(iel) + volume(iel)/dt(iel) & 398 *(cproa_sat(iel) - cpro_sat(iel)) 399 else 400 rhs0(iel) = 0.d0 401 endif 402enddo 403 404if (ncetsm.gt.0) then 405 !$omp parallel do private(iel) if(ncetsm > thr_n_min) 406 do iitsm = 1, ncetsm 407 iel = icetsm(iitsm) 408 rhs0(iel) = rhs0(iel) + volume(iel)*smacel(iitsm,ipr) 409 enddo 410endif 411 412do iel = 1, ncel 413 if (darcy_unsteady.eq.1) then 414 rhs(iel) = rhs0(iel) - rhs(iel) & 415 - volume(iel)/dt(iel)*cproa_capacity(iel)*cvar_pr(iel) 416 else 417 rhs(iel) = rhs0(iel) - rhs(iel) 418 endif 419enddo 420 421residu = sqrt(cs_gdot(ncel, rhs, rhs)) 422 423sinfo%rnsmbr = residu 424 425! We compute the normalisation residue, which is used as a stop criterion 426! in the loop of non-orthogonalities. We have to "normalize" the problem, 427! taking into account the boundary conditions. 428! This part is inspired from cs_equation_iterative_solve_scalar 429! (call to promav). 430 431allocate(w1(ncelet)) 432call promav(isym, ibsize, iesize, ivarfl(ipr), dam, xam, cvar_pr, w1) 433 434!$omp parallel do 435do iel = 1, ncel 436 w1(iel) = w1(iel) + rhs(iel) 437enddo 438rnormp = sqrt(cs_gdot(ncel, w1, w1)) 439 440! Free memory 441deallocate(w1) 442 443! Writing 444if (vcopt_p%iwarni.ge.2) then 445 write(nfecra,1400)chaine(1:16), isweep, residu, relaxp 446endif 447 448! Loop for non-orthogonalities 449! TODO: dynamic relaxation 450nswmpr = vcopt_p%nswrsm 451do while ( (isweep.le.nswmpr.and.residu.gt.vcopt_p%epsrsm*rnormp) & 452 .or. (isweep.eq.1) ) 453 ! We pass at least once to ensure exactness of mass flux 454 455 iwarnp = vcopt_p%iwarni 456 epsilp = vcopt_p%epsilo 457 ressol = residu 458 459 call sles_solve_native(ivarfl(ipr), chaine, & 460 isym, ibsize, iesize, dam, xam, & 461 epsilp, rnormp, niterf, ressol, rhs, dpvar) 462 463 if (isweep.le.nswmpr.and.residu.gt.vcopt_p%epsrsm*rnormp) then 464 do iel = 1, ncel 465 presa(iel) = cvar_pr(iel) 466 cvar_pr(iel) = presa(iel) + vcopt_p%relaxv*dpvar(iel) 467 enddo 468 469 ! If it is the last sweep, update with the total increment 470 else 471 do iel = 1, ncel 472 presa(iel) = cvar_pr(iel) 473 cvar_pr(iel) = presa(iel) + dpvar(iel) 474 enddo 475 endif 476 477 iccocg = 1 478 init = 1 479 inc = 1 480 481 if (darcy_anisotropic_permeability.eq.0) then 482 483 call itrgrp & 484 ( ivarfl(ipr), init , inc , imrgrp , iccocg , nswrgp , imligp , iphydr , & 485 iwgrp , iwarnp , & 486 epsrgp , climgp , extrap , & 487 rvoid , & 488 cvar_pr , & 489 coefa_p , coefb_p , & 490 coefaf_p , coefbf_p , & 491 viscf , viscb , & 492 cpro_permeability, & 493 rhs ) 494 495 else if (darcy_anisotropic_permeability.eq.1) then 496 497 call itrgrv & 498 ( ivarfl(ipr), init, inc, imrgrp, iccocg, nswrgp, imligp, ircflp, & 499 iphydr, iwgrp, iwarnp, epsrgp , climgp, extrap, rvoid, & 500 cvar_pr, coefa_p, coefb_p, coefaf_p, coefbf_p, & 501 viscf, viscb, & 502 cpro_permeability_6, weighf, weighb, rhs) 503 504 endif 505 506 do iel = 1, ncel 507 if (darcy_unsteady.eq.1) then 508 rhs(iel) = rhs0(iel) - rhs(iel) & 509 - volume(iel)/dt(iel)*cproa_capacity(iel)*cvar_pr(iel) 510 else 511 rhs(iel) = rhs0(iel) - rhs(iel) 512 endif 513 enddo 514 515 ! --- Convergence test 516 residu = sqrt(cs_gdot(ncel, rhs, rhs)) 517 518 ! Writing 519 sinfo%nbivar = sinfo%nbivar + niterf 520 521 if (vcopt_p%iwarni.ge.2) then 522 write(nfecra,1400)chaine(1:16), isweep, residu, relaxp 523 endif 524 525 if (iwarnp.ge.3) then 526 write(nfecra,1500) chaine(1:16), isweep, residu, rnormp, niterf 527 endif 528 529 isweep = isweep + 1 530 531enddo 532! --- Reconstruction loop (end) 533 534if (abs(rnormp).gt.epzero) then 535 sinfo%resvar = residu/rnormp 536 sinfo%dervar = - sinfo%rnsmbr/rnormp !FIXME 537else 538 sinfo%resvar = 0.d0 539 sinfo%dervar = 0.d0 540endif 541 542! Writing 543if (vcopt_p%iwarni.ge.2) then 544 if (isweep.gt.nswmpr) then 545 write(nfecra,1600) chaine(1:16), nswmpr 546 endif 547endif 548 549call sles_free_native(ivarfl(ipr), chaine) 550deallocate(dam, xam, rhs) 551 552!=============================================================================== 553! 3. Updating of mass fluxes 554!=============================================================================== 555 556iccocg = 1 557init = 1 558inc = 1 559imrgrp = vcopt_p%imrgra 560nswrgp = vcopt_p%nswrgr 561imligp = vcopt_p%imligr 562iwarnp = vcopt_p%iwarni 563epsrgp = vcopt_p%epsrgr 564climgp = vcopt_p%climgr 565extrap = 0 566 567! We compute the new mass flux, taking care not to reconstruct 568! the last increment of the lopp on isweep, to ensure an 569! exact conservation of the mass. 570 571if (darcy_anisotropic_permeability.eq.0) then 572 573 call itrmas & 574 ( f_id0 , init , inc , imrgrp , iccocg , nswrgp , imligp , iphydr , & 575 iwgrp , iwarnp , & 576 epsrgp , climgp , extrap , & 577 rvoid , & 578 presa , & 579 coefa_p , coefb_p , coefaf_p , coefbf_p , & 580 viscf , viscb , & 581 cpro_permeability, & 582 imasfl , bmasfl ) 583 584 ! The last increment is not reconstructed to fullfill exactly the continuity 585 ! equation (see theory guide). 586 iccocg = 0 587 nswrgp = 0 588 inc = 0 589 init = 0 590 591 call itrmas & 592 ( f_id0 , init , inc , imrgrp , iccocg , nswrgp , imligp , iphydr , & 593 iwgrp , iwarnp , & 594 epsrgp , climgp , extrap , & 595 rvoid , & 596 dpvar , & 597 coefa_p , coefb_p , coefaf_p , coefbf_p , & 598 viscf , viscb , & 599 cpro_permeability, & 600 imasfl , bmasfl ) 601 602else if (darcy_anisotropic_permeability.eq.1) then 603 604 call itrmav & 605 ( f_id0, init , inc , imrgrp , iccocg , nswrgp , imligp , ircflp , & 606 iphydr , iwgrp , iwarnp , & 607 epsrgp , climgp , extrap , & 608 rvoid , & 609 presa , & 610 coefa_p , coefb_p , coefaf_p , coefbf_p , & 611 viscf , viscb , & 612 cpro_permeability_6 , & 613 weighf , weighb , & 614 imasfl , bmasfl ) 615 616 ! The last increment is not reconstructed to fullfill exactly the continuity 617 ! equation (see theory guide). 618 iccocg = 0 619 nswrgp = 0 620 inc = 0 621 init = 0 622 ircflp = 0 623 624 call itrmav & 625 ( f_id0, init , inc , imrgrp , iccocg , nswrgp , imligp , ircflp , & 626 iphydr , iwgrp , iwarnp , & 627 epsrgp , climgp , extrap , & 628 rvoid , & 629 dpvar , & 630 coefa_p , coefb_p , coefaf_p , coefbf_p , & 631 viscf , viscb , & 632 cpro_permeability_6 , & 633 weighf , weighb , & 634 imasfl , bmasfl ) 635 636endif 637 638!=============================================================================== 639! 4. Updating of the velocity field and the pressure head 640!=============================================================================== 641 642! We compute the gradient of hydraulique head and multiply it by the 643! cpro_permeability in order to get the new velocities. These velocities will 644! only be used for post-processing, computation of the dispersion coefficients 645! of scalars and possibly convergence criterion of the Newton scheme. The 646! transport equation uses the mass fluxes at the center of faces, and not 647! the velocities at the center of cells. 648 649if (irangp.ge.0.or.iperio.eq.1) then 650 call synsca(cvar_pr) 651endif 652 653iccocg = 1 654inc = 1 655iprev = 0 656 657allocate(gradp(3,ncelet)) 658 659! We use gradient_scalar instead of potential because iphydr is 660! not taken into account in case of Darcy calculation. 661call field_gradient_scalar(ivarfl(ipr), iprev, 0, inc, & 662 iccocg, gradp) 663 664! Computation of velocity 665if (darcy_anisotropic_permeability.eq.1) then 666 667 !$omp parallel do 668 do iel = 1, ncel 669 cvar_vel(1,iel) = - ( cpro_permeability_6(1,iel)*gradp(1,iel) & 670 + cpro_permeability_6(4,iel)*gradp(2,iel) & 671 + cpro_permeability_6(6,iel)*gradp(3,iel)) 672 cvar_vel(2,iel) = - ( cpro_permeability_6(4,iel)*gradp(1,iel) & 673 + cpro_permeability_6(2,iel)*gradp(2,iel) & 674 + cpro_permeability_6(5,iel)*gradp(3,iel)) 675 cvar_vel(3,iel) = - (cpro_permeability_6(6,iel)*gradp(1,iel) & 676 + cpro_permeability_6(5,iel)*gradp(2,iel) & 677 + cpro_permeability_6(3,iel)*gradp(3,iel)) 678 enddo 679 680else 681 682 !$omp parallel do 683 do iel = 1, ncel 684 cvar_vel(1,iel) = - cpro_permeability(iel)*gradp(1,iel) 685 cvar_vel(2,iel) = - cpro_permeability(iel)*gradp(2,iel) 686 cvar_vel(3,iel) = - cpro_permeability(iel)*gradp(3,iel) 687 enddo 688 689endif 690 691! update pressure head (h = H - z) for post-processing 692! Only used when gravity is taken into account 693if (darcy_gravity.ge.1) then 694 call field_get_val_s(iprtot, cpro_prtot) 695 !$omp parallel do 696 do iel = 1, ncel 697 cpro_prtot(iel) = cvar_pr(iel) - xyzcen(1,iel)*darcy_gravity_x & 698 - xyzcen(2,iel)*darcy_gravity_y & 699 - xyzcen(3,iel)*darcy_gravity_z 700 enddo 701 702endif 703 704!=============================================================================== 705! 5. Checking of convergence criterion for the Newton scheme 706!=============================================================================== 707 708! Computation of the convergence criterion for the Newton sheme. If 709! the convergence is reached, we set icvrge=1, which will be used by tridim 710! for stopping the loop. 711 712if (nterup.gt.1) then 713 icvrge = 1 714 xnrmul = 0.d0 715 !$omp parallel do reduction(+:xnrmul) private(xdu, xdv, xdw) 716 do iel = 1,ncel 717 if (darcy_convergence_criterion.eq.0) then 718 xdu = cvar_pr(iel) - uvwk(1,iel) 719 xnrmul = xnrmul + xdu**2*volume(iel) 720 else 721 xdu = cvar_vel(1,iel) - uvwk(1,iel) 722 xdv = cvar_vel(2,iel) - uvwk(2,iel) 723 xdw = cvar_vel(3,iel) - uvwk(3,iel) 724 xnrmul = xnrmul + (xdu**2 + xdv**2 + xdw**2)*volume(iel) 725 endif 726 enddo 727 if (irangp.ge.0) call parsom (xnrmul) 728 xnrmu = sqrt(xnrmul) 729 ! Indicator of convergence for the Newton scheme 730 if (xnrmu.ge.epsup*xnrmu0) icvrge = 0 731endif 732 733!=============================================================================== 734! 6. Finalization 735!=============================================================================== 736 737! Save convergence info 738call field_set_key_struct_solving_info(ivarfl(ipr), sinfo) 739 740! Free memory 741deallocate(gradp) 742deallocate(uvwk) 743deallocate(dpvar) 744deallocate(presa) 745deallocate(rovsdt) 746if (allocated(weighf)) deallocate(weighf, weighb) 747deallocate(viscf, viscb) 748 749!-------- 750! Formats 751!-------- 752 753 1400 format(1X,A16,' : SWEEP = ',I5,' RIGHT HAND SIDE NORM = ',E14.6, & 754 ', RELAXP = ',E14.6) 755 1500 format ( & 756 1X,A16,' : Current reconstruction sweep = ',I5 ,/,& 757' sweep residual = ',E12.5,', norm = ',E12.5 ,/,& 758' number of sweeps for solver = ',I5) 759 1600 format( & 760'@' ,/,& 761'@ @@ WARNING: ', A16,' RICHARDS' ,/,& 762'@ ========' ,/,& 763'@ Maximum number of iterations ',I10 ,' reached' ,/,& 764'@' ) 765 766!-------- 767! End 768!-------- 769 770return 771end subroutine richards 772