From 3045d44bc458782bc0228e907419b22fb37ff3e6 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Sat, 14 Dec 2024 11:41:23 -0500 Subject: [PATCH] fix(prt): rework vertical tracking approach (#2066) Address #2014. This fixes a crash, avoids a potential infinite loop condition, and reworks vertical tracking behavior more generally. The aim is to produce results matching MF6.5 and MP7 by default (barring a few edge cases), while giving more control over how particles behave in dry conditions. Also some miscellaneous reorganization and cleanup. See the included dev notes document and/or mf6io for details on this change. --- .../test_prt_dry/test_mf6model[0-prtdry].npy | Bin 0 -> 18304 bytes .../test_mf6model[1-prtdry_drape].npy | Bin 0 -> 29504 bytes .../test_mf6model[2-prtdry_drop].npy | Bin 0 -> 26704 bytes .../test_mf6model[3-prtdry_stop].npy | Bin 0 -> 1504 bytes .../test_mf6model[4-prtdry_stay].npy | Bin 0 -> 14160 bytes .../test_mf6model[0-prtvor1l2r].npy | Bin 24384 -> 24384 bytes .../test_mf6model[1-prtvor1welp].npy | Bin 24384 -> 24384 bytes .../test_mf6model[2-prtvor1weli].npy | Bin 24384 -> 24384 bytes autotest/test_prt_dry.py | 474 ++++++++++++++++++ ..._prt_disv1.py => test_prt_newton_disv1.py} | 8 +- autotest/test_prt_voronoi1.py | 16 +- doc/ReleaseNotes/develop.tex | 5 +- doc/mf6io/mf6ivar/dfn/prt-prp.dfn | 9 + doc/mf6io/prt/prp.tex | 2 + doc/mf6io/prt/prt.tex | 36 ++ src/Model/ModelUtilities/TrackFile.f90 | 2 +- src/Model/ParticleTracking/prt-prp.f90 | 120 +++-- src/Model/ParticleTracking/prt.f90 | 8 +- src/Solution/ParticleTracker/CellDefn.f90 | 3 + src/Solution/ParticleTracker/Method.f90 | 130 +++-- .../ParticleTracker/MethodCellPassToBot.f90 | 13 +- .../ParticleTracker/MethodCellPollock.f90 | 26 +- .../ParticleTracker/MethodCellPollockQuad.f90 | 28 +- .../ParticleTracker/MethodCellTernary.f90 | 24 +- src/Solution/ParticleTracker/MethodDis.f90 | 128 +++-- src/Solution/ParticleTracker/MethodDisv.f90 | 206 ++++---- .../ParticleTracker/MethodSubcellPollock.f90 | 5 +- .../ParticleTracker/MethodSubcellTernary.f90 | 4 +- src/Solution/ParticleTracker/Particle.f90 | 85 ++-- src/Solution/ParticleTracker/vertical.md | 133 +++++ 30 files changed, 1134 insertions(+), 331 deletions(-) create mode 100644 autotest/__snapshots__/test_prt_dry/test_mf6model[0-prtdry].npy create mode 100644 autotest/__snapshots__/test_prt_dry/test_mf6model[1-prtdry_drape].npy create mode 100644 autotest/__snapshots__/test_prt_dry/test_mf6model[2-prtdry_drop].npy create mode 100644 autotest/__snapshots__/test_prt_dry/test_mf6model[3-prtdry_stop].npy create mode 100644 autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy create mode 100644 autotest/test_prt_dry.py rename autotest/{test_prt_disv1.py => test_prt_newton_disv1.py} (98%) create mode 100644 src/Solution/ParticleTracker/vertical.md diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[0-prtdry].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[0-prtdry].npy new file mode 100644 index 0000000000000000000000000000000000000000..5e8e184fe3129d7a79a9cb89877a64a40bc98455 GIT binary patch literal 18304 zcmciK%Wl&^7>417Gv$0fqwHLyEMQZi!XhhnpbEqSAp{n2Mb(f*c0f>s#67U#^0)vl z0;}neo#$tw%-^+L5d$;ws_0HPtZP~GS_}Yo$ox?U!HU4j%H*w3S@^!ojUo_l|6-)JoI&FQ#7guWRd3d7+-l3-wf9+Q#85p2`dLR9>j3@^U<|7wW0JP*3HBdMYo} zQ+cu9r`YE;jb5Mc6|Y~Yr}9EQl^5!%yiiZ&Z!b(3hae?DlgPid7+-l3-wf9P6zfvJ(U;gsk~56<%N1GFJ}UKp`OYM z^;BM{r}9EQm6x-Dy--i(g?cJ4)Kht(p32L)z+R}Q@WUZ|(?av`u6>Z!a?PvwPrDlgPidAS(a3-wf9sHgHmJ(U;gsk~eY?1g$NFVs_c zp`OYM^;BLi2lhfel^5!%yiiZ&g?cJ4R|0#Xp2`dLR9>j3@Z!a?Pvzx$U@z2Dd7+-l3-wf9sHgIBBd{0hsk~56<%N1G zFVs_cxf$3C^;BM{r}9EQl^5!%yxa=xg?cJ4)Kht(p2`dLR90(+sJ$_w>WUZ|(?LOqq2yMet>PvwPrDlgPid7)n0Uflhnzaz2vf|d9GXQ=D7 c?WL~gp7&bc+~32x*Ry@Iv28DPy|%sl1%eQ&O8@`> literal 0 HcmV?d00001 diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[1-prtdry_drape].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[1-prtdry_drape].npy new file mode 100644 index 0000000000000000000000000000000000000000..131af10426f9d2bcba71051bbfe23eb92412fba7 GIT binary patch literal 29504 zcmdU2d$3hi89(GB1T#}CB|}sSLsxE>eZ0P>e^K8>BXXY{@tKLM*KA%jIC1Xy4^Eo6eqhbu#Ph~aY*;sY z;yLGyPa8IG8K`Jj_p#L*DjNm{%_jy2wp4y%!^+Dm8@6uRxTdmY^OlucKC!v7ZE(%X z&6_q>v~3w&vmtH8m)GhTfyS#(%NKYBU|ww zHm^zBH-BX1z?$@W(zKaV&r7xNx$$-X8#N+_zgVds&++rT{*Q9oyhPf^7a?E#M#AsK zO7?3%&-M!G;H@3{d-%V{ds_tnFJ8%ydV2dn<5lyL#-aasUonZrV}Ins9_=;TDhrn2Bo^;zRgL9guk~YSC4}eG z=M@JEyT;8=k{t_DM;%uj_D_DE;|xjaji2N9*?#SwLz9+Ydo21N)rfrk(?as`p~SuC z8wo$p?^|i_vtRo;wikDeOR43O?+gq9==_&kk_g{4-fe%jAKs^Azx<5! zmn1Ky>b+uj7%!Axh|6&dO?#uB{#Rb)jmj_SLtH!*k7s!h$F`H-P&!+1%RKnYj^0GZ z6TvSO@5CiP+7H(=65;utVzm8;XL%6Ewv!+7%d|x7cyck?exdjyF8R@Za=on?aWS?S zZ9n2!9>lTjZ5ZAcQi-rEQHKso;SJv(l;Q(Atuop_Q@4$MnLzi0FV?`L5=<`>&9$7S1z zC*#^ZFF8XJXAf{A@Uq4q@dW0W^AYFA6r=fNjR)dMKHJ+0|2(l+#gBTUo~*Qbv-fG& zG;5|%?*&d8@6w*;YpIuhTDITe#*_K>C;vU3<$aC>!s*FphKSzWoIE%*-B0UFg!hlb zc=(lG@soM! zTW4Y4;K4$){m9$)OJ0=pPudTD*$$k)6{77I!n5O`e`3F20qV5e5%l_lFdqA(q~54c z+fF5y@k`CZH^Sci(Yy6QP>gzL*Uv4i%^UE3! z#FKgO&V{I#eyj7X!|&~C7}~{3t2f@a?Xq9pJq`J~X+Gla=iP{V$6tJS?2j_!FT@jn zdA1St^@?SfmmGib;f4H#JUI?oPrqxbu2*zE%5?xYLZ2Z#%bWZpuf44i`fOf^`ohUC zK7LvKlXz0!9y$&F`K3a%K12TL((p@A)p^-B-?^K0J<}l72`%QnFp_cH#5{*G;-!Sq3{g`m7btj)VLppKWf! zynf3v?DIz8WsL{oN&PvX@ji6{?C8XUk6$PreEO6+e%GncYr4NH>!S#MAv~Wxv0rxz z*6*td(fZ6957ejZvp&(O=SW(AI{G|L9gO<1(()!*Z>i@Wo`d~?L+yGF*oC;iw;b<2e(aBudh_WM_L~L$J)!3hyO#k!!hXb|w0O3i z{BTZ{M;-W##(M-_DBj7F?@ts4(Zg_KKq5@ z-SVUTz;_Jh+3rHL{X%#?emEC633{EY`6WW1q5NX`(SE{jTNl>1UoJ%3FNEjgC(lVe z*Mad2NceCxp!+uX5 z`g}(T@$dMj4=-!{NxnXFZh~<&X+Eothg$J`{3Op@(S&uwsq-Nt_+^a;>J#yBg6dP} zzX-gLf6{MgSLze>brF56B#S+C4$M!w!wf*qas^WlZ!!G|aNuB#IG<)52T7drXVhnLkq zEBK|)Ny@p^>b_`xS^b50GT)BwLmmI?zAAnv)X87qmz~I$7k6Rb_11Ps$6tJSv=62A z7y2pNW!-oE!bGmq=a*z&xUPi#y6>0ghg$LMILJ@ro@Vg8xdgwtCyW>JPs>m0&k+j| z=Qkw7JSCqx{^{cv!t?Qy{nqCb)Gv=F(f$&O2g{H4!}{Tkyx1{){!HrM2>XTbeEg*U z1%J=ud7r*$`-S4q@}vF8OX|@3G!7%|7sB)L!}D%^$e&-Wj`w%UfXa0V$4zPbWjm#Q zJturV(V56~x?ccU9q&g9p5-m)w3jc!ymwJI=E2+MBVRlD%*T)YQHp-jf2c>>PCSW= z#Vz_gyRMJ!F6S2?Ueo3G*yPR|U zZIP}UT7Y-At{)tIhVbY&?3Xxp9ONhY`qmb#qj#3D53{4}pFTX|*nWve$#Gy_Sd_?h zI=^E-rrJO2g@@<(dSSnpn}M6YCm`_{;V&V6mN)f@=a*Xa`3XHYj=&4WpVcS*1@Czr z%!@v6=>^XSe<2=uie1uwh~wiY`kQ_*FL8fyZ?wOJ;+?o02mM9n_gMw-O6v)+YlJ?t z;?Z9OKiwA;{f_L5_LmTzPoJ`mzrTRzeD3T;TsZvRr3#^bthD{w`_w1uY1Lz{2l`C+ z3#4vx^4XDsXL-x|=>8(+*{}3q{qXTt?1M(|V}G&nKs=d`zSaVN*;NAWUFH1g(Cljd#5UwnA%k9G{%k9aaqv}!&0fbNIfS@stn9&v~p ziwC{7>U=a8(%E0~;brv~@0N&ft75zkMw>hsPEKhl$|F{^+NopY$Kvf_!qJhBw|tEWwI=OEwI z0nhTLKBfM=p$+>kSI&;c%Np;*!~WS{P!EpPdMUzwp?LS%Px8t)i>l9_X!~Wwqy5BR zcI$hh_vkt>!hRtvPci^jya| zU-97)hjF39su~yDWqtChu4mJ{DeIPBlu-A6uN?2S;&D90BR`4zbDM$ls}gV{@W{*d z8`3A@Z$9$a!I_x<9Dnig3&kIC$&dQPeEZ_?eyi|};1|NPys1x-jr+A8)cwc^ysYs+ zJjpYyi}gJSeP7hk=P9ZX>c>jU!``PpMUS1!6S+>;e`s&Q`oW2ZBL&a$mUZ67OW>DJ z>bd6523$wzll?I+pdbJIPCQvhmlDM-VLrUQ9Pd7T5{I}Jk9K6ctiv9r{o1krdHaR% zn+Sd(JUb5Z!~EBdy1F|aFKhe}56`8%ke9miNj-mv&}S(AeEP)tc5xzlPR|7-zeMOW zgy++z)bn>O#=cTxHQu|`;2f`&wqJXn{v!Lyg%#kH&hKcSQ}sOGk%DJ=%lvxT3h;SC z-_zbR8+Ep`zvSb`{ume38}-TiY$u-N>-2pD$*0fgeEZ#Uy!-HoWBVl@CGo`n_g#wl zbH7%M9-+@r{Q2|= zJucDs>x|ZC2+yZasdvXM!94qvzR%_4m-nc_IbJJmzxKX<9q&IJ_u&NBOH|H4K6BnD zIa2T}Z}H2P4`baoYl)t(OozWX@8kLSu|LKk^i-bTsaLj3oqO*J)ngN$pMSQ5{gB(r z@m?z)$3r~wll9;7d8jA$m4NHSgAb3qY`^41Njy9sxl;LcK%ckh`;`&?8N##UAU~1& zJM?_|0@N`Pc%l4a`O$tdpC8zZ^+Q_sN?whyUkJ~~Pv-Z>p3P%Ftr%^;P=2xeXg@hW zTE=+S_Z%bi8N&1NlltuIn$J#ISY^Lc)xkJkD{a5_KK)bb&dycv&!k`dbPD3a`Tc?; z6;JdyZZ&wNzpIgxzX@1BIKNL+FMfv`&)QYupnJ8Rb1c>8HO9k!&hO><__06c8R*BK zuO)A8P=CqI!MyM%eSgKtpS9v~JRCY4Ew;wSTz^Srl@AN6Sa_32aMVapoT;|!bwyrlcf&Us_K z@PxVA0y9`$JZB_1WmA$jGw z)yQKzx>c`Ju%2-8r;lG&JX!a>K>PJz9d_Vc)Pv6Zu(je@-qfeq@sB-Nw+~Fld1nM( zD1Tag(m(Nh@^6r@)8B^{eMk6bR=i68d@L{ayif5W{DpW{Z{%g$sZZHQyXlvC`CW(; z`(eijJmwMGFYzd?f68@Q@8Z4D`O*HF6_5TYdi^}>>$YhB4B=Vc)Ti*eryu(-?`y5{ z&#|fy>W7m3QBrTbPsw)i%Ny4wxYl(;nqS0k9Q)M^PwMO65^n)r36b_~(0dfE+Jn*7%eBdB!@#U%GFGd42-?^DyAa`eDX8)PdLYe0l=v!^40l zaX06q;Fqp5Wj(!jJmw{5f2m&jlyl-!*8?xDlO)a$jm2}d&hOpV3s2_TpRCtzYi*2 zYCpN&oCPMUk-q6)jjOQX6;@P>*M zRWcF^F9r#@g19Xr+)PUg1!*NZndbdI30^?7*%TL825bnDooPGsZ6$f=RRF?eL`>(SZ-@1Af5CdZ)cBD@%kjvu{W$*i#^)w3{%JipeW|>w6u|KU(QX zCdU)v7vtG^6u6d`5I^#)-^igYaMam^&2^~1zl`-ysG{T@oYa(8ox{M_i^GYGrpehx0c6gEsw8fzgny(Q@Xy?4z%)< zcxWGhosR2^)gQT>2gd(zasABpfamI+^W*!hjq&*YB5yI`<<@Jg+Ye`QJV9U8dZ0cT ze{>UgUHg*t#dslo;BId5tUkWa`jc%Qm-o4U*iXxT*lRM_4X~r-df<3`zqLG$ajNqq zU%mR{c$UL9k6WJYmvQS$wca!M^UXc5f5`Kb@8kHi?59oT;dSefCo2wS()RL`yK zwH^xkW*&gfCA+cSBb*Zvu6+d`P(>;N&s{VuML^RVMGPiQ$0><3!!t{%`|t_P0G_tUaJ#Iya>C+qdB?zfYDh-X5+ z#Ojapvii{8TDksWypTTeYyX4LwOIE*Pv6x#LC)7&J1&pgehhOf@^NA_c#Vf)f0O4+ zh{y40>9-h9*0U+i&)!4uYl6O5zL3j#AU|(G{@1nsMt5d%yaPO|8{Au4Jki63MBEayqv)5mc+Uu*5SJZ}53|E<#e+@krp8}0LXeuj7) zj~4NXt=G4MGs}EAgtnJ2?egqA3i`U%BYsndV7>h)#EaEKh==oC?ngUH^&aA}-;NuL zpZIVD+q%R=+SzrW7&vt6FmjeNNi`O(z< zXlxhi;a71U$D^g+=ttX69{j5Lx%6;pJqzh;FF#eE+aJK+UaI#Huf6=d6MAO0fpI$8Eo?lYRGq*Y_ZJQ+=>r{6ajAM{DCpF8gu*d3T2P2Kcf6aH;;< z<=J^qpX5uIUf&KJ2H)EsLcCbMSbbQ>KCbc8^{iITmv(vc>st?hc5MUStM?EumY*Sg zh~IYAnVz39UPvF#uj$9VC4JXmLZ$m0oWC76)DO;m$j24E@aLO5(f(z=o`rbi(CV~W z^RoTq$+|M7`FZXz>d&k95HFUWA)d@@RPTQ-(S6X%mk_UN{0jQ!c4Wdmr~9qPyH*{= zxL9j-*yH;48tVOiaNPN)^t*pHe$6PP!V98G8=`Zs94C&)| zv@%ckPspFx?^f%&vN4m_PwRf05I=HgEzkB-pTw)N5jrpE{*b_n)dO{M9DbY-u(XT3-Ydm!%KBVf7x+E{e|-ucfBIRW;;+9UcHBSfiHch*vc~$rC%l{M7T%1YXtnkq1BMbxQwU_|xOPTOGx?SZn>pl6^Z*N)A>&5R%;n_TbE{#`H^Kwk{G(jK7qs4gn`GoB!PwcdQ z9@nqyi?^PI{29y7OUF+?SLxI9ms&YL+vSDyNuJKK9`tiEul_>3SUrUF!LPe3*F%gK z(uexnt^0#sk9qpuqlv)zT5HGUaoaEJ^_zMg@Af_9jb5*K`5EGIJX)RB{QIBd;k;MZ z?bEt`Ch%hUY3D(GV%>Tz>&lYNunGRG8b9(7zlWjof<7vjJ~p)p{v6l!Y%1U1L;7O*V)aS< z+<63SR@9jgzZfs1PjtHV8u_VyP2k1yGo(*=?sI<0)2a0QjPXMH#Gj}2Im~lq-g^<~ z*1ExY)7o*_PpfaQA#S`*p)J?LLi40fyn1MtXLU=wn!6CUIlcZ{lCF>A(b8}9qwO#7 z>iT){xU4pYJL`YbNXCeL+irhy|2yuwEmR5nAH5A z(a!zIcK!(R=cjr6LcCb~m^Yk<_;X617v%FqHNDT2;LjM( z>ZU$$c^=!$*5BqJue%}AOo=^NF&(9c-dBbtZvGXYCYvw%Ze>Kh{$Bs)L zE$1QQ{)hFxDMNobKSTOrJUb8SlYHsh4BqHE_}9zN5U*-IP@m}HeZj8u{EYEJ`h@qE z-j{Iejm*Q-_kMXnex1y?(Asg?PpfaQVSm`4p}kw}R`9%fxKemlH_kT?Vn3hM=WJXZ zieB&hGo+8>F;C&A-0wWj{sM0{uRmGGf1KA}h)0ecmvOWGMV?;w&8&r9?|ylRS2aJ$ z6WgcjT9-aw=kb=RLpX11tqvZy{gR*KhmnsH%4z6%niszikK@s=j1zgZ><3Suw|4h= zByL=IB&*brzwV5x}S@7r5M>1)<^A@r0JWcrC`gx`BZ2Uyez$1uPAJ6C0^>IArjg061 zp&xmi{p3j<%ihW^L&jaiFIj6Utg?Qw!E_7P0dGR>= z$&>b|UPrt7lREifTpu}h+?YOzpSx~HJ?isJ-t%rDeO2Q}eZo1jAM2XCZy;@t_rW&f zN1X+Kt{ceU&%LnjISSG1)!&uEvvC!_H3pz(oYxQfoNIzUj>o*A-{?mkXFqwUKfO=q z@>BfzQl6h7f0ARzjq!w^@pG1Y!SkLk5AmwzC-q@H)9a+^oh9Cl>JZM`+EjhQ+430b zV0fSI^J`(fdI<5zpf`J)QF;e|`|> zQI{Pzq)*0uX#nd5J6^{VICgKlCTs2lTouy&kIagyZ}veq4MYleT*vQ*?Xv-mboR{W(CMes1jf^Fwez zYu~1yXzjS{r`5N^y!}`v|8eJMVr%+*rl;>p;aT0XzxN%*`ZdP={XVo4{K@fH7w}K` zdWAfBovR0~XOcgY`T7#_CpqL=o{gW>!}g;YUdQLz?@HDe^ttZ^NP8xoSGE4A z4;-xrauebg<(}wOl zq0jq%f2Z^IAbA^a0rUs^nVze>+{zn`x!A%BuXuH`XL z*)RFoY-pY9K2+2ACFEy}XXinG;=E;ta@C&+yjZ=5^od@VpQ5+^P;&jnc=PJ()!!}X zD)VmA0C2w6+Hvi1#!u|PXEL#)hhS^BLEj(qd|7Co_;mqzSf^E=`(B6WtS!U&+H~jj zg?O^=c==h?pUflH3v|?KA)Xxfo@WU0m?w7JP`+TleO&#xy>@0eZFU=2YS8VV+r{)#^utk-0z}D+$Qusy~q2oOd#LaIDd{uzv|0c*EKxOeyNAop3r>J z=N8@l0@3@OJYPaQa_qRA4{d=rb}U0%?_1BTD#fo;JkfFH3Gm$KM~LI?XnX5fyF9C# z{sf25QQVfytD2v}bN&?GEYAni^Q9_}{**l3ulpaLn@aF!j5ohO^>-)gn@av%Uc(>r zjh5pvPxWQNdXRcMe;n(y`&|^#>FPo3xAS@k@yH?9@;EQ{OZ=7(BmZB~&o|xqg!tL> zXGeKUp3H>ro`1nSZUKJ*c;e4B$H8gc54jobKV7~*t)Ce`(dB-(6?J$MYyvNqpP~2( z&z%nnZ&E)$O7Le}JjQRD=bsyr<5!heh+j9)2iGLWFUAYSPq^-PAmG!@_GpIs;)1xZryqUJoh`JqIcp(sXt0Xi3eC;F~E34gY(|E@vX`+d`n>RYZn zcRhynOn+b1yWdll$GT)b)$|YhQJE)P@A@dj52Mc*v->HVIB{H)3=4Ae+l^+;;}!J zpOz)0=TyZ)(6+WcOC(5Mez8$-mmEQMW0(B-ql|MkM~>0wGJ1G_wO%*=k`CzlYYHE znphxx(_aDaoUWhT2hVn%SG)N4!dJm_oQusGxeLG(|NZl8sDoaOoAayG$u~N$PsV+9 k9K18ib9p0q`}ND`aXk>f@t*_FLXXQ!srO0!y?e*|AAmBGQ2+n{ literal 0 HcmV?d00001 diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[3-prtdry_stop].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[3-prtdry_stop].npy new file mode 100644 index 0000000000000000000000000000000000000000..ff9b202367552d3d592d4ac0d0af7a7f76575104 GIT binary patch literal 1504 zcmb`F!Aiq09L8In2z!d$GNA{($e12-+HGSf9zXOQCNSosg>$y_Z7Kg3#0Tm>v|w2HHpztRKn> z2T}zSN*X*-@fRFynp)sQ8*Uyo&XnM~sZnN>sBG#MulyNZ8@T6o0XSC-Zu9f)+Cr`x zWfZTMjVdkmgiLv*h0XOHm%{qd(P3dP;f?(NGt4vh#Ht_9U>+4J6 oOXwT;vxYpO2PCB1W5<>7fA;mo@um0e{e_RQWWK(>D87!qZ@Hm+2LJ#7 literal 0 HcmV?d00001 diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy new file mode 100644 index 0000000000000000000000000000000000000000..5a921e5f6f30448d3e8ba37466e6c13a9f0cecb5 GIT binary patch literal 14160 zcmcgyZA?~G9KWRY!WVNg(=tQH0vxeSN+ghjA&9)h13Z8)EGQn5NxVD{HN|wS#V6OS z)t0%|aMtEpvs!DpT1~0e)^z2nQ7ui%%+!^wX`A_>^ZbA3ulF8&yyrYzxF0y2bMJG1 z_xFE0zjN+=u_jzyyEesnDrIYCdnEQ)G_xRmeOBhC-bgevCq1*UGe0vsCtVK2`g=VC zolmrPc?Wu<+7;2>e(x1sZCktp13ldl??|k_t^diGcPtudi}iGS#`>d?E;))Hcc^Cs z8uvES54ZsSrY+x>ufKr@Jni{n>L!85_v8EXZS_8n zds@b8trs+Iz|Rkl&wFqjjh_YHG<74LkKXsFQB+K6+nW z=4}Qa5Hw&BpMH4P=Vy;MLp`z1NAG)ysh$vm{9&cA+tmj>q-FQHNwvl=^t9h|?)boh7K={M96g_s>`H zZd2FmbxQB^eG`pel(W%DOPe5n4gFVKcVNNL*%z7x{qB2nJe>#Eah|XEIVOt7ILwFd z%YANn{H{2$*A~N$GnJlq8qjyu>qPdn#Ur;Z69V|lMckaB{IEDD(_zR}uYm=YI zUP1~l7je>~^3PSTgW^a2OQ{!qRKT)<%k|7vuYUB_f z>@SCoG|aP4H)B3>)$5?}sE&QN4)^akRgAnauV4-zTYiqM zafe-6ESAQP_f$O-wcu+|&usBXr`#kECodHWJVy=U#h{+q;=L|sz0qkEs-6`AZ?F=0 z2KCGqkK(_-Lh)LGwzU#`4eE<69_g^VL&*1ciWP1J)(?|_M|JsTRi96*d1QPaXoE zoZ$C=0kB)q}G3O{cftyZ1Gs1fnP4_%4c=P`pg!O^%?LoV8^fOjrEx= z9_usWd32EgZiBHtv&Cb5hV{`r;IZB{s4up7tk1wJ4S9YvjJjq}Uu^MMpFt1K1V5>- zNpBnEyDc8;Gn$t=9NoZ$C)@Q&k6@>p|$XK7*;;}vhzcSRbXBD1-KC{L1=(FK60o;(W zKC{JReFl402;w_CWUS9@@mQaMS0(159ks^#%odOJ8P@Mru%}!pt>Zc_IAivS+e* zIzO6ccXXhhtSJYt&kJBzm!3$3=Vwpdf6Nd2JZj%g?ibL!YrxMBFA+b|p*3~LzxUUS z<9B;N{-gP5&PMR5UkyLX_l^`lE_)`6r}M-7-if@syc)RwEJJ-U!1J@G?my;7=923$ z%Kx`2jQRQDCE`c%F|8i+yWDT0`sw0#M?n6g``*&2{8pjXnd9x5ES}Df^6#Tf;4`lh zdhk*n^3i~wpFMT|F+bE%^&OF8ZcjhFMEq!8KU5E0eX496zxe_Ak91N&7yKvd3hAfK zdC(Oudj^F^dA+X-`%A~Q_a3hIjmhHa_M~;@fo`==Pz^hNxde7K!1ME;9xrT9(o171 zUgW$){xqwV z%)xr|gg-yq_+_1yP_8;L&Ji?XxB2-5W{O5<4 z$bX~bs%<@O1wXzptZ+$KiU!=g3054?O12&o5r|`H1~T>*;S* ze9HAL<)eZBZ1LECRA(ZM@T1i8G#?t+(-x2YN4WeQ*xmWS4FE5Z|7e~0MI-9yw)?O6 zPv3uFe{idxr{|V`_&m2fPX4>B@NAr0%K z(ctlD{r91=r;JbR)8DG>8ay7&Lx&Zfe9uAq2N$MeoqHsCJkCqDH`_B29@U!<{PFm{ zx}UhuE#q-MP!H8do*BmZV2j7`NqVeVy@zwsjPoH89*-x=2OW>|;eW)M&rbjV literal 0 HcmV?d00001 diff --git a/autotest/__snapshots__/test_prt_voronoi1/test_mf6model[0-prtvor1l2r].npy b/autotest/__snapshots__/test_prt_voronoi1/test_mf6model[0-prtvor1l2r].npy index e52948316f1b0df5420c5130a5102a1955525df6..ff1e38b14d337b44cb7610524080d0a767ab7cee 100644 GIT binary patch literal 24384 zcmb`Pe|(i?6~OP968=DEz(7E~in5_ouU{6b>0(onTFm^)K)h3E~u)m^(tAmV19M)XX+a&8y3~~E~~4qtY0v{ zS6M?{buE-S&}T` zJmJW5a_#hEG2w}n@?&#{Fh86t_2q{0!@sHL zb~(<0W+!7;vJ7wMzFd*mtL0+(tznja;;M_$cW;ORTBNyYh=CZB; z<%Sc5n^JR6H*&F^_g~I_)sMMtx2U;AMlPyVT*2IvrGIZ#bI%9r3%x^^U<2271dDML4PBr(;2;3cKaUAXEI127kbB7waI1U|v`3B|Zze)VgQgeqH zxyU;S#$QrD%-mFknmgRcMSEvLp2YDXE#pG0Qq4Wf$i-oI{W!LFAdd?zdhXdCxAZv1 zg+w{qTfEWDJIZIP?LEiHrE&jCo^PD8(fQXKdTu-dcPq@V$e;h$?Cf1p6}alc*D2>3 zxu}-CklR`2vGF;9t1g^-o_jD>z*5`4=*E3veuXl9!bUg0ds)vNVdSD(%T?UJr?Wp- zR0poQaC^^p4<5I@7)P&9;Q8}V=El-$?gbuqw>JRDrAU69#8LLP=(!_}T#BOv=MCsz zswTMG$yj(hml(O|PNQL6hB7Vw+^Of58oAiNh39eq%DOu{S8eYoBNx@WVZBf5GU3K+ z)Z7b=Tx{NA$X950Qu5U70Kdx zYVO5GE{*#TN7TQ4*xs#r?j;ep+oXRBBz|kv_Fih_()f4r5x=5pI`!Oyk&Aw{8s>Gx z|BuWk?^WA-na8bH&fAgqMk#au#oYLOHTQBO7u8b!dZA6vO&Yn_&Ie(=K=C8%g^C4% zt1f(8xWdS#ejU!^(OcZkv5cB~rICwj!7%3jj=3p4ceIg<`C`|Z%}(B28@GLO;*nx zZ{(66!?;i5FU&Wnnsz;Rg2$cZIj&raaq`9)WBwh=Z$NPT(q~= zf1gsJ=T0(m(cUENZ&1F1{SB(7RnNW7z zk7;j8&z)@K9T(tKnf9}?Erx>}%jNf2!;HnFcqZ^D|w0FQz_A8m^H0ik? zHFD8^$HKas{P<5Y{_43O^SCp;fRRg){5XvZBmDR5i6w!%AKuP4M&Kq#aDM&^+nd#M zZ!&Umyx$A^D%j7-KQ}x1Gd=faBNzL1y5GO?rGcw1yq!O8M`7|rKkTp;cs<`y*u?tXYXf5OPc_N_W-(IFfZ)i=O*QBNutM0+;;xV6$7VysPI<^SJZ8fRRg){5Xlg&(Fc~z}*jzqv>2x z0ZT0x?cH=X+xvU&-zGix4kH))_f&t}@6>Z=7`d3Q77yb1-N*JORs^oP@OGYQmzWt!}}xTC*?hRho0MD zCLx3gZaR$vAb+&!wyQI-I-MDIPbM)<=*h$&cmy zMoQ0JV&r1G7Q%TI*S>7;QayL6k&AX7l5;Y0-mFE>ZS=Tr_Td5+g`VtrH_;=0seJdU zz}*k;-(~K>T>l;qw*}@o)W3h_d0n@jyF3E72~fF+y{(YZ0ACMAF1DifvYaupI3R@-$V(FxX3R)WP;9 z^xV}ExCaKXKTH2k({mq;!0mwZNR+RB%l2mV+=q-@RBM;{q?||EqUSywf!lTg_wV1? z-gZ5AjggD~ydCCs7*A>WKIJ_bZ}5 z+&4W3l1q{Ncnmq^kT-}2d6M$S+j{O3MlQB<7o10;IFj}6z=s1@UAR9#Y2>0>6YL`q z7tUXioS^46Mc`)P+%@rFeMEA$p8F{y7u&hY-$z=j=YHDAMf)bhcWvZP19+X%s^_lv zxc&KBg!r*YDMxPK@f^QHn7dEUeagtib{;XB=U4ri8(R~&>cai`Ge$0|1^#|rsh<0^ zk&F3aCd3it4LR4B(sMr>f!opJeHpfUa*3Y%jFF3K?ZSoglf->a&)wj0mlSXT3-M!- zau@yiX~|d6&J<6t=(*2E;EpWe_=WfER86;@`<#)BYE3ZTp!S9Dj>+CZYXetZcpN?N z9?Y%xwe5_&D%b}me}?l{#2v5ael7yHUd9EOPZGCU&uunx(az`Ld=ByCJKV?h+>H^q zGp}LpF#ay}MLoC0EU>bajcaP^XRU?=B6V`*2CuJY$89jHKk&AwHXc%+5c)ro5 z=YGxO{*tdn$h=CVl%s6*+xZ7dgsqlK?F{?RBxV1(=&`_67w*sBh`^nC8n^Rn zyw5O3&uxvs-3<0p-V%E&^xSVo;O_AIcdefLEh87>XAH$0Sk#Ak#ZONcQ1T@ME)#)f3j20{h^VIYOQ|! z_Io05)rF4>KXMP|#s#gt*uRg+c_jIJstG;!*jZ=fwX5BBELw literal 24384 zcmb`Pf0Wc!6~Jdj5fx=+bzOdlBd80CFgk)jn>@`>{6S+1fuioZ?8uK5hMh$LQAbu0 zb=5TsSwq>P2#abBBTTC`i|Cx~98Q*DQSE6vL}j&QdBEt$z4zTa?7m~?d%0iu=keZr z@B8fi@$UQaeKXI@n0b9>PD$eKl6x!X)#Ym&DyJlG8C|(3SJzNEHd$G0#fsTiquvXKf4vw&Ms+aY^Yg$`;vxvV23HmX*KovI%v<& zugTRxee}ht)L5{~eaV>rPAp0If7+1-x* z^Ue#tOibW1+Rpt+Fp&^_kA&MXlDU5o?i4lmcq2Cr|MHdH4&$r!pTeE0<_<7&f85#G zxoS+eLw@?-!o5PxJs}2nZ>rm6@!zpLehXKsxhERADAqE8{X0&I? z-tAz#rr`VZUL z-Ko~M(#S>L)^R)^m2kf5(sM`0;0}a*g?>)=V{T!-THi57E{f$Y3i5Muftq{1k&9}+ z8SX(aa!lBgSevc193&=tc7ZQ#~Qg5$LDgslKkAG=U!;!qMf(F z{)qGhE|u#Rsr4NfgF6bi$WQ-`xzp;^-0?;(=Bpoxo&V0A9B$ z+?$p0jNGy_dH($k=c{3P=evuxGrzzr#+ zcj~#9d&3+t7wwz}|I&OU{;gzWo+l4vVZ$8bN=&aJAafrD&bPgrMz(luiGK6NIawGPBU`R zzx9wOF}~80C)@PgtBhO}%k)}5Q#U<&?$t&v=?UiyRLc8bcenH1Mf-QUk&E$Mo@W0J z;C*+Co_mdvi*a>uGV5C==L~x8wE;J!oEIoBPvX2$%G}gF&UY7W=j)7I${TRbK;r$G zThMc_H*#?t&4qOx#TBgUX#Vfib7#ch9)fip>YaU$^VNvuUjIb2oj+#eqJQ6p_hqD~ zyf4e?xfvrDd3!HqI}h6E=l|FB+#8HsH1AHwb{@Xb+q-;)+Rh)3!7aav*Ey9NeSFpG zxi`k(PMg8pG>@ZpJ@=-7dz&(Epq=;4u~^GqWb<7(CA99L&CH>2myGICL@ zwAVQybdB#Hl;*Bd>0?rb9$#X8_zh4SRv%ysuV-(7URn&TfVuAjAW zh52fQ>=)pD3hJMIQqR5F$i=ui1p8m&75Be-?kz?x&A;+~@BrI6waWSKqU}7_$VELT z%KOmb`vpDsRwEbV>CK58&!_PC-L2=|7H~_0je%69Bd22C%X;p-fGhWOx%vOz_r`eI3Hx8N@6UMO zUGjkQa?y6KGjfqP8T2J_Ip4_WxmhC@#d_piMdraaJ@-x{7wx=R&iUouW{;jb-^j%{ znhX1ViYJNZWWo9FqW!zT$i;C~2KN}TO#g|;(banHLL-;rIlM2L({mRE+_GQ-kgAmY zcmm^jNO(VXy`EcdMMlOma!*$LiJvSF{2L}NpRVn##Dq7&aDC(U(aigP#YCU%;cT~cqmW#Ol zu-_;D%6*l!dTv7u?wW8N+pg#4V{j+-T5q8K*|+rEMk5!;Q5EbLa2%)q%;PAr+WGFH z^Ybzz7jb!6*U7o)Xg&9CBNz4E8m`+jdhR_&F1mA{oac<-akN~|U2f!}Sn^sPzvGzO zrsu8*xc&NZ0~b=2Qa&YpflKuib9b%wbw!?BNrH#JU#8`vSm83(cPMjvRy%A@D7g3f z2aEgtBjF|{@Hmot<`qrOcNaa5R>k19!#O{Wr}PP|uhVnyi@~ji{UFWza_&7}&%Hke zw-nCXi7Vf0YteHbh{5d$&!2bbxrG?q9NgO^eG}~8E%rbRZ4z5fqHI$^)tmYeBX-7@p^8Pk&9woa36`p;XV@Cd$ykY5D6a6S6Y2BPp*Oc zgtX3im$?N!_hBO!#a#Gaa*Lk3#>hpx4)1l2MBVJxb00BsF^FHK{onyTcWuBO z5(KVPrR2w{$UuHZKWF7VZ1Q2}KaU4NE z>ztnZxRFcxh2Ljiuje)!xj1g;K)#}UA?NKKdhQcOE{YYxd11ev`=pVJ`o0M7_etNu zJRg;=alX6gcz(*rMX}vi@;Ew&xs&wVmVkRQ*QJngE2Xq!-e`h#ESBl<%&pdQpEhzS zu3*1F+%uWGR?q#UkxM*S$I^Tx>)3WZ_ZcG>#X7=!6npgCPZ_zGH@1h|#3Rmk7ah-^ zHgahkh1}74?m8nE$L$8-l7HoUiy1xlStA$Y>Y0$cT+dx^98ymWi0# zC0|I%%?Y&cE8+L8cImm##o$iwwVtGIdi31S7`Yf%TjBhQ@`s$qR;+ctyXbMW!N^6i zQkZ|y&smv&ou2!=k&C!#;re;Lp8JB4OL5gJKhrpVLeJf3OlE{gSp=Y#_ub-ugkc>bJ!uy~H7)fdOnK*(39cky?7 zQhM&D7~HWJ@Vf4OJ_nwy=YHPE#rT~bzQ-x(xnD4HF@Cp&_kp+Qxtjy-v&wizu@rD= z{79bMt>=Ex$VEF(hWB)oFXZ0!0X=t14DOgAtS@|Tg5rAQW6pOMJ&s;9a*?wS?(@?= z1-?&5dQaDLzhva1SQ75llYf87+@*T%myKN1Z-v+yz86D!uh(<82HZEfE`_Xbq?C5# z9)x=Zq_6yLK!=|D6(bkL+RNF$a(`pLp8M4p+|*#^{)FTCl*gU#E;^pSX5^B82Qjyk zxfAu=b|V-4nuq&26i=WpjiVYp_a!42#hlz1mG3=0rssY=26s2)XUbFZedl&PcU!>S zbPP9eA>&6%`4sK@j?BOCdlb~ZN6&qkJCbQg%SEwFc+Q__cD}plc>aceu(;0Aa#7#m z@IHm|2E5;=xEiD9ZZ~q#ujS$IQDpSoZyLF%@2YSgx?IoQVdSD%7o4|~edTu%+w|O5 z0&a^07`ICtOUca%^y?xxRV~C7iOc*uThINzk&Aw<8qNNN@6%B?1wHo% zMlOmK;GB@e<##Z)=(#&d@NhhbaeEwp9=4;t8x}BE{JUGv{o&E%9?)}l9ZhcXDd)TE Gn)*M7HDdz+ diff --git a/autotest/__snapshots__/test_prt_voronoi1/test_mf6model[1-prtvor1welp].npy b/autotest/__snapshots__/test_prt_voronoi1/test_mf6model[1-prtvor1welp].npy index 3a4f3890ddd06160dc10cbd4d9ab0cf8ababef66..c3a6a65497d530ff1045552acb9c263b398c8b72 100644 GIT binary patch literal 24384 zcmb`Odt6mj7RRreju}m)Wztw(M-xpllfp+jp0Y4ai^fOMLq*_*iikWk4JRh?QIat_ zmZ*WRAt@|P%j~IRW%jY62n_2}<^ywT%BjYw(Q?k(>$f4!r{*r!iGOxybJzZ|*7~ir z&N1)jnKd;dwee@NbCPm$vm5Var6y%h&1kfjla-q8 z?y_+Re>L;%)!oLxUh$t|uV&T%?oQY3(A2D)tfY)_Q?pXsIP?$3C#7eny8GE_Ntvnc zb+_KVdd0YLnH}``)23N~f28ZIRxrTQ!m~F1t-ZR!;r}`|h}Zpt!tt8J)PBw%w*&a$ zxz?Xwhuc-mJ?&TZx0jgPQpsIjOh4Zbilv{;S2%X%7;`M(HneZgkl1LCme|E#|gXa*e++w>j+hHQ-$)=AN$Pntbk@fV^IZICm3s+jzKX zHsRQ}y0i6$?SSFVNdxbC2_{@0m)j z(Nn}U=4Ic+x{i|T`L}+!>vB2HuO}eSc0Ppf8)y09Zn=m3+;%hSSl%}xm0V+XYd+#T z65J~tJ}$h!XM4E0o;tGiV^@4Uz`B})IFEv!FBfy$D!Hb9hs}S73_!`-)p<6Cwl?!uQi;`>1eyrd+>X?T7%I{+@R&q_B3_i*JRe~EOuivgpt}%<< z$o0E+BJ%_C`i=H*r+9Q_`;=Yr@c=(&Z{qw;pM*He{q3gY8nc>B%o{qH^E_NWkGd`rk1kPiO&t|%VQ$wc?01yh-yR7yR&ZKi0$T;mPm+t>3qHaGY<>W`3RgIpqo^m(23x zcQUtGPLuO{rIJf=t}A_&xtnr$ZNhzBoQIp}iL>pGcE!g7R7csBud%L@bHSD0=k!r> zsg7#f>~f5GHv zK74=17nt?_F0%sy0L!y0wgV==efBZu+I;pmB%cdx7^vq)cywjm_i|0&c=rJNIchHJ zdS?j7h4=SrIAGnBa?Sn5y`L~QW-hN4a_%)sF8Q0d?L+3B1Ap7yCH8l)l1u)6bFh|m z^+N&k^W@xXm0U9iqd!3X6kwHeuTyfVek*$JWj`wm*zcif;yA}Exn#B~p@#i@9&xUa zbB8Fol;3uHtGRw2gTGPs`8#|aU9aSlo})%obDY-~IOdp*a_&$im-3t5zMA9wcmex8 zcoOpkbAOl*Zegox*3)aq?;bh#20z^Mz&(sOJCnu!4p(x?EI)+&2GOrt$hjkwTF zDscZJ=iZ{^l38Z|I*U-!lw2~~e-q-o z61o~C=iaI0Qaa-=MVwu(IVMHU9joM0oCilB&MUypk#ol>xnwr}^Fut3mZOg5$hqT{ zT#9qt8;2ZY?sDlmv{=qfQgX>Ga>F6kb1`(aO3t0&gL`Zx{9T7QKPBfTE4ier7K@PI zr=hDC<=hk{m-=MCnTL2T>_FexE$5~xxzr~+rW|5DRYF%aa&Ab;HToS5?#IaQK{@v> zAKX?qAkP1Uo{!17X+F52YoO~ zVc$L{*cICW%5U82I>%gpUBdn2h6Q{rcwU%{gDg(UC9|T1xWD?og#8^W=cX&U)Hl}7 zuH(57DCIhyCg)C3aw)%wxwwyQUCR6dIX6SerTjjejyRuH>X>7e$+=UNT=KU(6>)A~ z%DQ?&&dv0}JvtiP3&DLw&YkArPO!qpzO5g-V%sbmhakV_gIgu%X8GU->~qSwrCdjK za&ESgOLa83C!SNzMjd@G=jJH6GzUNGg8W9Hjw0@N__*+WGS>&UFbZ6^eo5cO7s$EO zm0YUdu4h75r@`M|a_$Tt-1r|4v!ACh7Y56@cPqJM7CeUMg%g+yW8~bKN-pVn=fT6= zmyeh5nkMJYQgTVpCu{LM=+0krOrD&3kCJQV(L3<>OT>ARoSUcQ8h>{ozyHF&gI+D? z&Q@|Mzu)YDp1*;fH_N$y@^I6HbAkLV-0I=(l5_7>aw)&NH#Fqe^s8DqH($x6e)ak~ z==mS0qc7y#Kl|XuFGrjYBfqEQ-20SVGK+i&dj1r3)b;^~j|)GK=I{sW-xEr8MRAUw zhd6(PICqhA3w&_P-1lb2`~$?fkDNPK$tC?Boq;&-MVyDpx$~43GWGkdocoZHOJR+_8uO?MeWUe4hmQ+C7mECF`@!Ez=qgIi z{fm-IX5)J!zk4tjddRtpm0Y8vp78f=_&ZR}U83ZYp1b#eziyn3t`g+jzxv=NMq>`T zx-$MJ%ehOHT=FlbGxS`IxiDSMeOSpQ|JDZK?}v!<0y*~)CD)j_T(j08&MV~HWj?sa zE&}%${H>64i#^=-&9Q;S&I7yJ-%@^ucS4-q=We6t?Q-sNC6~;KTz{!whE8JszAxvl zP;$*Y>Igmm06o{qxhs`ia&vkV{QVJe{!Y$)RLLc?gbwi6zRx+M(Bb34&!ZClVEuDS zQ*@Qr9{#%T+sXgXg>r7Gl1qJKPCNM94DahMmvdJsx#Zrnk+{EVj`wxf%DJnRT=I8` z`#n49y0``2Ki?wfu2FJ{cY161>%N~P|3j1H+_g$BVai&fPqspw^W@xhKDfu5!(W$c zjwzCJ|EA+&CLuor(L&x8&Ssm0aR=w%jvtU-yZeyH&|0 zf8(v5?fu5ra_)0}xVCS&_gAd1MGhYqelBcNa>=adBh-<5A8U^3Am=`B3>x$-)fJgJ z`%p*jbBfVd4>|V*C70q_`vLUqJ`b8W50rCXRB}mIaeEu?ugv^SkaJ&Ba>*>R8uRzR zn7_$#?#oIp<@4Zs(A9T{bFQ5GijqtD9QrPF^>5tgERb`zE4h^4>F+>S_Ivgfa_$Zv z+%c66^VjHkgPi+!57+yCtR;T$-KXt@D&*IFzhU&eUC!O9O|-nvd44D2{e1e|hv)77b3EtXbH4ZaJ{-qgn-~wYs;*xL|f_R*!)(V|w#VlF@7|Co7ol zj`HL*d~4?0-gGm;_V}yV-ppU$yTdg*ZAw;7R!YX5Q?kanIt+{%nv$L!bjP#Dr(_1* ze$Re=`oy}r+!y2X)v;B8U()uqc9?*sji+t?uI*jo@V|~OY33E>JdsbYOC9dViMj2S_8MOFAScee&=&U-b5Brm z4X*flws({Ro}Al3$u+vxCxSNwI`>6CbUcSCxn`P^1pi*+94MuiJgz!=xX;)?!?^9* ze8X_S`1krC^TMx#pXJ;W{cx-90`~@PhsonPT*)=lGABXb8y#LI9&KpB*{N)AwEP@< zmXd39eqvzQp=2 zn~40CpC@}Nxh9YHzRdRJPJ*5CdeuwGHTISlBc9V+%)LO#CI1$;dzE?Uo)-;zll^=4G`711&tv;5xiqb0>syYQUSqlP^VI+|vFZF=;?eo{;+b@A z*uS6Ss^A`+bL9RVsN_|RsK_7+V?{>sm>gFM_CkG(dYY>SU`v*R7a^9*q1=c~a= zF2&WxZDp+Qu9>JC^8M&?4|li6-rw^t)s63y4l?(YS+G~WAN|?G{m{d;`q~yB519Nd zI>fvW?{j!u&b`9J-6G_Y&egkWSl`wUvfa`0`iw97>pgUphx>av)7M>k>@dgG(z(o! zACKii-@C7*1>Bc?zR?nHAoml;OddCv`RQ`*Rc4}&>-4T>{|=gqx*_KdQF2Wlb^Czz z?F{aIIrnNG++o4@nOi@Hb#9b%uTgSoT2aP(%qyA0I`^{o;Lve3RLQ0MEx)yrxzB)` zFXvwCgIoOeA?B{0gL*a5;pIZxdz}w%!+=Aq@1t{=S0?9PujJA`*uJcSdF#PVw)gMQ z_THf6n)Br93f6ZgxX0z(8+dE9jHTLE~=Wmd|;r8Ainmb&{ zHF>li_MSAC^=rzay;aFIan%4BIW0rYE z&P`Nusa{Q6@;Tf2B)EI!+`lNfr1R0Z&)H7Pt&np^D7mC_!~S~Kc`fX%m2*cbxil?e zUOnqvw37J^a_((DxVbmfJ7&7&wwcF$K|J5C%z>KrrOa{I}-W6Z>+bNkik+s4)@+uCu8=Qj7(aXhVoy;sS(Nj|vc zW1#cum0T}xk#p}*a%o!3#n9R1QXMHxmUI8=hZ|PMcCJF*$dq&MRC39`k;iJ;&L`mC z`{mqZC6{zfDzA0SbeBu{oWEGkP4U64ezBH$`LMS@&K;}dlFkjwp|e{j%`zM1+*Bo( z>ea}HYPn8rg1y`1+;K`S`7>`S?A?L-{HB~6RC3AQeW|rvuMR-x3OP4T$u;}YXvEcr zs2eqM?p;bQ*&9C`{{0epbWF}2ujEoZmtPM5eh2^lAm>g{a>?Em7eVI}@b^$19(H)S z(D$;5N-ou_#Q7`Jh3TMP%NPWHMM z_MTqIadkk>%~Enn=g2{@Hwtl8Bj;u-x#nEc7u?S9Z-bnhqvVo%yPw|_&t!kx33(1L z7rIVP^}*d74gZEAu1=G4r}^O4hvPoj4*rdmbMH}d$-iGT;+_&f-EidG={~s8U*n$A zh<)K!Id_JVYtHST;hyq8>4AQ zxfEBCnc%wiQeLHvmUACga>>7gCxZI{bWWFZ=PSA7-%&yMw;FkrE9d4ZxyIfU*!wZK z%jMh!KDZm*=ds4GHKHcAPj#!=bZl zFU5D<+j8z=C70rAz)*1AxT1dkM{@2GC6~-=b0xU8ULBWnA5n72zfpt1t$@yL7CO9K z=zU?S5ANo8aLZ7yqU791{czpi7fc+NBA$E7xsNHire0l$bIt+Ojlpv6G9{O$mGngZ zTAgo_bC>(zMw}0w401;-P~@-BYWRJDCe$Ha>>qyu88L!p>wUAyIRSmxC-k6dw)V6eJAIxQF2WlodSE^ z=ianGrFC57@N%K|qXH$DFx?|yuYHcxP0lU!!Hqf@_O`O#Ic>X*>&V5qJCI61>0DIfw`Q&&xcb$?;`Makr>~)`uQa|oqId{F1 zOVf&e#_xmf?>MA;{$e@zDL>pF@jG1y=v*Y{KCR@M{B4B2F4yQ>Eaz_UaF+_7n-Mqa zd;I?E>P-E#eRA$IKDd?N;CFU+Uoifyl5?L`a;a`qe+8ZGbHbx??sG~m`7`wh;_4{& zg`ecy=apR2b<|<}-u)T&--yKyFBf`0+Nk83`uquWcJJM0ne*h_7kqHHe}s5;xh9_D z<=nq3xh8F@q4TG(_c}TEMJ1PX4*wW=^f}IxW8~aTez-O8ul-$Pik$nBl1ncMx&)f7Iu(a_$Z#mvqi6 zLtMFk=Vrb~j-0zw$)&j3_YUIfztDNUoLi#g(zJvFh%5K+87QyfR>`@$d~j=AdyT&C z{!6;&zaZ!C_Hez=soLRv@3^tgjVpc|b)!Vi-J|5vw6b@bxaND5%ei}%T#Boazu_94 z>*U-wlw7iJb|jvQy3dh}pBv@eeM&BIhIMY@lKp9&9&vcN(EH$Jan9-6#Qg=2A;Y}@ diff --git a/autotest/__snapshots__/test_prt_voronoi1/test_mf6model[2-prtvor1weli].npy b/autotest/__snapshots__/test_prt_voronoi1/test_mf6model[2-prtvor1weli].npy index a0e91bdb042c158ab85ae7f65c29a314106d736a..9ec69de0ee02a177b877f280766e8e09e2c7210f 100644 GIT binary patch literal 24384 zcmb{2dt6mj76NL zlu{OHnPnD5Sy`0Xgp+2PmNjOEpk)u!$;oQY+UvJ9&gWBi-F4%i?LFRozO2XI=bU?4 zJmBVQ2M%fED`+$$E-gJUHZQJ6^so+b6LQn@;^L#@dS!Nx>lh#H9th;;HaL(uF)gd% zf!sXLJMwb#8@?keby~v%W7D%dCz+m;o!;<0f&A3`$$^H)^3qcSIoS=4<>#emxkq`j zAHOyG_SIA~!M@^O#lD(;;{QgtVf%CP^7B%&N9E+Dxq0Xj-7hsOknSE2WTfV%yRSQR z?b0RQ&C9IlFn^ji^6`&!oF0J*7#e%VhW~3{b#wS%-}K>ePoH4kDVS>C7oOV|esI^u z^T)WcV)xWPG~O;^w~5kST1G#|`XtcLx^9jg`E)k=@Hutz0w%|B! z`_8U69QK)b0~zcMxf))3jKTZrbezccmDDwDv~dEvyRJdJ9mQ@7rE8{b$mTr9UJFlN z-o5?-@$h^zQ{T(9hTN=$G#st;2Dz zEaZ5P-GF`aa^yRByk~mc$)0s&*N+|Xb|0_TdF+(m#PQZ%CU(zKx~9IZZ%17RVc(Kp zkIq)QyieZA-s-`e?=ty1YNK?`IvTJD^}Pl8y%b~c`u<7jnmn&8VaFMY>tVduja0fE z@4cMotYOIW#bWmykK5bR&%bqDF#Y`SVs?v%!|fz?+bUfXZ_54HH%4In%J<1ArEBWD zyOiVGFcR^~_p3iET{A6v3CFo@6#Ioec=QGD=W~^=S-(w}qQ0Zi&ppNNc}mwzD}9LF z;i){PB#GU29(RJLpY6VEN4(u<`f6x7_f_nekn8&wrE8{bdzjrFaHFnp_;kVRd%n^& z>vzoxcFM={SeU@0ps%8R{A2Iyq8RWSEA(Cg;=F)^1N#e$GbWMeI>u2j8nQat#I@+oZrk$_EY3;2ah|^ z(^qzXv?Jc`qj=ZltYt5Le8};3RJx|V=RWJ0>EDcJzot9C7VK9SD_zsiBi5n56WA-1 zUk5Kyy0kw>t$mKYk}U4)b@JNv-alD`BaK5AD&)r>=u35jgz0CPd%f{ma$%6U$n^Fp{$^8Go{<6h#a@9*s!u}N>RdwKx< zEPpP$Lg|`y)VG@RTNFT^<@e9sm9CjK@=abxE%M>Y*Kd;2HSy-W#cpvv=Q}~Zj(T|9 zA)fl$>wz8dcAvQ(-M5kJ+G;Y~IEVi~`1>I|2DgKP3E{$%H?(}0-NY#&*B4)4`s{P4 z)TMX>y=&O3o)Yppau8RjF5Tx;wBO9r2TtX_s*&%{S1Mg|pVQ+#j`#3X#4G=t@G6fx zO1OT>9nx(J*SB;UT>0m*$ziy=KH0+g?KYkB?UO$b{!Qu9v<;7LW%q;WoNxF0SUz3w z{i-j&@kBp=WET#O+n3&TM0u{td7twsZN9B2M) z^p*U&e~rgoVmB@vx8Zx^rG0Xnvy;8#TOEG8{B!<)C|xtH^DfSFvqIEYzK&9quGuF? zf5hu358o4pa36cE#~mBWKAE`NG1Kb`dGyKqYJk$EX@!-aa-K8iaNK_R{(POsP4UG0 zTYW|SJm{=Fj+wr14v#+hefvPAOZ6=)|BU0EI2ZBC?_;m`xT8Js+Vx{cyxnJ>2j}i( z_sCqVUwM6R@VLqL50r5H_wroPW}jo;Hh&(EKKbVxH{$$uy9mWC7*NacUNoP_GWql1 zO&)iO{RME`=GTrmZ1P;({9CSTi{*|vUUCPIg3k$W#)&K^b*Y~t4t~vbJqPY;xjRVd z(*At*%0ryrcFTF6TrYPAd)zTLa2&V!wIdFjJg?_Ai zZ@Zl1{zUE$3ByhP^h}0bak~!JR30 z$0}W_Z?oLB`g zzuXM3{UC|&C3_y^$H zJolHobCoW|n^**QEAo7c+?^MOTR9W%HuSS!?#@@bW?DX6yMFWK?rln!`gzs({f?RL z?pL&rB+Qe$3zRPP)$VkxBYXY5NABJph8sU>KR++H&t>K`<#P89rAz%BIT-uZx9F=i za<@q7QhlRR&{xNBU05%7@ASBlr{Dkvd!J!P`@iHCC8Mv7qp#kTy9<>rwQo#M)Yo;X z&uVtc-D0In&UuN*v-`Oy#b5S?++Cz}Dc%+FSVup@JtBAS3d2o^ML+)nw@I1 z%sTpB?v^WEGwoUA`2g-G8!vMBbiwa89#*;(UvUNc`E%UooGEu7QMxA2<;b)9IU)6T z!ufJ{h0>*EQ~D70XV<0QqpXwMeN^e1X=SLd{amz{+xQhdkm$2xL97d7$rm%9~8 zm&&Vof z+fiS8pR-!-KBaUi-s20fPr83^X5xKG?ym8;*~0q`vp>&AeOuu^XOrB0TItfX=s9q0 zygTIXGhw*SEUcsExPRU!ch@Ris%Pyq)YrX_HTC^Y?mnw@sh?LAz_rf{zPlVgUGV#j zbxN1wTM)oHvd?9$YL!YJRR$| zyWD+2=~BFzRGrEBs$9DQZ~esZPUeJu>P+Krd|gdY)arQF@%alPN0G{OISk5gVFMxNZSOugTf zyRV1gmb!W7>A3#BFL&Pv!%h5^YvSD_cdL~y?axKSJn!S8c6DqV`Vy5i*bIF$FYc5?S^rAtnd*X<~GHz{3;H{n+< K<=vln()~YhF0gd~ literal 24384 zcmb`NeRver6~-5VD2fu)qI}h~v;-t#bOHHNrs5YE6bhmd5MoI9%Ey`wfgnf_Oo;>% z1kor0f>2P@XcR;dBT$7Jd{7h#U>mgt0#vbJ1eLaP?s+emc^=)9Gwwe(GrMwY3FDYx#48!8bRP$=#JW)cgt_~ju#oQ=gdV4uR^TJqmwU~Rhl56y> z>gVjJa9;@ycyjJJO0J1-aO@6k^9y8ez>`r z9N+K&)Nyw)_W~u?_$n=%InF@Dxto~VOv&Xs%4MCEgX*2%XeHO^9X*|SX@kL)`?oj@g2;qNBv%i3)yu+$~DuX^O>7B6!j~wql=VWli%&LnYZ9h z=(YQf$aQqFl55t5N9J(duDT1l<$nH+l56ztn9F|Nkcc?T>$im;?x`Z?mEDc82kQ<* zj`Jl-uJLouBj61~etU|!mwLFhcEg49wyeB$Wb)g6A;-Dp2yo?kTPnF`+6#*uGcLOa zy5+uVrR178w<^IpG?KYNc^$R(!;M?a@!dEI@s{s%+9T-uK7unBk zJ=`p7B%HT?vLj9hjIY8^v(DUm;j2y#zaN?VTU^NYXDQd@cj|J)HyL{6>(OOOuF*Si z1@qQ{TP@!|$0)fbzlqN=e5Y|7{s}=!&?Bw!+lfOuMD@-diGVj3~=Rrva^y)eKO&j^_<_@4AilF-+rZ%YwGBM4b0o0$^2^h zdK9nZ8b9Cu2J#z1-;l4rS1GwPE$@?!%-tSB-;n$HY7h5okFV^0&W`wa!02tUiFNk6 z&oSroZCR@8dB+9_KDfu37h!%UN&hERJ_exu3f#xkhidw>hu*S?ZYh`aRu;d<+%=QXujGD4VCu~!KZghHV4XuI^JVhn`TdiR zI47KVpYvKf8FB6a8IkjQy@y-n$*-*+JL2O3@^f0Noy;rDWnP~AKD4h7Zs9g?+e~3@ zK;9>BP;$+*1s|}V-kZYluE%qg8$I03o;d$p9nn4~D*PeGIqiPbliXJ~DY?eaIlGwK zW-5;f@;-U9l52d`_ao+2Om)m<^5pCOEgtSZPn@lf?TC*DNbiEok2${cr!lu4_o4li zT(j=4+0A|$`ar$8{XN_}ZNtR*pYxT`yXq76bL;8Mtw)@1RdQ)s`M>rs_t{Gj~#^8gR`HKDI)T1>P3?C1S?@U^^tZ}Y(&^~$Hr8#06SmdV$_ z+dbTso`U%~U(x#8@%$?0R?Xmemvnaca*@}gfl4m*$=$;~W1SNoM4V&A+(900KVg2& zek1*J<~E!OU&;IApFP~Z!Z@40@xm9(E1JnV>ruahm0XH*{_-!mZm*k#`jtNy9pd3G z_vp2UcEra66z2_@2iZ@jWE4+qc1CuAd*5u&-~Ib4U5$ z&iwcg^G>u_T7Y|W3G04J&P`Ersc-Z*Zov+4kHc3R z<=inI?qkCBmwZ+GTn+atSFgEDg`7K9$)#!WMbK+~wNK6+=Yv}{9^6AqSoa}0cf68o z_U+f!aKAbLU!9P1Cn&izt+Z_o*UuNo?>P@Ue7VT$ej0!9bib;PuezOG!<;JgjZ5U* zbS0Pi#^9p|nY$nPZ7=6$D7iE(e)mE4(;novi=3OOT$t8JP1|YxxM4T7Pxj9NM>79D*L5}ZdsH3%V?j$9b zrd4%Bem_T?Uzc+yd$_xBYzO;%!H)Jk^^L?9XY`tUZ_ zcOE;y`Q3*&SIfEgE4dWsf&&K})(Jl!mvg5oxs=bfA0Oa;^%43;!#Vst1kaJCDY>Ng z_-)0JGB7X2LbR=`)C<=i|am-1V<99&ngso&mm z?hGZ@#JL3AorrUSoco}XOY8okB9E^U<=mM{uBoG0;JS5W;yhN)ou%ZG-WIvwy8Aj4 z_epZ@LrN~yZ*Cg;q}^XVB+Z)4I4VQA?7$D~^RB~xr^r>pcOm~0BG2?JK_fa$9%yWKe zpQLG}$E(@LF4y$QadPe=C71kL_Ek0W8bj}7IrlLom!{>_RP#C*wUozvIrniTm-

HN7ubjKY2RC{X@>`2_|Dc?^ zRLQ0Ol=vFB?(cLa&d254WlAo^wS6h#Tmx>Sxei}0@_yrKC71N|UIVWCdy2{Lg>vq4 zC71l1yb^JCpA(um$H=)Wd~n-6js9$(Pj;4bpHXs2-~JNxNvpTFoV!xVrD?^FqHnnO zp(elm<=j~ki{xoebM zs-t$rsNe6}P`H1rc(7R2}eNoA!xE9U=w<-2lyXD+c zC71ks(+v3eT->L8Dd)bV*NBz3b>&)@Za_(zNE=>!K2iIMXOn)9C=WbAPN#BaGsNe6QH(Abo zUCA}|n}Rwz0lgtP_YEbN^d={R>)w~qad5hvyHUv{z3WDyj@vhhu$hVcZ-rs({{SP;^|m_zm{{i z^2F2sKk3cmFXi{hh%;Q{=O5+VzbUyC=N0aA2CJ9)Na_VeIF7u}c}K}L>;8%}T>2it g)^hHwxil^G6RznWUC(g;59Qq0VgLXD diff --git a/autotest/test_prt_dry.py b/autotest/test_prt_dry.py new file mode 100644 index 00000000000..abc11fa7fe5 --- /dev/null +++ b/autotest/test_prt_dry.py @@ -0,0 +1,474 @@ +""" +Tests particle tracking in "dry" conditions. + +The PRP package provides the `DRY` option +to specify how particles should behave in +dry conditions when the flow model enables +the Newton formulation. + +This test case is adapted from the example +simulation provided by @javgs-bd in +https://github.com/MODFLOW-USGS/modflow6/issues/2014. +""" + +import os +from warnings import warn + +import flopy +import matplotlib.colors as clt +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pytest +from framework import TestFramework +from prt_test_utils import get_model_name + +simname = "prtdry" +cases = [ + # expect termination with status 7 immediately in 1st time step + simname, + # expect all particles to be draped to water table and tracked + simname + "_drape", + # the rest of the test cases activate newton and test the behavior + # of the DRY option. with drop, expect all particles to be dropped + # to the highest active cell below and then to be tracked as usual. + simname + "_drop", + # expect termination with status 7 immediately in 1st time step + simname + "_stop", + # expect particles to remain in release positions until the water + # table rises to meet them + simname + "_stay", +] + +nper = 3 +perlen = [1, 300, 1000] +nstp = [1, 30, 1] +tsmult = [1, 1.1, 1] + +Lx = 100.0 +Ly = 100.0 + +nlay = 3 +ncol = 20 +nrow = 20 +arot = 0 + +delr = Lx / ncol +delc = Ly / nrow + +strt = 1600 + +z_top = 1600 +thickness = 10 +xorigin, yorigin = (0, 0) + +z_bot = np.zeros((nlay, nrow, ncol)) +for i in range(nlay): + z_bot[i, :, :] = z_top - (i + 1) * thickness + +period_data = [] +for i in range(nper): + period_data.append((perlen[i], nstp[i], tsmult[i])) + + +def build_gwf_sim(name, gwf_ws, mf6, newton=False): + sim = flopy.mf6.MFSimulation( + sim_name=simname, exe_name=mf6, version="mf6", continue_=True, sim_ws=gwf_ws + ) + + gwf_name = get_model_name(name, "gwf") + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwf_name, + model_nam_file=gwf_name + ".nam", + newtonoptions="NEWTON" if newton else None, + save_flows=True, + ) + + tdis = flopy.mf6.ModflowTdis( + sim, time_units="days", nper=nper, perioddata=period_data + ) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + length_units="METERS", + xorigin=xorigin, + yorigin=yorigin, + angrot=arot, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=z_top, + botm=z_bot, + filename=gwf.name + ".dis", + pname="dis", + ) + + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + csv_outer_output_filerecord="ims_outer.csv", + csv_inner_output_filerecord="ims_inner.csv", + outer_dvclose=1e-5, + outer_maximum=100, + under_relaxation="DBD", + under_relaxation_gamma=0.01, + under_relaxation_theta=0.7, + under_relaxation_kappa=0.01, + under_relaxation_momentum=0.0, + inner_maximum=100, + inner_dvclose=1e-6, + rcloserecord=0.1, + linear_acceleration="BICGSTAB", + relaxation_factor=0.99, + number_orthogonalizations=2, + reordering_method="NONE", + pname=gwf.name, + ) + + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=gwf.name + ".cbb", + head_filerecord=gwf.name + ".hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + ic = flopy.mf6.modflow.mfgwfic.ModflowGwfic(gwf, strt=strt) + + bdry_name = "upstr.chd" + + filename = gwf.name + "." + bdry_name + pname = bdry_name + + stress_period_data = {} + lay = 1 + row, col = (0, 0) + h_upstr = 1587 + + for i in range(nper): + stress_period_data[i] = [lay, row, col, h_upstr, bdry_name] + + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=stress_period_data, + filename=filename, + pname=pname, + boundnames=True, + save_flows=True, + ) + + bdry_name = "dwstr.chd" + + filename = gwf.name + "." + bdry_name + pname = bdry_name + + stress_period_data = {} + lay = 1 + row, col = (nrow - 1, ncol - 1) + h_dwstr = 1582 + + for i in range(nper): + stress_period_data[i] = [lay, row, col, h_dwstr, bdry_name] + + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=stress_period_data, + filename=filename, + pname=pname, + boundnames=True, + save_flows=True, + ) + + bdry_name = "inj.wel" + + filename = gwf.name + "." + bdry_name + pname = bdry_name + + stress_period_data = {} + + lay = 1 + row, col = (int(nrow / 4), int(ncol / 4)) + q = 100 + + for i in range(1, nper): + stress_period_data[i] = [lay, row, col, q, bdry_name] + + wel = flopy.mf6.ModflowGwfwel( + gwf, + stress_period_data=stress_period_data, + pname=pname, + boundnames=True, + save_flows=True, + ) + + icelltype = 1 + + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=icelltype, + k=0.5, + k33=0.1, + save_specific_discharge=True, + save_saturation=True, + ) + + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=0.0001, + sy=0.1, + steady_state={0: True, 2: True}, + transient={1: True}, + ) + return sim + + +def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False): + sim = flopy.mf6.MFSimulation( + sim_name=name, exe_name=mf6, version="mf6", continue_=True, sim_ws=prt_ws + ) + + gwf_ws = gwf.model_ws + + prt_name = get_model_name(name, "prt") + prt = flopy.mf6.ModflowPrt( + sim, + modelname=prt_name, + model_nam_file=prt_name + ".nam", + ) + + tdis = flopy.mf6.ModflowTdis( + sim, time_units="days", nper=nper, perioddata=period_data + ) + + dis = flopy.mf6.ModflowPrtdis( + prt, + length_units="METERS", + xorigin=xorigin, + yorigin=yorigin, + angrot=arot, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=z_top, + botm=z_bot, + filename=prt_name + ".dis", + pname="dis", + ) + + porosity = 0.1 + + flopy.mf6.ModflowPrtmip(prt, porosity=porosity, pname="mip") + + # particles are released in layer 0 + offsets = [ + (-1, 0, 0), + (-1, -1, 0), + (-1, 1, 0), + (-1, 0, -1), + (-1, 0, 1), + ] + + lay = 1 + row, col = (int(nrow / 4), int(ncol / 4)) + prp_cells = [(lay + k, row + i, col + j) for k, i, j in offsets] + + prp_coords = [] + prp_data = [] + + for i, (lay, row, col) in enumerate(prp_cells): + x, y, z = ( + gwf.modelgrid.xyzcellcenters[0][row, col], + gwf.modelgrid.xyzcellcenters[1][row, col], + gwf.modelgrid.xyzcellcenters[2][lay, row, col], + ) + + prp_lst = [i, (lay, row, col), x, y, z] + + prp_data.append(prp_lst) + prp_coords.append((x, y, z)) + + flopy.mf6.ModflowPrtprp( + prt, + nreleasepts=len(prp_data), + packagedata=prp_data, + nreleasetimes=1, + releasetimes=[(0.0,)], + exit_solve_tolerance=1e-7, + drape=drape, + dry_tracking_method=dry_tracking_method, + pname="prp", + filename="tracking_1.prp", + print_input=True, + ) + + budgetfile = prt.name + ".bud" + trackfile = prt.name + ".trk" + trackcsvfile = prt.name + ".csv" + + flopy.mf6.ModflowPrtoc( + prt, + budget_filerecord=[budgetfile], + track_filerecord=[trackfile], + trackcsv_filerecord=[trackcsvfile], + saverecord=[("BUDGET", "ALL")], + pname="oc", + ) + + rel_prt_folder = os.path.relpath(gwf_ws, start=prt_ws) + + packagedata = [ + ("GWFHEAD", f"{rel_prt_folder}/{gwf.name}.hds"), + ("GWFBUDGET", f"{rel_prt_folder}/{gwf.name}.cbb"), + ] + + flopy.mf6.ModflowPrtfmi( + prt, + packagedata=packagedata, + ) + + ems = flopy.mf6.ModflowEms(sim, pname="ems", filename=prt.name + ".ems") + + return sim + + +def build_models(idx, test, newton, drape=False, dry_tracking_method=False): + gwf_sim = build_gwf_sim( + test.name, test.workspace / "gwf", test.targets["mf6"], newton=newton + ) + prt_sim = build_prt_sim( + test.name, + gwf_sim.get_model(), + test.workspace / "prt", + test.targets["mf6"], + drape=drape, + dry_tracking_method=dry_tracking_method, + ) + return gwf_sim, prt_sim + + +def check_output(idx, test, snapshot): + name = test.name + gwf_ws = test.workspace / "gwf" + prt_ws = test.workspace / "prt" + gwf_name = get_model_name(name, "gwf") + prt_name = get_model_name(name, "prt") + gwf_sim = test.sims[0] + gwf = gwf_sim.get_model(gwf_name) + hds_file = gwf.name + ".hds" + trackcsv_file = prt_name + ".csv" + trackcsv_path = prt_ws / trackcsv_file + pls = pd.read_csv(trackcsv_path) + strtpts = pls[pls.ireason == 0] + + # compare to expected results + decimals = 1 if "drop" in name else 2 + actual = pls.drop(["name", "icell"], axis=1).round(decimals).reset_index(drop=True) + # ignore particle 4, it terminates early with optimization=2 when built with ifort + if "drop" in name: + actual = actual.drop(actual[actual.irpt == 4].index) + assert snapshot == actual.to_records(index=False) + + plot_pathlines = False + if plot_pathlines: + + def plot_pathlines_and_timeseries( + ax, mg, ibd, pathlines, timeseries, plottitle, layer=1 + ): + ax.set_aspect("equal") + mm = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=layer) + mm.plot_grid(color=(0.4, 0.4, 0.4, 0.5), lw=0.2) + mm.plot_grid(color=(0.4, 0.4, 0.4, 0.5), lw=0.2) + mm.plot_bc("WEL", plotAll=True) + mm.plot_bc("CHD", plotAll=True) + mm.plot_grid(lw=0.5) + # mm.plot_array(gwf.output.head().get_data()) + v = mm.plot_array(ibd, cmap=cmapbd, edgecolor="gray") + plt.scatter(pathlines["x"], pathlines["y"]) + mm.plot_pathline(pathlines, layer="all", colors=["blue"], lw=0.75) + ax.set_title(plottitle, fontsize=12) + ax.scatter(strtpts.x, strtpts.y) + from shapely.geometry import Polygon + + cellids = [105, 85, 125, 104, 106] + polys = [Polygon(gwf.modelgrid.get_cell_vertices(ic)) for ic in cellids] + mm.plot_shapes(polys, alpha=0.2) + plt.show() + + fig, ax = plt.subplots(1, 1, figsize=(8, 8)) + cmapbd = clt.ListedColormap(["b", "r"]) + chdcells = [] + for pckge in gwf.get_package("chd"): + chd_nodes = pckge.stress_period_data.get_data(0).cellid + for item in chd_nodes: + chdcells.append(item) + + welcells = [] + pckg = gwf.get_package("wel") + wel_nodes = pckg.stress_period_data.get_data(1).cellid + for item in wel_nodes: + welcells.append(item) + + # identify the boundary locations + ibd = np.zeros((nlay, nrow, ncol), dtype=int) + ilay, irow, icol = zip(*chdcells) + ibd[ilay, irow, icol] = 1 + ilay, irow, icol = zip(*welcells) + ibd[ilay, irow, icol] = 2 + ibd = np.ma.masked_equal(ibd, 0) + + plot_pathlines_and_timeseries(ax, gwf.modelgrid, ibd, pls, None, name) + + plot_3d = False + if plot_3d: + try: + import pyvista as pv + from flopy.export.vtk import Vtk + except: + warn("Couldn't make 3D plots, need pyvista/vtk") + return + + vert_exag = 1 + vtk = Vtk(model=gwf, binary=False, vertical_exageration=vert_exag, smooth=False) + vtk.add_model(gwf) + vtk.add_pathline_points(pls.to_records(index=False)) + gwf_mesh, prt_mesh = vtk.to_pyvista() + axes = pv.Axes(show_actor=False, actor_scale=2.0, line_width=5) + p = pv.Plotter(window_size=[700, 700]) + p.enable_anti_aliasing() + p.add_mesh(gwf_mesh, opacity=0.025, style="wireframe") + p.add_mesh( + prt_mesh, + point_size=8, + line_width=2.5, + smooth_shading=True, + color="blue", + ) + p.show() + + +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets, array_snapshot): + dry_tracking_methods = ["drop", "stop", "stay"] + if any(t in name for t in dry_tracking_methods): + dry_tracking_method = name[-4:] + else: + dry_tracking_method = None + newton = any(t in name for t in dry_tracking_methods) + drape = "drape" in name + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t, newton, drape, dry_tracking_method), + check=lambda t: check_output(idx, t, array_snapshot), + targets=targets, + compare=None, + ) + test.run() diff --git a/autotest/test_prt_disv1.py b/autotest/test_prt_newton_disv1.py similarity index 98% rename from autotest/test_prt_disv1.py rename to autotest/test_prt_newton_disv1.py index d9d62e7793a..bbba8e6ad04 100644 --- a/autotest/test_prt_disv1.py +++ b/autotest/test_prt_newton_disv1.py @@ -6,13 +6,7 @@ refined cells, instead of the new ternary method which applies more generally to polygonal cells. -The simulation includes a single stress period -with multiple time steps. This serves to test -whether PRT properly solves trajectories over -"internal" time steps, i.e. within the step's -slice of simulation time, as well as extending -tracking to termination or particle stop times -during the simulation's final time step. +The simulation uses the Newton formulation. Several cases are provided: - default: No user-specified tracking times, MP7 in pathline mode. diff --git a/autotest/test_prt_voronoi1.py b/autotest/test_prt_voronoi1.py index 16e0ebf1d65..a36eee6c4b4 100644 --- a/autotest/test_prt_voronoi1.py +++ b/autotest/test_prt_voronoi1.py @@ -10,7 +10,6 @@ """ from pathlib import Path -from platform import processor, system import flopy import matplotlib as mpl @@ -267,6 +266,7 @@ def plot_output(idx, test): ax = plt.subplot(1, 1, 1, aspect="equal") pmv = flopy.plot.PlotMapView(model=gwf, ax=ax) pmv.plot_grid(alpha=0.25) + pmv.plot_ibound(alpha=0.5) headmesh = pmv.plot_array(head, alpha=0.25) cv = pmv.contour_array(head, levels=np.linspace(0, 1, 9), colors="black") @@ -331,7 +331,7 @@ def plot_output(idx, test): plt.savefig(prt_ws / f"{name}.png") -def check_output(idx, test, snapshot): +def check_output(idx, test): name = test.name prt_ws = test.workspace / "prt" prt_name = get_model_name(name, "prt") @@ -349,22 +349,20 @@ def check_output(idx, test, snapshot): pls = pd.read_csv(prt_ws / prt_track_csv_file, na_filter=False) endpts = pls[pls.ireason == 3] # termination - # compare pathlines with snapshot. particles shouldn't - # have moved vertically. round for cross-platform error. - # skip macos-14 in CI because grid is slightly different - if not (system() == "Darwin" and processor() == "arm"): - assert snapshot == endpts.drop("name", axis=1).round(1).to_records(index=False) + assert np.allclose(endpts.z, 0.5) + assert np.isclose(endpts.y.min(), 1, atol=4) + assert np.isclose(endpts.y.max(), 996, atol=4) @requires_pkg("syrupy") @pytest.mark.slow @pytest.mark.parametrize("idx, name", enumerate(cases)) -def test_mf6model(idx, name, function_tmpdir, targets, benchmark, array_snapshot, plot): +def test_mf6model(idx, name, function_tmpdir, targets, benchmark, plot): test = TestFramework( name=name, workspace=function_tmpdir, build=lambda t: build_models(idx, t), - check=lambda t: check_output(idx, t, array_snapshot), + check=lambda t: check_output(idx, t), plot=lambda t: plot_output(idx, t) if plot else None, targets=targets, compare=None, diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index e222b010532..4c5e63e92cc 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -51,8 +51,9 @@ \item The number of characters used to represent integers and floating point numbers in MODFLOW input files was restricted to 30. The program was modified to accept any number of characters provided the number is valid. This may be useful for parameter estimation programs that use character substitution to create new input files. \item A string to character array conversion function in the BMI interface could fail on Apple silicon macOS with recent versions of GNU Fortran, producing array extent errors at runtime. This has been fixed by properly specifying the intent of a dummy argument in the relevant function. \item The UZF Package will facilitate UZF objects with areas less than the area of the host cell. However, within the GWE model type, the UZE package will not work properly when the area of UZF objects is not equal to the area of the host cell. New code was added to ensure that the area of each UZF object is equal to that of the host grid cell. When this condition is not satisfied, the new code will stop the simulation with an error message indicating which cell is in violation. - % \item xxx - + \item With a flow model using the Newton formulation, the PRT model could crash upon a particle's entry into a dry cell. This has been fixed. + \item With a flow model using the Newton formulation, the PRT model could enter an endless loop upon a particle's entry to a dry cell if that cell contains a boundary package (e.g. a pumping well). Where the particle should be captured and terminate, it would instead be passed back and forth between the cell bottom and the top of the cell below. To avoid this, particles are forbidden from backtracking (reentering the previous cell) within the same time step. + \item The PRT model now allows more control over vertical particle motion in dry conditions. In addition to the existing DRAPE option, which controls release-time behavior, the PRP package now provides a DRY\_TRACKING\_METHOD option which configures how dry particles (particles in dry cells, or above the water table in partially saturated cells) behave at tracking time. This option is relevant only when the Newton formulation is used, in which case dry cells remain active; otherwise, dry cells are inactive and particles will terminate. See the MF6IO document for a detailed explanation of DRY\_TRACKING\_METHOD. \end{itemize} %\underline{INTERNAL FLOW PACKAGES} diff --git a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn index bb5edd48604..76aa68560c8 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-prp.dfn +++ b/doc/mf6io/mf6ivar/dfn/prt-prp.dfn @@ -167,6 +167,15 @@ optional true longname drape description is a text keyword to indicate that if a particle's release point is in a cell that happens to be inactive at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is inactive at release time. +block options +name dry_tracking_method +type string +valid drop stop stay +reader urword +optional true +longname what to do in dry-but-active cells +description is a string indicating how particles should behave in dry-but-active cells (as can occur with the Newton formulation). The value can be ``DROP'', ``STOP'', or ``STAY''. The default is ``DROP'', which passes particles vertically and instantaneously to the water table. ``STOP'' causes particles to terminate. ``STAY'' causes particles to remain stationary but active. + block options name dev_forceternary type keyword diff --git a/doc/mf6io/prt/prp.tex b/doc/mf6io/prt/prp.tex index 473e1c930b7..5b3ac83e6b9 100644 --- a/doc/mf6io/prt/prp.tex +++ b/doc/mf6io/prt/prp.tex @@ -2,6 +2,8 @@ The PRP Package offers multiple ways to specify particle release times. Particle release times may either be provided explicitly, relative to the simulation start time, or configured relative to the time discretization of stress periods via period block settings. When multiple ways of specifying release times are used together, the resulting set of release times is the union of the times specified by each method. +This package contains options to configure how particles behave in dry conditions. Each of these is described in detail below. Their interaction is summarized in the PRT chapter introduction. + \vspace{5mm} \subsubsection{Structure of Blocks} \lstinputlisting[style=blockdefinition]{./mf6ivar/tex/prt-prp-options.dat} diff --git a/doc/mf6io/prt/prt.tex b/doc/mf6io/prt/prt.tex index 196eba092c6..385f45f4260 100644 --- a/doc/mf6io/prt/prt.tex +++ b/doc/mf6io/prt/prt.tex @@ -21,6 +21,42 @@ \subsection{Specifying Cell Face Flows using IFLOWFACE} \subsection{Particle Mass Budget} A summary of all inflow (sources) and outflow (sinks) of particle mass is called a mass budget. The particle mass budget is printed to the PRT Model Listing File for selected time steps. In the current implementation, each particle is assigned unit mass, and the numerical value of the flow can be interpreted as particles per time. +\subsection{Vertical Tracking} + +When a particle is in the flow field, vertical motion can be solved in the same way as lateral motion. Special handling is necessary above the water table. + +A dry cell is either an inactive cell or an active but dry cell, as can occur with the Newton formulation. + +Normally, an inactive cell might be dry or explicitly disabled (idomain). With Newton, dry cells remain active. + +Release-time and tracking-time considerations are described (and implemented) separately. + +\subsubsection{Release} + +At release time, PRT decides whether to release each particle or to terminate it unreleased. + +If the release cell is active, the particle will be released at the user-specified location. + +If the release cell is inactive, behavior is determined by the DRAPE option. If the DRAPE option is enabled, the particle will be released from the top-most active cell beneath it, if any. If there is no active cell underneath the particle in any layer, or if DRAPE is not enabled, the particle will terminate unreleased (with status code 8). + +Since under the Newton formulation dry cells can remain active, the DRAPE option has no effect when Newton is on (assuming particles are not released into disabled grid regions). Vertical tracking behavior with Newton can be configured using the DRY\_TRACKING\_METHOD option discussed below. + +\subsubsection{Tracking} + +With the Newton formulation, particles can be released into dry-but-active cells. + +Particle trajectories are solved over the same time discretization used by the flow model. A particle may be immersed in the flow field in one time step, and find that the water table has dropped below it in the next. + +A particle which finds itself in an inactive cell will terminate with status code 7. This is consistent with MODPATH 7's behavior. + +A particle in a dry but active cell, or above the water table in a partially saturated cell, need not terminate. We call such a particle dry. The PRP package provides an option DRY\_TRACKING\_METHOD determining how dry particles should behave. Supported values are DROP (default), STOP, and STAY. + +If DROP is selected, or if a DRY\_TRACKING\_METHOD is unspecified, a dry particle is passed vertically and instantaneously to the water table (if the cell is partially saturated) or to the bottom of the cell (if the cell is dry). This repeats (i.e. the particle may drop through multiple cells) until it reaches the water table. Tracking then proceeds as usual. If the vertical column containing the particle is entirely dry, the particle will terminate upon reaching the bottom of the model grid. + +If STOP is selected, dry particles will be terminated. + +If STAY is selected, a dry particle will remain stationary until a) the water table rises and tracking can continue, b) the particle terminates due to reaching its STOPTIME or STOPTRAVELTIME, or c) the simulation ends. + \subsection{Particle Track Output} The PRT Model supports both binary and CSV particle track output files. A particle track CSV file contains the output data in tabular format. The first line of the CSV file contains column names. Each subsequent line in the file contains a row of data for a single particle track record, with the following fields: diff --git a/src/Model/ModelUtilities/TrackFile.f90 b/src/Model/ModelUtilities/TrackFile.f90 index 4ba94b71e9c..a09be74b7a8 100644 --- a/src/Model/ModelUtilities/TrackFile.f90 +++ b/src/Model/ModelUtilities/TrackFile.f90 @@ -35,7 +35,7 @@ module TrackFileModule !! 1: particle transitioned between cells !! 2: current time step ended**** !! 3: particle terminated - !! 4: particle exited weak sink + !! 4: particle in weak sink !! 5: user-specified tracking time !! !! Each record has an "istatus" property, which is the tracking status; diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90 index fe69f888cec..96ebf145415 100644 --- a/src/Model/ParticleTracking/prt-prp.f90 +++ b/src/Model/ParticleTracking/prt-prp.f90 @@ -49,6 +49,7 @@ module PrtPrpModule integer(I4B), pointer :: istopweaksink => null() !< weak sink option: 0 = no stop, 1 = stop integer(I4B), pointer :: istopzone => null() !< optional stop zone number: 0 = no stop zone integer(I4B), pointer :: idrape => null() !< drape option: 0 = do not drape, 1 = drape to topmost active cell + integer(I4B), pointer :: idrymeth => null() !< dry tracking method: 0 = drop, 1 = stop, 2 = stay integer(I4B), pointer :: itrkout => null() !< binary track file integer(I4B), pointer :: itrkhdr => null() !< track header file integer(I4B), pointer :: itrkcsv => null() !< CSV track file @@ -89,6 +90,7 @@ module PrtPrpModule procedure :: release procedure :: log_release procedure :: validate_release_point + procedure :: initialize_particle procedure, public :: bnd_obs_supported => prp_obs_supported procedure, public :: bnd_df_obs => prp_df_obs end type PrtPrpType @@ -158,6 +160,7 @@ subroutine prp_da(this) call mem_deallocate(this%istopweaksink) call mem_deallocate(this%istopzone) call mem_deallocate(this%idrape) + call mem_deallocate(this%idrymeth) call mem_deallocate(this%nreleasepoints) call mem_deallocate(this%nreleasetimes) call mem_deallocate(this%nparticles) @@ -247,6 +250,7 @@ subroutine prp_allocate_scalars(this) call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath) call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath) call mem_allocate(this%idrape, 'IDRAPE', this%memoryPath) + call mem_allocate(this%idrymeth, 'IDRYMETH', this%memoryPath) call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath) call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath) call mem_allocate(this%nparticles, 'NPARTICLES', this%memoryPath) @@ -270,6 +274,7 @@ subroutine prp_allocate_scalars(this) this%istopweaksink = 0 this%istopzone = 0 this%idrape = 0 + this%idrymeth = 0 this%nreleasepoints = 0 this%nreleasetimes = 0 this%nparticles = 0 @@ -413,47 +418,52 @@ subroutine release(this, ip, trelease) integer(I4B), intent(in) :: ip !< particle index real(DP), intent(in) :: trelease !< release time ! local - integer(I4B) :: irow, icol, ilay, icpl - integer(I4B) :: ic, icu integer(I4B) :: np - real(DP) :: x, y, z - real(DP) :: top, bot, hds type(ParticleType), pointer :: particle + call this%initialize_particle(particle, ip, trelease) + ! Increment cumulative particle count np = this%nparticles + 1 this%nparticles = np - ! Get reduced and user node numbers - ic = this%rptnode(ip) - icu = this%dis%get_nodeuser(ic) + ! Save the particle to the store + call this%particles%save_particle(particle, np) + deallocate (particle) - ! Load coordinates and transform if needed - x = this%rptx(ip) - y = this%rpty(ip) - if (this%ilocalz > 0) then - top = this%fmi%dis%top(ic) - bot = this%fmi%dis%bot(ic) - hds = this%fmi%gwfhead(ic) - z = bot + this%rptz(ip) * (hds - bot) - else - z = this%rptz(ip) - end if + ! Accumulate mass for release point + this%rptm(ip) = this%rptm(ip) + DONE - ! Make sure release point is in the grid/cell - call this%validate_release_point(ic, x, y, z) + end subroutine release + + subroutine initialize_particle(this, particle, ip, trelease) + class(PrtPrpType), intent(inout) :: this !< this instance + type(ParticleType), pointer, intent(inout) :: particle !< the particle + integer(I4B), intent(in) :: ip !< particle index + real(DP), intent(in) :: trelease !< release time + ! local + integer(I4B) :: irow, icol, ilay, icpl + integer(I4B) :: ic, icu, ic_old + real(DP) :: x, y, z + real(DP) :: top, bot, hds + + ic = this%rptnode(ip) + icu = this%dis%get_nodeuser(ic) - ! Initialize the particle call create_particle(particle) + if (size(this%boundname) /= 0) then particle%name = this%boundname(ip) else particle%name = '' end if + particle%irpt = ip particle%istopweaksink = this%istopweaksink particle%istopzone = this%istopzone + particle%idrymeth = this%idrymeth particle%icu = icu + select type (dis => this%dis) type is (DisType) call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay) @@ -462,28 +472,49 @@ subroutine release(this, ip, trelease) end select particle%ilay = ilay particle%izone = this%rptzone(ic) + particle%istatus = 0 - ! Handle inactive cells + ! If the cell is inactive, either drape the particle + ! to the top-most active cell beneath it if drape is + ! enabled, or else terminate permanently unreleased. if (this%ibound(ic) == 0) then - ! If drape option activated, release in highest active - ! cell vertically below release point. - if (this%idrape /= 0) & + ic_old = ic + if (this%idrape > 0) then call this%dis%highest_active(ic, this%ibound) - ! If returned cell is inactive, do not release particle - if (this%ibound(ic) == 0) & - particle%istatus = 8 ! permanently unreleased + if (ic == ic_old .or. this%ibound(ic) == 0) & + particle%istatus = 8 + else + particle%istatus = 8 + end if + end if + + ! Load coordinates and transform if needed + x = this%rptx(ip) + y = this%rpty(ip) + if (this%ilocalz > 0) then + top = this%fmi%dis%top(ic) + bot = this%fmi%dis%bot(ic) + hds = this%fmi%gwfhead(ic) + z = bot + this%rptz(ip) * (hds - bot) + else + z = this%rptz(ip) end if + + call this%validate_release_point(ic, x, y, z) + particle%x = x particle%y = y particle%z = z particle%trelease = trelease - ! Set stopping time to earlier of times specified by STOPTIME and STOPTRAVELTIME + + ! Set stop time to earlier of STOPTIME and STOPTRAVELTIME if (this%stoptraveltime == huge(1d0)) then particle%tstop = this%stoptime else particle%tstop = particle%trelease + this%stoptraveltime if (this%stoptime < particle%tstop) particle%tstop = this%stoptime end if + particle%ttrack = particle%trelease particle%idomain(1) = 0 particle%iboundary(1) = 0 @@ -495,15 +526,7 @@ subroutine release(this, ip, trelease) particle%iexmeth = this%iexmeth particle%iextend = this%iextend particle%extol = this%extol - - ! Store the particle's state and deallocate it - call this%particles%save_particle(particle, np) - deallocate (particle) - - ! Accumulate particle mass for release point - this%rptm(ip) = this%rptm(ip) + DONE - - end subroutine release + end subroutine initialize_particle !> @ brief Read and prepare period data for particle input subroutine prp_rp(this) @@ -691,6 +714,25 @@ subroutine prp_options(this, option, found) case ('DRAPE') this%idrape = 1 found = .true. + case ('DRY_TRACKING_METHOD') + call this%parser%GetStringCaps(keyword) + select case (keyword) + case ('DROP') + this%idrymeth = 0 + case ('STOP') + this%idrymeth = 1 + case ('STAY') + this%idrymeth = 2 + case default + write (errmsg, '(a, a)') & + 'Unknown dry tracking method: ', trim(keyword) + call store_error(errmsg) + write (errmsg, '(a, a)') & + 'DRY must be "DROP", "STOP" or "STAY"' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end select + found = .true. case ('TRACK') call this%parser%GetStringCaps(keyword) if (keyword == 'FILEOUT') then @@ -891,7 +933,7 @@ subroutine prp_read_packagedata(this) if (n < 1 .or. n > this%nreleasepoints) then write (errmsg, '(a,1x,i0,a)') & - 'Release point number must be greater than 0 and less than ', & + 'Release point number must be greater than 0 and less than', & 'or equal to', this%nreleasepoints, '.' call store_error(errmsg) cycle diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90 index e47f8ff6a0b..0730d8f6b02 100644 --- a/src/Model/ParticleTracking/prt.f90 +++ b/src/Model/ParticleTracking/prt.f90 @@ -932,7 +932,7 @@ subroutine prt_solve(this) ! -- If particle is permanently unreleased, record its initial/terminal state if (particle%istatus == 8) & - call this%method%save(particle, reason=3) ! reason=3: termination + call this%method%save(particle, reason=3) ! If particle is inactive or not yet to be released, cycle if (particle%istatus > 1) cycle @@ -940,7 +940,7 @@ subroutine prt_solve(this) ! If particle released this time step, record its initial state particle%istatus = 1 if (particle%trelease >= totimc) & - call this%method%save(particle, reason=0) ! reason=0: release + call this%method%save(particle, reason=0) ! Maximum time is the end of the time step or the particle ! stop time, whichever comes first, unless it's the final @@ -957,6 +957,10 @@ subroutine prt_solve(this) ! Get and apply the tracking method call this%method%apply(particle, tmax) + ! Reset previous cell, exit face, and zone numbers + particle%icp = 0 + particle%izp = 0 + ! Update particle storage call packobj%particles%save_particle(particle, np) end do diff --git a/src/Solution/ParticleTracker/CellDefn.f90 b/src/Solution/ParticleTracker/CellDefn.f90 index 6075a38a399..836f0421633 100644 --- a/src/Solution/ParticleTracker/CellDefn.f90 +++ b/src/Solution/ParticleTracker/CellDefn.f90 @@ -1,5 +1,7 @@ module CellDefnModule use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: DZERO + use MathUtilModule, only: is_close implicit none private @@ -15,6 +17,7 @@ module CellDefnModule integer(I4B), public :: npolyverts !< number of vertices for cell polygon real(DP), public :: porosity !< cell porosity real(DP), public :: retfactor !< cell retardation factor + integer(I4B), public :: ilay !< layer number integer(I4B), public :: izone !< cell zone number integer(I4B), public :: iweaksink !< weak sink indicator integer(I4B), public :: inoexitface !< no exit face indicator diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90 index 38049b5d51a..061a627e1ad 100644 --- a/src/Solution/ParticleTracker/Method.f90 +++ b/src/Solution/ParticleTracker/Method.f90 @@ -2,6 +2,7 @@ module MethodModule use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: DZERO use ErrorUtilModule, only: pstop use SubcellModule, only: SubcellType use ParticleModule @@ -11,6 +12,7 @@ module MethodModule use CellDefnModule, only: CellDefnType use TrackControlModule, only: TrackControlType use TimeSelectModule, only: TimeSelectType + use MathUtilModule, only: is_close implicit none private @@ -27,7 +29,7 @@ module MethodModule !! depending on cell geometry (implementing the strategy pattern). !< type, abstract :: MethodType - character(len=40), pointer, public :: type !< method name + character(len=40), pointer, public :: name !< method name logical(LGP), public :: delegates !< whether the method delegates type(PrtFmiType), pointer, public :: fmi => null() !< ptr to fmi class(CellType), pointer, public :: cell => null() !< ptr to the current cell @@ -50,7 +52,7 @@ module MethodModule procedure :: save procedure :: track procedure :: try_pass - procedure :: update + procedure :: check end type MethodType abstract interface @@ -94,7 +96,8 @@ subroutine init(this, fmi, cell, subcell, trackctl, tracktimes, & if (present(retfactor)) this%retfactor => retfactor end subroutine init - !> @brief Track particle through subdomains + !> @brief Track the particle over domains of the given + ! level until the particle terminates or tmax elapses. recursive subroutine track(this, particle, level, tmax) ! dummy class(MethodType), intent(inout) :: this @@ -106,6 +109,7 @@ recursive subroutine track(this, particle, level, tmax) integer(I4B) :: nextlevel class(methodType), pointer :: submethod + ! Advance the particle over subdomains advancing = .true. nextlevel = level + 1 do while (advancing) @@ -115,28 +119,28 @@ recursive subroutine track(this, particle, level, tmax) end do end subroutine track - !> @brief Try passing the particle to the next subdomain + !> @brief Try passing the particle to the next subdomain. subroutine try_pass(this, particle, nextlevel, advancing) class(MethodType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle integer(I4B) :: nextlevel logical(LGP) :: advancing - ! tracking submethod marked tracking complete? - ! reset domain boundary flag and don't advance + ! if the particle is done advancing, reset the domain boundary flag. if (.not. particle%advancing) then particle%iboundary = 0 advancing = .false. else - ! otherwise pass particle to next subdomain - ! and if it's on a boundary, stop advancing + ! otherwise pass the particle to the next subdomain. + ! if that leaves it on a boundary, stop advancing. call this%pass(particle) - if (particle%iboundary(nextlevel - 1) .ne. 0) & + if (particle%iboundary(nextlevel - 1) .ne. 0) then advancing = .false. + end if end if end subroutine try_pass - !> @brief Load subdomain tracking method (submethod) + !> @brief Load the subdomain tracking method (submethod). subroutine load(this, particle, next_level, submethod) class(MethodType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle @@ -145,14 +149,14 @@ subroutine load(this, particle, next_level, submethod) call pstop(1, "load must be overridden") end subroutine load - !> @brief Pass a particle to the next subdomain, internal use only + !> @brief Pass the particle to the next subdomain. subroutine pass(this, particle) class(MethodType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle call pstop(1, "pass must be overridden") end subroutine pass - !> @brief Save a particle's current state. + !> @brief Save the particle's state to output files. subroutine save(this, particle, reason) use TdisModule, only: kper, kstp, totimc ! dummy @@ -179,50 +183,106 @@ subroutine save(this, particle, reason) end if end if - ! Save the particle's state to any registered tracking output files - call this%trackctl%save(particle, kper=per, & - kstp=stp, reason=reason) + call this%trackctl%save(particle, kper=per, kstp=stp, reason=reason) end subroutine save - !> @brief Update particle state and check termination conditions + !> @brief Check reporting/terminating conditions before tracking. !! - !! Update the particle's properties (e.g. advancing flag, zone number, - !! status). If any termination conditions apply, the particle's status - !! will be set to the appropriate termination value. If any reporting - !! conditions apply, save particle state with the proper reason code. + !! Check a number of conditions determining whether to continue + !! tracking the particle or terminate it, as well as whether to + !! record any output data as per selected reporting conditions. !< - subroutine update(this, particle, cell_defn) + subroutine check(this, particle, cell_defn) + ! modules + use TdisModule, only: endofsimulation, totim ! dummy class(MethodType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle type(CellDefnType), pointer, intent(inout) :: cell_defn + ! local + logical(LGP) :: dry_cell, dry_particle, no_exit_face, stop_zone, weak_sink + + dry_cell = this%fmi%ibdgwfsat0(cell_defn%icell) == 0 + dry_particle = particle%z > cell_defn%top + no_exit_face = cell_defn%inoexitface > 0 + stop_zone = cell_defn%izone > 0 .and. particle%istopzone == cell_defn%izone + weak_sink = cell_defn%iweaksink > 0 particle%izone = cell_defn%izone - if (cell_defn%izone .ne. 0) then - if (particle%istopzone .eq. cell_defn%izone) then - particle%advancing = .false. - particle%istatus = 6 - call this%save(particle, reason=3) ! particle terminated - return - end if + if (stop_zone) then + particle%advancing = .false. + particle%istatus = 6 + call this%save(particle, reason=3) + return end if - if (cell_defn%inoexitface .ne. 0) then + + if (no_exit_face .and. .not. dry_cell) then particle%advancing = .false. particle%istatus = 5 - call this%save(particle, reason=3) ! particle terminated + call this%save(particle, reason=3) return end if - if (cell_defn%iweaksink .ne. 0) then - if (particle%istopweaksink .ne. 0) then + + if (weak_sink) then + if (particle%istopweaksink > 0) then particle%advancing = .false. particle%istatus = 3 - call this%save(particle, reason=3) ! particle terminated + call this%save(particle, reason=3) return else - call this%save(particle, reason=4) ! particle exited weak sink + call this%save(particle, reason=4) + end if + end if + + if (dry_cell) then + if (particle%idrymeth == 0) then + no_exit_face = .false. + else if (particle%idrymeth == 1) then + ! stop + particle%advancing = .false. + particle%istatus = 7 + call this%save(particle, reason=3) + return + else if (particle%idrymeth == 2) then + ! stay + no_exit_face = .false. + particle%advancing = .false. + particle%ttrack = totim + ! terminate if last period/step + if (endofsimulation) then + particle%istatus = 5 + call this%save(particle, reason=3) + return + end if + call this%save(particle, reason=2) + end if + else if (dry_particle .and. this%name /= "passtobottom") then + ! dry particle + if (particle%idrymeth == 0) then + ! drop to water table + particle%z = cell_defn%top + call this%save(particle, reason=1) + else if (particle%idrymeth == 1) then + ! terminate + particle%advancing = .false. + particle%istatus = 7 + call this%save(particle, reason=3) + return + else if (particle%idrymeth == 2) then + ! stay + no_exit_face = .false. + particle%advancing = .false. return end if end if - end subroutine update + + if (no_exit_face) then + particle%advancing = .false. + particle%istatus = 5 + call this%save(particle, reason=3) + return + end if + + end subroutine check end module MethodModule diff --git a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 b/src/Solution/ParticleTracker/MethodCellPassToBot.f90 index 33825b5b7b6..433b3ca9650 100644 --- a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 +++ b/src/Solution/ParticleTracker/MethodCellPassToBot.f90 @@ -28,15 +28,15 @@ module MethodCellPassToBotModule subroutine create_method_cell_ptb(method) type(MethodCellPassToBotType), pointer :: method allocate (method) - allocate (method%type) - method%type = "passtobottom" + allocate (method%name) + method%name = "passtobottom" method%delegates = .false. end subroutine create_method_cell_ptb !> @brief Deallocate the pass-to-bottom tracking method subroutine deallocate (this) class(MethodCellPassToBotType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine deallocate !> @brief Pass particle vertically and instantaneously to the cell bottom @@ -46,10 +46,15 @@ subroutine apply_ptb(this, particle, tmax) type(ParticleType), pointer, intent(inout) :: particle real(DP), intent(in) :: tmax - call this%update(particle, this%cell%defn) + ! Check termination/reporting conditions + call this%check(particle, this%cell%defn) if (.not. particle%advancing) return + + ! Pass to bottom face particle%z = this%cell%defn%bot particle%iboundary(2) = this%cell%defn%npolyverts + 2 + + ! Save datum call this%save(particle, reason=1) end subroutine apply_ptb diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90 index 71fc327fef7..86b1d09f4ca 100644 --- a/src/Solution/ParticleTracker/MethodCellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollock.f90 @@ -36,7 +36,7 @@ subroutine create_method_cell_pollock(method) allocate (method) call create_cell_rect(cell) method%cell => cell - method%type => method%cell%type + method%name => method%cell%type method%delegates = .true. call create_subcell_rect(subcell) method%subcell => subcell @@ -45,7 +45,7 @@ end subroutine create_method_cell_pollock !> @brief Destroy the tracking method subroutine destroy_mcp(this) class(MethodCellPollockType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine destroy_mcp !> @brief Load subcell tracking method @@ -129,19 +129,12 @@ subroutine apply_mcp(this, particle, tmax) select type (cell => this%cell) type is (CellRectType) - ! Update particle state, return early if done advancing - call this%update(particle, cell%defn) + ! Check termination/reporting conditions + call this%check(particle, cell%defn) if (.not. particle%advancing) return - ! If the particle is above the top of the cell (presumed water table) - ! pass it vertically and instantaneously to the top - if (particle%z > cell%defn%top) then - particle%z = cell%defn%top - call this%save(particle, reason=1) - end if - - ! Transform particle location into local cell coordinates - ! (translated and rotated but not scaled relative to model). + ! Transform model coordinates to local cell coordinates + ! (translated/rotated but not scaled relative to model) xOrigin = cell%xOrigin yOrigin = cell%yOrigin zOrigin = cell%zOrigin @@ -150,14 +143,13 @@ subroutine apply_mcp(this, particle, tmax) call particle%transform(xOrigin, yOrigin, zOrigin, & sinrot, cosrot) - ! Track the particle across the cell. + ! Track the particle over the cell call this%track(particle, 2, tmax) - ! Transform particle location back to model coordinates, then - ! reset transformation and eliminate accumulated roundoff error. + ! Transform cell coordinates back to model coordinates call particle%transform(xOrigin, yOrigin, zOrigin, & sinrot, cosrot, invert=.true.) - call particle%transform(reset=.true.) + call particle%reset_transform() end select end subroutine apply_mcp diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 index 8fe7e41525f..c5f61f2c221 100644 --- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 +++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 @@ -37,7 +37,7 @@ subroutine create_method_cell_quad(method) allocate (method) call create_cell_rect_quad(cell) method%cell => cell - method%type => method%cell%type + method%name => method%cell%type method%delegates = .true. call create_subcell_rect(subcell) method%subcell => subcell @@ -46,7 +46,7 @@ end subroutine create_method_cell_quad !> @brief Deallocate the Pollock quad-refined cell method subroutine deallocate (this) class(MethodCellPollockQuadType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine deallocate !> @brief Load subcell into tracking method @@ -208,19 +208,12 @@ subroutine apply_mcpq(this, particle, tmax) select type (cell => this%cell) type is (CellRectQuadType) - ! Update particle state, return early if done advancing - call this%update(particle, cell%defn) + ! Check termination/reporting conditions + call this%check(particle, cell%defn) if (.not. particle%advancing) return - ! If the particle is above the top of the cell (presumed water table) - ! pass it vertically and instantaneously to the top - if (particle%z > cell%defn%top) then - particle%z = cell%defn%top - call this%save(particle, reason=1) ! reason=1: cell transition - end if - - ! Transform particle location into local cell coordinates - ! (translated and rotated but not scaled relative to model). + ! Transform model coordinates to local cell coordinates + ! (translated/rotated but not scaled relative to model) xOrigin = cell%xOrigin yOrigin = cell%yOrigin zOrigin = cell%zOrigin @@ -229,14 +222,13 @@ subroutine apply_mcpq(this, particle, tmax) call particle%transform(xOrigin, yOrigin, zOrigin, & sinrot, cosrot) - ! Track the particle across the cell. + ! Track the particle over the cell call this%track(particle, 2, tmax) - ! Transform particle location back to model coordinates, then - ! reset transformation and eliminate accumulated roundoff error. + ! Transform cell coordinates back to model coordinates call particle%transform(xOrigin, yOrigin, zOrigin, & sinrot, cosrot, invert=.true.) - call particle%transform(reset=.true.) + call particle%reset_transform() end select end subroutine apply_mcpq @@ -299,7 +291,7 @@ subroutine load_subcell(this, particle, subcell) qintl2 = cell%qintl(isc + 1) qextl1 = cell%qextl1(isc) qextl2 = cell%qextl2(isc) - ! + subcell%dx = dx subcell%dy = dy subcell%dz = dz diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90 index fac00b00845..a04a65d231f 100644 --- a/src/Solution/ParticleTracker/MethodCellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodCellTernary.f90 @@ -62,7 +62,7 @@ subroutine create_method_cell_ternary(method) allocate (method) call create_cell_poly(cell) method%cell => cell - method%type => method%cell%type + method%name => method%cell%type method%delegates = .true. call create_subcell_tri(subcell) method%subcell => subcell @@ -71,7 +71,7 @@ end subroutine create_method_cell_ternary !> @brief Destroy the tracking method subroutine destroy_mct(this) class(MethodCellTernaryType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine destroy_mct !> @brief Load subcell into tracking method @@ -156,22 +156,16 @@ subroutine apply_mct(this, particle, tmax) ! local integer(I4B) :: i - ! Update particle state, return early if done advancing - call this%update(particle, this%cell%defn) - if (.not. particle%advancing) return - - ! If the particle is above the top of the cell (presumed water table) - ! pass it vertically and instantaneously to the top - if (particle%z > this%cell%defn%top) then - particle%z = this%cell%defn%top - call this%save(particle, reason=1) - end if - ! (Re)allocate type-bound arrays select type (cell => this%cell) type is (CellPolyType) + ! Check termination/reporting conditions + call this%check(particle, this%cell%defn) + if (.not. particle%advancing) return + ! Number of vertices this%nverts = cell%defn%npolyverts + ! (Re)allocate type-bound arrays if (allocated(this%xvert)) then deallocate (this%xvert) @@ -185,6 +179,7 @@ subroutine apply_mct(this, particle, tmax) deallocate (this%xvertnext) deallocate (this%yvertnext) end if + allocate (this%xvert(this%nverts)) allocate (this%yvert(this%nverts)) allocate (this%vne(this%nverts)) @@ -195,15 +190,18 @@ subroutine apply_mct(this, particle, tmax) allocate (this%iprev(this%nverts)) allocate (this%xvertnext(this%nverts)) allocate (this%yvertnext(this%nverts)) + ! Cell vertices do i = 1, this%nverts this%xvert(i) = cell%defn%polyvert(1, i) this%yvert(i) = cell%defn%polyvert(2, i) end do + ! Top, bottom, and thickness this%ztop = cell%defn%top this%zbot = cell%defn%bot this%dz = this%ztop - this%zbot + ! Shifted arrays do i = 1, this%nverts this%iprev(i) = i diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90 index ad427c9e411..6167837dc19 100644 --- a/src/Solution/ParticleTracker/MethodDis.f90 +++ b/src/Solution/ParticleTracker/MethodDis.f90 @@ -1,6 +1,6 @@ module MethodDisModule - use KindModule, only: DP, I4B + use KindModule, only: DP, I4B, LGP use ConstantsModule, only: DONE, DZERO use MethodModule, only: MethodType use MethodCellPoolModule @@ -11,6 +11,7 @@ module MethodDisModule use PrtFmiModule, only: PrtFmiType use DisModule, only: DisType use GeomUtilModule, only: get_ijk, get_jk + use MathUtilModule, only: is_close implicit none private @@ -40,23 +41,23 @@ module MethodDisModule !> @brief Create a new structured grid (DIS) tracking method subroutine create_method_dis(method) - ! -- dummy + ! dummy type(MethodDisType), pointer :: method - ! -- local + ! local type(CellRectType), pointer :: cell allocate (method) - allocate (method%type) + allocate (method%name) call create_cell_rect(cell) method%cell => cell - method%type = "dis" + method%name = "dis" method%delegates = .true. end subroutine create_method_dis !> @brief Destructor the tracking method subroutine deallocate (this) class(MethodDisType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine deallocate subroutine load_cell(this, ic, cell) @@ -81,11 +82,14 @@ subroutine load_cell(this, ic, cell) select type (dis => this%fmi%dis) type is (DisType) icu = dis%get_nodeuser(ic) + call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, & irow, jcol, klay) + dx = dis%delr(jcol) dy = dis%delc(irow) dz = cell%defn%top - cell%defn%bot + cell%dx = dx cell%dy = dy cell%dz = dz @@ -95,41 +99,54 @@ subroutine load_cell(this, ic, cell) cell%yOrigin = cell%defn%polyvert(2, 1) cell%zOrigin = cell%defn%bot cell%ipvOrigin = 1 - areax = dy * dz - areay = dx * dz - areaz = dx * dy + factor = DONE / cell%defn%retfactor factor = factor / cell%defn%porosity + + areaz = dx * dy + term = factor / areaz + + cell%vz1 = cell%defn%faceflow(6) * term + cell%vz2 = -cell%defn%faceflow(7) * term + + if (this%fmi%ibdgwfsat0(ic) == 0) then + cell%vx1 = DZERO + cell%vx2 = DZERO + cell%vy1 = DZERO + cell%vy2 = DZERO + cell%vz1 = DZERO + cell%vz2 = DZERO + return + end if + + areax = dy * dz term = factor / areax cell%vx1 = cell%defn%faceflow(1) * term cell%vx2 = -cell%defn%faceflow(3) * term + + areay = dx * dz term = factor / areay cell%vy1 = cell%defn%faceflow(4) * term cell%vy2 = -cell%defn%faceflow(2) * term - term = factor / areaz - cell%vz1 = cell%defn%faceflow(6) * term - cell%vz2 = -cell%defn%faceflow(7) * term + end select end subroutine load_cell !> @brief Load the cell geometry and tracking method subroutine load_dis(this, particle, next_level, submethod) - ! -- dummy + ! dummy class(MethodDisType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle integer(I4B), intent(in) :: next_level class(MethodType), pointer, intent(inout) :: submethod - ! -- local + ! local integer(I4B) :: ic select type (cell => this%cell) type is (CellRectType) - ! -- Load cell definition and geometry ic = particle%idomain(next_level) call this%load_celldefn(ic, cell%defn) call this%load_cell(ic, cell) - - ! -- If cell is active but dry, Initialize instant pass-to-bottom method if (this%fmi%ibdgwfsat0(ic) == 0) then call method_cell_ptb%init( & fmi=this%fmi, & @@ -138,7 +155,6 @@ subroutine load_dis(this, particle, next_level, submethod) tracktimes=this%tracktimes) submethod => method_cell_ptb else - ! -- Otherwise initialize Pollock's method call method_cell_plck%init( & fmi=this%fmi, & cell=this%cell, & @@ -176,7 +192,7 @@ subroutine load_particle(this, cell, particle) select type (dis => this%fmi%dis) type is (DisType) - ! compute and set reduced/user node numbers and layer + ! compute reduced/user node numbers and layer inface = particle%iboundary(2) inbr = cell%defn%facenbr(inface) idiag = this%fmi%dis%con%ia(cell%defn%icell) @@ -185,6 +201,25 @@ subroutine load_particle(this, cell, particle) icu = dis%get_nodeuser(ic) call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, & irow, icol, ilay) + + ! if returning to a cell through the bottom + ! face after previously leaving it through + ! that same face, we've entered a cycle + ! as can occur e.g. in wells. terminate + ! in the previous cell. + if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then + particle%advancing = .false. + particle%idomain(2) = particle%icp + particle%istatus = 2 + particle%izone = particle%izp + call this%save(particle, reason=3) + return + else + particle%icp = particle%idomain(2) + particle%izp = particle%izone + end if + + ! update node numbers and layer particle%idomain(2) = ic particle%icu = icu particle%ilay = ilay @@ -218,7 +253,6 @@ subroutine load_particle(this, cell, particle) end if particle%z = z end select - end subroutine load_particle !> @brief Update cell-cell flows of particle mass. @@ -239,13 +273,12 @@ subroutine update_flowja(this, cell, particle) idiag = this%fmi%dis%con%ia(cell%defn%icell) ipos = idiag + inbr - ! -- leaving old cell + ! leaving old cell this%flowja(ipos) = this%flowja(ipos) - DONE - ! -- entering new cell + ! entering new cell this%flowja(this%fmi%dis%con%isym(ipos)) & = this%flowja(this%fmi%dis%con%isym(ipos)) + DONE - end subroutine update_flowja !> @brief Pass a particle to the next cell, if there is one @@ -267,15 +300,17 @@ subroutine pass_dis(this, particle) particle%advancing = .false. call this%save(particle, reason=3) else - ! Otherwise, load cell properties into the - ! particle and update intercell mass flows + ! Update old to new cell properties call this%load_particle(cell, particle) + if (.not. particle%advancing) return + + ! Update intercell mass flows call this%update_flowja(cell, particle) end if end select end subroutine pass_dis - !> @brief Apply the method to a particle + !> @brief Apply the structured tracking method to a particle. subroutine apply_dis(this, particle, tmax) class(MethodDisType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle @@ -299,36 +334,35 @@ end function get_top !> @brief Loads cell definition from the grid subroutine load_celldefn(this, ic, defn) - ! -- dummy + ! dummy class(MethodDisType), intent(inout) :: this integer(I4B), intent(in) :: ic type(CellDefnType), pointer, intent(inout) :: defn - ! -- Load basic cell properties + ! Load basic cell properties call this%load_properties(ic, defn) - ! -- Load cell polygon vertices + ! Load cell polygon vertices call this%fmi%dis%get_polyverts( & defn%icell, & defn%polyvert, & closed=.true.) - - ! -- Load cell face neighbors call this%load_neighbors(defn) - ! -- Load 180 degree face indicators + ! Load 180 degree face indicators defn%ispv180(1:defn%npolyverts + 1) = .false. - ! -- Load face flows (assumes face neighbors already loaded) call this%load_flows(defn) end subroutine load_celldefn subroutine load_properties(this, ic, defn) - ! -- dummy + ! dummy class(MethodDisType), intent(inout) :: this integer(I4B), intent(in) :: ic type(CellDefnType), pointer, intent(inout) :: defn + ! local + integer(I4B) :: irow, icol, ilay, icu defn%icell = ic defn%npolyverts = 4 ! rectangular cell always has 4 vertices @@ -341,19 +375,24 @@ subroutine load_properties(this, ic, defn) defn%sat = this%fmi%gwfsat(ic) defn%porosity = this%porosity(ic) defn%retfactor = this%retfactor(ic) + select type (dis => this%fmi%dis) + type is (DisType) + icu = dis%get_nodeuser(ic) + call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay) + defn%ilay = ilay + end select defn%izone = this%izone(ic) defn%can_be_rect = .true. defn%can_be_quad = .false. - end subroutine load_properties !> @brief Loads face neighbors to cell definition from the grid. !! Assumes cell index and number of vertices are already loaded. subroutine load_neighbors(this, defn) - ! -- dummy + ! dummy class(MethodDisType), intent(inout) :: this type(CellDefnType), pointer, intent(inout) :: defn - ! -- local + ! local integer(I4B) :: ic1 integer(I4B) :: ic2 integer(I4B) :: icu1 @@ -371,7 +410,7 @@ subroutine load_neighbors(this, defn) select type (dis => this%fmi%dis) type is (DisType) - ! -- Load face neighbors + ! Load face neighbors defn%facenbr = 0 ic1 = defn%icell icu1 = dis%get_nodeuser(ic1) @@ -387,7 +426,7 @@ subroutine load_neighbors(this, defn) call get_ijk(icu2, dis%nrow, dis%ncol, dis%nlay, & irow2, jcol2, klay2) if (klay2 == klay1) then - ! -- Edge (polygon) face neighbor + ! Edge (polygon) face neighbor if (irow2 > irow1) then ! Neighbor to the S iedgeface = 4 @@ -403,15 +442,15 @@ subroutine load_neighbors(this, defn) end if defn%facenbr(iedgeface) = int(iloc, 1) else if (klay2 > klay1) then - ! -- Bottom face neighbor + ! Bottom face neighbor defn%facenbr(defn%npolyverts + 2) = int(iloc, 1) else - ! -- Top face neighbor + ! Top face neighbor defn%facenbr(defn%npolyverts + 3) = int(iloc, 1) end if end do end select - ! -- List of edge (polygon) faces wraps around + ! List of edge (polygon) faces wraps around defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1) end subroutine load_neighbors @@ -442,7 +481,6 @@ subroutine load_flows(this, defn) else defn%iweaksink = 0 end if - end subroutine load_flows subroutine load_face_flows_to_defn(this, defn) @@ -467,10 +505,10 @@ end subroutine load_face_flows_to_defn !> @brief Add boundary flows to the cell definition faceflow array. !! Assumes cell index and number of vertices are already loaded. subroutine load_boundary_flows_to_defn(this, defn) - ! -- dummy + ! dummy class(MethodDisType), intent(inout) :: this type(CellDefnType), intent(inout) :: defn - ! -- local + ! local integer(I4B) :: ioffset ioffset = (defn%icell - 1) * MAX_POLY_CELLS diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90 index 4303594af9f..de65bbbf134 100644 --- a/src/Solution/ParticleTracker/MethodDisv.f90 +++ b/src/Solution/ParticleTracker/MethodDisv.f90 @@ -46,16 +46,16 @@ module MethodDisvModule !> @brief Create a new vertex grid (DISV) tracking method subroutine create_method_disv(method) - ! -- dummy + ! dummy type(MethodDisvType), pointer :: method - ! -- local + ! local type(CellPolyType), pointer :: cell allocate (method) - allocate (method%type) + allocate (method%name) call create_cell_poly(cell) method%cell => cell - method%type = "disv" + method%name = "disv" method%delegates = .true. call create_defn(method%neighbor) end subroutine create_method_disv @@ -63,7 +63,7 @@ end subroutine create_method_disv !> @brief Destroy the tracking method subroutine deallocate (this) class(MethodDisvType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine deallocate !> @brief Load the cell and the tracking method @@ -72,12 +72,12 @@ subroutine load_disv(this, particle, next_level, submethod) use CellRectModule use CellRectQuadModule use CellUtilModule - ! -- dummy + ! dummy class(MethodDisvType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle integer(I4B), intent(in) :: next_level class(MethodType), pointer, intent(inout) :: submethod - ! -- local + ! local integer(I4B) :: ic class(CellType), pointer :: base type(CellRectType), pointer :: rect @@ -85,53 +85,47 @@ subroutine load_disv(this, particle, next_level, submethod) select type (cell => this%cell) type is (CellPolyType) - ! load cell definition ic = particle%idomain(next_level) call this%load_cell_defn(ic, cell%defn) if (this%fmi%ibdgwfsat0(ic) == 0) then - ! -- Cell is active but dry, so select and initialize pass-to-bottom - ! -- cell method and set cell method pointer call method_cell_ptb%init( & fmi=this%fmi, & cell=this%cell, & trackctl=this%trackctl, & tracktimes=this%tracktimes) submethod => method_cell_ptb + else if (particle%ifrctrn > 0) then + call method_cell_tern%init( & + fmi=this%fmi, & + cell=this%cell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) + submethod => method_cell_tern + else if (cell%defn%can_be_rect) then + call cell_poly_to_rect(cell, rect) + base => rect + call method_cell_plck%init( & + fmi=this%fmi, & + cell=base, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) + submethod => method_cell_plck + else if (cell%defn%can_be_quad) then + call cell_poly_to_quad(cell, quad) + base => quad + call method_cell_quad%init( & + fmi=this%fmi, & + cell=base, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) + submethod => method_cell_quad else - ! -- Select and initialize cell method and set cell method pointer - if (particle%ifrctrn > 0) then - call method_cell_tern%init( & - fmi=this%fmi, & - cell=this%cell, & - trackctl=this%trackctl, & - tracktimes=this%tracktimes) - submethod => method_cell_tern - else if (cell%defn%can_be_rect) then - call cell_poly_to_rect(cell, rect) - base => rect - call method_cell_plck%init( & - fmi=this%fmi, & - cell=base, & - trackctl=this%trackctl, & - tracktimes=this%tracktimes) - submethod => method_cell_plck - else if (cell%defn%can_be_quad) then - call cell_poly_to_quad(cell, quad) - base => quad - call method_cell_quad%init( & - fmi=this%fmi, & - cell=base, & - trackctl=this%trackctl, & - tracktimes=this%tracktimes) - submethod => method_cell_quad - else - call method_cell_tern%init( & - fmi=this%fmi, & - cell=this%cell, & - trackctl=this%trackctl, & - tracktimes=this%tracktimes) - submethod => method_cell_tern - end if + call method_cell_tern%init( & + fmi=this%fmi, & + cell=this%cell, & + trackctl=this%trackctl, & + tracktimes=this%tracktimes) + submethod => method_cell_tern end if end select end subroutine load_disv @@ -156,7 +150,6 @@ subroutine load_particle(this, cell, particle) select type (dis => this%fmi%dis) type is (DisvType) - ! compute and set reduced/user node numbers and layer inface = particle%iboundary(2) idiag = dis%con%ia(cell%defn%icell) inbr = cell%defn%facenbr(inface) @@ -164,13 +157,31 @@ subroutine load_particle(this, cell, particle) ic = dis%con%ja(ipos) icu = dis%get_nodeuser(ic) call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay) + + ! if returning to a cell through the bottom + ! face after previously leaving it through + ! that same face, we've entered a cycle + ! as can occur e.g. in wells. terminate + ! in the previous cell. + if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then + particle%advancing = .false. + particle%idomain(2) = particle%icp + particle%istatus = 2 + particle%izone = particle%izp + call this%save(particle, reason=3) + return + else + particle%icp = particle%idomain(2) + particle%izp = particle%izone + end if + particle%idomain(2) = ic particle%icu = icu particle%ilay = ilay - ! map/set particle entry face and z coordinate z = particle%z call this%map_neighbor(cell%defn, inface, z) + particle%iboundary(2) = inface particle%idomain(3:) = 0 particle%iboundary(3:) = 0 @@ -221,9 +232,11 @@ subroutine pass_disv(this, particle) particle%advancing = .false. call this%save(particle, reason=3) else - ! Otherwise, load cell properties into the - ! particle and update intercell mass flows + ! Update old to new cell properties call this%load_particle(cell, particle) + if (.not. particle%advancing) return + + ! Update intercell mass flows call this%update_flowja(cell, particle) end if end select @@ -252,14 +265,14 @@ subroutine map_neighbor(this, defn, inface, z) real(DP) :: bot real(DP) :: sat - ! -- Map to shared cell face of neighbor + ! Map to shared cell face of neighbor inbr = defn%facenbr(inface) if (inbr .eq. 0) then - ! -- Exterior face; no neighbor to map to + ! Exterior face; no neighbor to map to ! todo AMP: reconsider when multiple models allowed inface = -1 else - ! -- Load definition for neighbor cell (neighbor with shared face) + ! Load definition for neighbor cell (neighbor with shared face) icin = defn%icell j = this%fmi%dis%con%ia(icin) ic = this%fmi%dis%con%ja(j + inbr) @@ -268,13 +281,13 @@ subroutine map_neighbor(this, defn, inface, z) npolyvertsin = defn%npolyverts npolyverts = this%neighbor%npolyverts if (inface .eq. npolyvertsin + 2) then - ! -- Exits through bot, enters through top + ! Exits through bot, enters through top inface = npolyverts + 3 else if (inface .eq. npolyvertsin + 3) then - ! -- Exits through top, enters through bot + ! Exits through top, enters through bot inface = npolyverts + 2 else - ! -- Exits and enters through shared polygon face + ! Exits and enters through shared polygon face j = this%fmi%dis%con%ia(ic) do m = 1, npolyverts + 3 inbrnbr = this%neighbor%facenbr(m) @@ -283,7 +296,7 @@ subroutine map_neighbor(this, defn, inface, z) exit end if end do - ! -- Map z between cells + ! Map z between cells topfrom = defn%top botfrom = defn%bot zrel = (z - botfrom) / (topfrom - botfrom) @@ -295,34 +308,26 @@ subroutine map_neighbor(this, defn, inface, z) end if end subroutine map_neighbor - !> @brief Apply the DISV-grid method + !> @brief Apply the DISV tracking method to a particle. subroutine apply_disv(this, particle, tmax) class(MethodDisvType), intent(inout) :: this type(ParticleType), pointer, intent(inout) :: particle real(DP), intent(in) :: tmax + call this%track(particle, 1, tmax) end subroutine apply_disv !> @brief Load cell definition from the grid subroutine load_cell_defn(this, ic, defn) - ! -- dummy + ! dummy class(MethodDisvType), intent(inout) :: this integer(I4B), intent(in) :: ic type(CellDefnType), pointer, intent(inout) :: defn - ! -- Load basic cell properties call this%load_properties(ic, defn) - - ! -- Load polygon vertices call this%load_polygon(defn) - - ! -- Load face neighbors call this%load_neighbors(defn) - - ! -- Load 180-degree indicator call this%load_indicators(defn) - - ! -- Load flows (assumes face neighbors already loaded) call this%load_flows(defn) end subroutine load_cell_defn @@ -336,6 +341,7 @@ subroutine load_properties(this, ic, defn) real(DP) :: top real(DP) :: bot real(DP) :: sat + integer(I4B) :: icu, icpl, ilay defn%icell = ic defn%iatop = get_iatop(this%fmi%dis%get_ncpl(), & @@ -350,7 +356,12 @@ subroutine load_properties(this, ic, defn) defn%porosity = this%porosity(ic) defn%retfactor = this%retfactor(ic) defn%izone = this%izone(ic) - + select type (dis => this%fmi%dis) + type is (DisvType) + icu = dis%get_nodeuser(ic) + call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay) + defn%ilay = ilay + end select end subroutine load_properties subroutine load_polygon(this, defn) @@ -363,16 +374,15 @@ subroutine load_polygon(this, defn) defn%polyvert, & closed=.true.) defn%npolyverts = size(defn%polyvert, dim=2) - 1 - end subroutine load_polygon !> @brief Loads face neighbors to cell definition from the grid !! Assumes cell index and number of vertices are already loaded. subroutine load_neighbors(this, defn) - ! -- dummy + ! dummy class(MethodDisvType), intent(inout) :: this type(CellDefnType), pointer, intent(inout) :: defn - ! -- local + ! local integer(I4B) :: ic1 integer(I4B) :: ic2 integer(I4B) :: icu1 @@ -392,14 +402,14 @@ subroutine load_neighbors(this, defn) integer(I4B) :: nfaces integer(I4B) :: nslots - ! -- expand facenbr array if needed + ! expand facenbr array if needed nfaces = defn%npolyverts + 3 nslots = size(defn%facenbr) if (nslots < nfaces) call ExpandArray(defn%facenbr, nfaces - nslots) select type (dis => this%fmi%dis) type is (DisvType) - ! -- Load face neighbors. + ! Load face neighbors. defn%facenbr = 0 ic1 = defn%icell icu1 = dis%get_nodeuser(ic1) @@ -420,14 +430,14 @@ subroutine load_neighbors(this, defn) dis%javert(istart2:istop2), & isharedface) if (isharedface /= 0) then - ! -- Edge (polygon) face neighbor + ! Edge (polygon) face neighbor defn%facenbr(isharedface) = int(iloc, 1) else if (k2 > k1) then - ! -- Bottom face neighbor + ! Bottom face neighbor defn%facenbr(defn%npolyverts + 2) = int(iloc, 1) else if (k2 < k1) then - ! -- Top face neighbor + ! Top face neighbor defn%facenbr(defn%npolyverts + 3) = int(iloc, 1) else call pstop(1, "k2 should be <> k1, since no shared edge face") @@ -435,9 +445,8 @@ subroutine load_neighbors(this, defn) end if end do end select - ! -- List of edge (polygon) faces wraps around + ! List of edge (polygon) faces wraps around defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1) - end subroutine load_neighbors !> @brief Load flows into the cell definition. @@ -476,7 +485,6 @@ subroutine load_flows(this, defn) else defn%iweaksink = 0 end if - end subroutine load_flows subroutine load_face_flows_to_defn_poly(this, defn) @@ -501,10 +509,10 @@ end subroutine load_face_flows_to_defn_poly !> @brief Load boundary flows from the grid into a rectangular cell. !! Assumes cell index and number of vertices are already loaded. subroutine load_boundary_flows_to_defn_rect(this, defn) - ! -- dummy + ! dummy class(MethodDisvType), intent(inout) :: this type(CellDefnType), intent(inout) :: defn - ! -- local + ! local integer(I4B) :: ioffset ! assignment of BoundaryFlows to faceflow below assumes clockwise @@ -523,16 +531,15 @@ subroutine load_boundary_flows_to_defn_rect(this, defn) this%fmi%BoundaryFlows(ioffset + 9) defn%faceflow(7) = defn%faceflow(7) + & this%fmi%BoundaryFlows(ioffset + 10) - end subroutine load_boundary_flows_to_defn_rect !> @brief Load boundary flows from the grid into rectangular quadcell. !! Assumes cell index and number of vertices are already loaded. subroutine load_boundary_flows_to_defn_rect_quad(this, defn) - ! -- dummy + ! dummy class(MethodDisvType), intent(inout) :: this type(CellDefnType), intent(inout) :: defn - ! -- local + ! local integer(I4B) :: m integer(I4B) :: n integer(I4B) :: nn @@ -546,7 +553,7 @@ subroutine load_boundary_flows_to_defn_rect_quad(this, defn) ioffset = (defn%icell - 1) * 10 - ! -- Polygon faces in positions 1 through npolyverts + ! Polygon faces in positions 1 through npolyverts do n = 1, 4 if (n .eq. 2) then nbf = 4 @@ -569,23 +576,23 @@ subroutine load_boundary_flows_to_defn_rect_quad(this, defn) if (m2 .lt. m1) m2 = m2 + defn%npolyverts mdiff = m2 - m1 if (mdiff .eq. 1) then - ! -- Assign BoundaryFlow to corresponding polygon face + ! Assign BoundaryFlow to corresponding polygon face defn%faceflow(m1) = defn%faceflow(m1) + qbf else - ! -- Split BoundaryFlow between two faces on quad-refined edge + ! Split BoundaryFlow between two faces on quad-refined edge qbf = 5d-1 * qbf defn%faceflow(m1) = defn%faceflow(m1) + qbf defn%faceflow(m1 + 1) = defn%faceflow(m1 + 1) + qbf end if end do - ! -- Wrap around to 1 in position npolyverts+1 + ! Wrap around to 1 in position npolyverts+1 m = defn%npolyverts + 1 defn%faceflow(m) = defn%faceflow(1) - ! -- Bottom in position npolyverts+2 + ! Bottom in position npolyverts+2 m = m + 1 defn%faceflow(m) = defn%faceflow(m) + & this%fmi%BoundaryFlows(ioffset + 9) - ! -- Top in position npolyverts+3 + ! Top in position npolyverts+3 m = m + 1 defn%faceflow(m) = defn%faceflow(m) + & this%fmi%BoundaryFlows(ioffset + 10) @@ -595,10 +602,10 @@ end subroutine load_boundary_flows_to_defn_rect_quad !> @brief Load boundary flows from the grid into a polygonal cell. !! Assumes cell index and number of vertices are already loaded. subroutine load_boundary_flows_to_defn_poly(this, defn) - ! -- dummy + ! dummy class(MethodDisvType), intent(inout) :: this type(CellDefnType), intent(inout) :: defn - ! -- local + ! local integer(I4B) :: ic integer(I4B) :: npolyverts integer(I4B) :: ioffset @@ -624,13 +631,14 @@ subroutine load_boundary_flows_to_defn_poly(this, defn) end subroutine load_boundary_flows_to_defn_poly !> @brief Load 180-degree vertex indicator array and set flags - !! indicating how cell can be represented. - !! Assumes cell index and number of vertices are already loaded. + !! indicating how cell can be represented. Assumes cell index + !! and number of vertices are already loaded. + !< subroutine load_indicators(this, defn) - ! -- dummy + ! dummy class(MethodDisvType), intent(inout) :: this type(CellDefnType), pointer, intent(inout) :: defn - ! -- local + ! local integer(I4B) :: npolyverts integer(I4B) :: m integer(I4B) :: m0 @@ -655,12 +663,9 @@ subroutine load_indicators(this, defn) ic = defn%icell npolyverts = defn%npolyverts - ! -- expand ispv180 array if needed if (size(defn%ispv180) < npolyverts + 3) & call ExpandArray(defn%ispv180, npolyverts + 1) - ! -- Load 180-degree indicator. - ! -- Also, set flags that indicate how cell can be represented. defn%ispv180(1:npolyverts + 1) = .false. defn%can_be_rect = .false. defn%can_be_quad = .false. @@ -668,6 +673,7 @@ subroutine load_indicators(this, defn) num90 = 0 num180 = 0 last180 = .false. + ! assumes non-self-intersecting polygon do m = 1, npolyverts m1 = m @@ -708,9 +714,10 @@ subroutine load_indicators(this, defn) end if end do - ! -- List of 180-degree indicators wraps around for convenience + ! List of 180-degree indicators wraps around for convenience defn%ispv180(npolyverts + 1) = defn%ispv180(1) + ! Set rect/quad flags if (num90 .eq. 4) then if (num180 .eq. 0) then defn%can_be_rect = .true. @@ -718,7 +725,6 @@ subroutine load_indicators(this, defn) defn%can_be_quad = .true. end if end if - end subroutine load_indicators end module MethodDisvModule diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 index de9157b8b9e..1cbd0cce45b 100644 --- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90 @@ -36,14 +36,14 @@ subroutine create_method_subcell_pollock(method) allocate (method) call create_subcell_rect(subcell) method%subcell => subcell - method%type => method%subcell%type + method%name => method%subcell%type method%delegates = .false. end subroutine create_method_subcell_pollock !> @brief Deallocate the Pollock's subcell method subroutine deallocate (this) class(MethodSubcellPollockType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine deallocate !> @brief Apply Pollock's method to a rectangular subcell @@ -73,6 +73,7 @@ subroutine apply_msp(this, particle, tmax) call particle%transform(xOrigin, yOrigin) call this%track_subcell(subcell, particle, tmax) call particle%transform(xOrigin, yOrigin, invert=.true.) + call particle%reset_transform() end select end subroutine apply_msp diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 index 057eca00050..97637b4c59b 100644 --- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 +++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90 @@ -38,14 +38,14 @@ subroutine create_method_subcell_ternary(method) allocate (method) call create_subcell_tri(subcell) method%subcell => subcell - method%type => method%subcell%type + method%name => method%subcell%type method%delegates = .false. end subroutine create_method_subcell_ternary !> @brief Deallocate the ternary subcell tracking method. subroutine deallocate (this) class(MethodSubcellTernaryType), intent(inout) :: this - deallocate (this%type) + deallocate (this%name) end subroutine deallocate !> @brief Apply the ternary subcell tracking method. diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90 index 3bb1471345c..347ae4b1a2b 100644 --- a/src/Solution/ParticleTracker/Particle.f90 +++ b/src/Solution/ParticleTracker/Particle.f90 @@ -40,12 +40,15 @@ module ParticleModule ! stop criteria integer(I4B), public :: istopweaksink !< weak sink option (0: do not stop, 1: stop) integer(I4B), public :: istopzone !< stop zone number + integer(I4B), public :: idrymeth !< dry tracking method ! state - integer(I4B), allocatable, public :: idomain(:) !< tracking domain hierarchy + integer(I4B), allocatable, public :: idomain(:) !< tracking domain hierarchy ! TODO: rename to itdomain? idomain integer(I4B), allocatable, public :: iboundary(:) !< tracking domain boundaries - integer(I4B), public :: icu !< user cell (node) number + integer(I4B), public :: icp !< previous cell number (reduced) + integer(I4B), public :: icu !< user cell number integer(I4B), public :: ilay !< grid layer - integer(I4B), public :: izone !< zone number + integer(I4B), public :: izone !< current zone number + integer(I4B), public :: izp !< previous zone number integer(I4B), public :: istatus !< tracking status real(DP), public :: x !< x coordinate real(DP), public :: y !< y coordinate @@ -68,6 +71,7 @@ module ParticleModule procedure, public :: get_model_coords procedure, public :: load_particle procedure, public :: transform => transform_coords + procedure, public :: reset_transform end type ParticleType !> @brief Structure of arrays to store particles. @@ -81,12 +85,15 @@ module ParticleModule ! stopping criteria integer(I4B), dimension(:), pointer, public, contiguous :: istopweaksink !< weak sink option: 0 = do not stop, 1 = stop integer(I4B), dimension(:), pointer, public, contiguous :: istopzone !< stop zone number + integer(I4B), dimension(:), pointer, public, contiguous :: idrymeth !< stop in dry cells ! state integer(I4B), dimension(:, :), pointer, public, contiguous :: idomain !< array of indices for domains in the tracking domain hierarchy integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< array of indices for tracking domain boundaries - integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user, not reduced) + integer(I4B), dimension(:), pointer, public, contiguous :: icp !< previous cell number (reduced) + integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user) integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number + integer(I4B), dimension(:), pointer, public, contiguous :: izp !< previous zone number integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status real(DP), dimension(:), pointer, public, contiguous :: x !< model x coord of particle real(DP), dimension(:), pointer, public, contiguous :: y !< model y coord of particle @@ -126,9 +133,11 @@ subroutine allocate_particle_store(this, np, mempath) call mem_allocate(this%irpt, np, 'PLIRPT', mempath) call mem_allocate(this%iprp, np, 'PLIPRP', mempath) call mem_allocate(this%name, LENBOUNDNAME, np, 'PLNAME', mempath) + call mem_allocate(this%icp, np, 'PLICP', mempath) call mem_allocate(this%icu, np, 'PLICU', mempath) call mem_allocate(this%ilay, np, 'PLILAY', mempath) call mem_allocate(this%izone, np, 'PLIZONE', mempath) + call mem_allocate(this%izp, np, 'PLIZP', mempath) call mem_allocate(this%istatus, np, 'PLISTATUS', mempath) call mem_allocate(this%x, np, 'PLX', mempath) call mem_allocate(this%y, np, 'PLY', mempath) @@ -138,6 +147,7 @@ subroutine allocate_particle_store(this, np, mempath) call mem_allocate(this%ttrack, np, 'PLTTRACK', mempath) call mem_allocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath) call mem_allocate(this%istopzone, np, 'PLISTOPZONE', mempath) + call mem_allocate(this%idrymeth, np, 'PLIDRYMETH', mempath) call mem_allocate(this%ifrctrn, np, 'PLIFRCTRN', mempath) call mem_allocate(this%iexmeth, np, 'PLIEXMETH', mempath) call mem_allocate(this%extol, np, 'PLEXTOL', mempath) @@ -155,9 +165,11 @@ subroutine deallocate (this, mempath) call mem_deallocate(this%iprp, 'PLIPRP', mempath) call mem_deallocate(this%irpt, 'PLIRPT', mempath) call mem_deallocate(this%name, 'PLNAME', mempath) + call mem_deallocate(this%icp, 'PLICP', mempath) call mem_deallocate(this%icu, 'PLICU', mempath) call mem_deallocate(this%ilay, 'PLILAY', mempath) call mem_deallocate(this%izone, 'PLIZONE', mempath) + call mem_deallocate(this%izp, 'PLIZP', mempath) call mem_deallocate(this%istatus, 'PLISTATUS', mempath) call mem_deallocate(this%x, 'PLX', mempath) call mem_deallocate(this%y, 'PLY', mempath) @@ -167,6 +179,7 @@ subroutine deallocate (this, mempath) call mem_deallocate(this%ttrack, 'PLTTRACK', mempath) call mem_deallocate(this%istopweaksink, 'PLISTOPWEAKSINK', mempath) call mem_deallocate(this%istopzone, 'PLISTOPZONE', mempath) + call mem_deallocate(this%idrymeth, 'PLIDRYMETH', mempath) call mem_deallocate(this%ifrctrn, 'PLIFRCTRN', mempath) call mem_deallocate(this%iexmeth, 'PLIEXMETH', mempath) call mem_deallocate(this%extol, 'PLEXTOL', mempath) @@ -177,7 +190,7 @@ end subroutine deallocate !> @brief Reallocate particle arrays subroutine resize(this, np, mempath) - ! -- dummy + ! dummy class(ParticleStoreType), intent(inout) :: this !< particle store integer(I4B), intent(in) :: np !< number of particles character(*), intent(in) :: mempath !< path to memory @@ -187,9 +200,11 @@ subroutine resize(this, np, mempath) call mem_reallocate(this%iprp, np, 'PLIPRP', mempath) call mem_reallocate(this%irpt, np, 'PLIRPT', mempath) call mem_reallocate(this%name, LENBOUNDNAME, np, 'PLNAME', mempath) + call mem_reallocate(this%icp, np, 'PLICP', mempath) call mem_reallocate(this%icu, np, 'PLICU', mempath) call mem_reallocate(this%ilay, np, 'PLILAY', mempath) call mem_reallocate(this%izone, np, 'PLIZONE', mempath) + call mem_reallocate(this%izp, np, 'PLIZP', mempath) call mem_reallocate(this%istatus, np, 'PLISTATUS', mempath) call mem_reallocate(this%x, np, 'PLX', mempath) call mem_reallocate(this%y, np, 'PLY', mempath) @@ -199,6 +214,7 @@ subroutine resize(this, np, mempath) call mem_reallocate(this%ttrack, np, 'PLTTRACK', mempath) call mem_reallocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath) call mem_reallocate(this%istopzone, np, 'PLISTOPZONE', mempath) + call mem_reallocate(this%idrymeth, np, 'PLIDRYMETH', mempath) call mem_reallocate(this%ifrctrn, np, 'PLIFRCTRN', mempath) call mem_reallocate(this%iexmeth, np, 'PLIEXMETH', mempath) call mem_reallocate(this%extol, np, 'PLEXTOL', mempath) @@ -219,7 +235,7 @@ subroutine load_particle(this, store, imdl, iprp, ip) integer(I4B), intent(in) :: iprp !< index of particle release package particle originated in integer(I4B), intent(in) :: ip !< index into the particle list - call this%transform(reset=.true.) + call this%reset_transform() this%imdl = imdl this%iprp = iprp this%irpt = store%irpt(ip) @@ -227,9 +243,12 @@ subroutine load_particle(this, store, imdl, iprp, ip) this%name = store%name(ip) this%istopweaksink = store%istopweaksink(ip) this%istopzone = store%istopzone(ip) + this%idrymeth = store%idrymeth(ip) + this%icp = store%icp(ip) this%icu = store%icu(ip) this%ilay = store%ilay(ip) this%izone = store%izone(ip) + this%izp = store%izp(ip) this%istatus = store%istatus(ip) this%x = store%x(ip) this%y = store%y(ip) @@ -261,9 +280,12 @@ subroutine save_particle(this, particle, ip) this%name(ip) = particle%name this%istopweaksink(ip) = particle%istopweaksink this%istopzone(ip) = particle%istopzone + this%idrymeth(ip) = particle%idrymeth + this%icp(ip) = particle%icp this%icu(ip) = particle%icu this%ilay(ip) = particle%ilay this%izone(ip) = particle%izone + this%izp(ip) = particle%izp this%istatus(ip) = particle%istatus this%x(ip) = particle%x this%y(ip) = particle%y @@ -285,9 +307,14 @@ subroutine save_particle(this, particle, ip) this%extend(ip) = particle%iextend end subroutine save_particle - !> @brief Apply the given global-to-local transformation to the particle. + !> @brief Transform particle coordinates. + !! + !! Apply a translation and/or rotation to particle coordinates. + !! No rescaling. It's also possible to invert a transformation. + !! Be sure to reset the transformation after using it. + !< subroutine transform_coords(this, xorigin, yorigin, zorigin, & - sinrot, cosrot, invert, reset) + sinrot, cosrot, invert) use GeomUtilModule, only: transform, compose class(ParticleType), intent(inout) :: this !< particle real(DP), intent(in), optional :: xorigin !< x coordinate of origin @@ -296,44 +323,33 @@ subroutine transform_coords(this, xorigin, yorigin, zorigin, & real(DP), intent(in), optional :: sinrot !< sine of rotation angle real(DP), intent(in), optional :: cosrot !< cosine of rotation angle logical(LGP), intent(in), optional :: invert !< whether to invert - logical(LGP), intent(in), optional :: reset !< whether to reset - - ! Reset if requested - if (present(reset)) then - if (reset) then - this%xorigin = DZERO - this%yorigin = DZERO - this%zorigin = DZERO - this%sinrot = DZERO - this%cosrot = DONE - this%cosrot = DONE - this%transformed = .false. - return - end if - end if - ! Otherwise, transform coordinates call transform(this%x, this%y, this%z, & this%x, this%y, this%z, & xorigin, yorigin, zorigin, & sinrot, cosrot, invert) - ! Modify transformation from model coordinates to particle's new - ! local coordinates by incorporating this latest transformation call compose(this%xorigin, this%yorigin, this%zorigin, & this%sinrot, this%cosrot, & xorigin, yorigin, zorigin, & sinrot, cosrot, invert) - ! Set isTransformed flag to true. Note that there is no check - ! to see whether the modification brings the coordinates back - ! to model coordinates (in which case the origin would be very - ! close to zero and sinrot and cosrot would be very close to 0. - ! and 1., respectively, allowing for roundoff error). this%transformed = .true. end subroutine transform_coords - !> @brief Return the particle's model (global) coordinates. + subroutine reset_transform(this) + class(ParticleType), intent(inout) :: this !< particle + + this%xorigin = DZERO + this%yorigin = DZERO + this%zorigin = DZERO + this%sinrot = DZERO + this%cosrot = DONE + this%cosrot = DONE + this%transformed = .false. + end subroutine reset_transform + + !> @brief Return the particle's model coordinates. subroutine get_model_coords(this, x, y, z) use GeomUtilModule, only: transform, compose class(ParticleType), intent(inout) :: this !< particle @@ -342,12 +358,11 @@ subroutine get_model_coords(this, x, y, z) real(DP), intent(out) :: z !< z coordinate if (this%transformed) then - ! Transform back from local to model coordinates + ! Untransform coordinates call transform(this%x, this%y, this%z, x, y, z, & this%xorigin, this%yorigin, this%zorigin, & - this%sinrot, this%cosrot, .true.) + this%sinrot, this%cosrot, invert=.true.) else - ! Already in model coordinates x = this%x y = this%y z = this%z diff --git a/src/Solution/ParticleTracker/vertical.md b/src/Solution/ParticleTracker/vertical.md new file mode 100644 index 00000000000..f0afbc109ba --- /dev/null +++ b/src/Solution/ParticleTracker/vertical.md @@ -0,0 +1,133 @@ +# Vertical tracking + +This document describes the approach PRT takes to vertical particle motion. + +## Legend + +Diagrams use the following conventions. + +* Stadium-shaped boxes represent steps or processes. +* Square boxes represent outcomes. +* Diamond boxes represent conditions (i.e. runtime state). +* Round-corner boxes represent user options. +* Thin lines represent decisions made by the program on the basis of runtime state, e.g. particle, cell, flows. +* Thick lines represent decisions made by the user by way of options. +* Green outcome boxes indicate the particle remains active. +* Red outcome boxes indicate the particle terminates. + +```mermaid +flowchart LR + OPTION[Outcome] + OPTION(OPTION) ==> |Yes| STEP([Step]) + STEP --> ACTIVE + OPTION ==> |No| CONDITION{Condition} + CONDITION --> |Yes| TERMINATE + CONDITION --> |No| ACTIVE + ACTIVE[Active]:::active + TERMINATE[Termination]:::terminate + + classDef active stroke:#98fb98 + classDef terminate stroke:#f08080 +``` + +## The problem + +When a particle is in the flow field, vertical motion can be solved in the same way as lateral motion. Special handling is necessary above the water table. + +A "dry" cell is either 1) an inactive cell or 2) an active-but-dry cell, as can occur with the Newton formulation. + +Normally, an inactive cell might be dry or explicitly disabled (idomain). With Newton, dry cells remain active. + +## The approach + +Release-time and tracking-time considerations are described (and implemented) separately. + +Each particle is either released into the simulation, or terminates unreleased. In the former case the particle's first record will be reason 0 (release), status 1 (active). In the latter reason 3 (termination), status 8 (permanently unreleased). + +At each time step, the PRT model applies the tracking method to each particle. The particle's trajectory is solved over the model grid until the end of the time step, or until the particle terminates (due e.g. to stop time or encountering a termination condition), whichever occurs first. + +Particles may traverse an arbitrary number of cells in a time step, in the lateral as well as vertical dimensions. + +Sometimes it is convenient to avoid "stranding" particles — rather than terminating dry particles, it is often convenient instead to move them down to the saturated zone and continue tracking. PRT allows particles (and indeed configures them by default) to move instantaneously down to the water table in dry conditions. + +### Release + +At release time, PRT decides whether to release each particle or to terminate it unreleased. + +If the release cell is active, the particle will be released at the specified coordinates. + +If the release cell is inactive, behavior is determined by the `DRAPE` option. If the `DRAPE` option is enabled, the particle will be "draped" down to and released from the top-most active cell beneath it, if any. If there is no active cell underneath the particle in any layer, or if `DRAPE` is not enabled, the particle will terminate unreleased (with status code 8). + +Since under the Newton formulation dry cells can remain active, the `DRAPE` option has no effect when Newton is on (assuming particles are not released into disabled grid regions). Vertical tracking behavior with Newton can be configured with tracking-time settings. + +```mermaid +flowchart LR + ACTIVE --> |No| DRAPE(DRAPE) + ACTIVE{Cell active?} --> |Yes| RELEASE[Release in specified cell]:::release + DRAPE ==> |No| TERMINATE:::terminate + DRAPE ==> |Yes| ACTIVE_UNDER{Active under?} + ACTIVE_UNDER --> |Yes| RELEASE_AT_TABLE[Drape to active cell]:::release + ACTIVE_UNDER --> |No| TERMINATE[Terminate] + + classDef release stroke:#98fb98 + classDef terminate stroke:#f08080 +``` + +### Tracking + +A particle might find itself above the water table for one of two reasons: + +1. It was released above the water table. + + With the Newton formulation, particles can be released into dry-but-active cells. + +2. The water table has receded. + + Particle trajectories are solved over the same time discretization used by the flow model. A particle may be immersed in the flow field in one time step, and find that the water table has dropped below it in the next. + +Tracking and termination decisions are made on the basis of information like + +1) a cell's active status +2) whether the cell is dry +3) whether the cell has outgoing flow across any face +4) whether the particle is dry (above the water table) +5) the particle's prior path + +A particle which finds itself in an inactive cell will terminate with status code 7. This is consistent with MODPATH 7's behavior. + +A particle in a dry-but-active cell, or above the water table in a partially saturated cell, need not terminate. We call such a particle dry. MODFLOW version 6.6.0 introduces a new option `DRY_TRACKING_METHOD` for the PRP package, determining how dry particles should behave. Supported values are: + +- `DROP` (default) +- `STOP` +- `STAY` + +If `DROP` is selected, or if a `DRY_TRACKING_METHOD` is unspecified, a particle in a dry position is passed vertically and instantaneously to the water table (if the cell is partially saturated) or to the bottom of the cell (if the cell is dry). This repeats (i.e. the particle may drop through multiple cells) until it reaches the water table. Tracking then proceeds as usual. + +**Note**: A divide-by-zero crash has been fixed for `gfortran`, which could occur upon a particle's entry into a dry cell in a structured grid. + +If `STOP` is selected, dry particles will be terminated. + +If `STAY` is selected, a dry particle will remain stationary until a) the water table rises and tracking can continue or b) the simulation ends. + +```mermaid +flowchart LR + ACTIVE{Cell active?} --> |No| TERMINATE{Terminate} + ACTIVE{Cell active?} --> |Yes| PARTICLE_DRY + PARTICLE_DRY{Particle dry?} --> |Yes| DRY_TRACKING_METHOD(DRY_TRACKING_METHOD) + DRY_TRACKING_METHOD ==> |STOP| TERMINATE[Terminate]:::terminate + DRY_TRACKING_METHOD ==> |DROP| CELL_DRY{Cell dry?} + CELL_DRY --> |Yes| DROP_BOTTOM[Pass to cell bottom]:::track + CELL_DRY --> |No| DROP_TABLE([Pass to water table]) + DRY_TRACKING_METHOD ==> |STAY| STAY[Stationary]:::track + DROP_TABLE --> TRACK:::track + PARTICLE_DRY --> |No| TRACK[Track normally] + + classDef track stroke:#98fb98 + classDef terminate stroke:#f08080 +``` + +#### Caveat + +In MF6.5, behavior was as described by `DROP`, with one major exception: lack of an exit face (i.e. any face with outgoing flow) took precedence over cell saturation; a particle finding itself in a dry cell with no outgoing flow would previously terminate, where if `DROP` is selected (or a dry tracking method unspecified) the pass-to-bottom method will now be applied instead. + +With this change, it also becomes necessary to prohibit backtracking between vertically adjacent pairs of cells within the same time step, in order to avoid the possibility of infinite loops — a particle might otherwise be passed endlessly between e.g. the bottom face of a cell containing a pumping well and the top face of the cell below. Note that this limitation applies only to vertically adjacent cells, and only to the immediately previous cell — a particle may re-enter a cell after entering a third cell. \ No newline at end of file