Coverage for C:\src\imod-python\imod\wq\slv.py: 34%

92 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2024-04-08 10:26 +0200

1import pathlib 

2 

3import jinja2 

4import pandas as pd 

5 

6from imod.wq.pkgbase import Package 

7 

8 

9class PreconditionedConjugateGradientSolver(Package): 

10 """ 

11 The Preconditioned Conjugate Gradient Solver is used to solve the finite 

12 difference equations in each step of a MODFLOW stress period. 

13 

14 Parameters 

15 ---------- 

16 max_iter: int 

17 is the maximum number of outer iterations - that is, calss to the 

18 solutions routine (MXITER). For a linear problem max_iter should be 1, unless 

19 more than 50 inner iterations are required, when max_iter could be as 

20 large as 10. A larger number (generally less than 100) is required for a 

21 nonlinear problem. 

22 inner_iter: int 

23 is the number of inner iterations (iter1). For nonlinear problems, 

24 inner_iter usually ranges from 10 to 30; a value of 30 will be 

25 sufficient for most linear problems. 

26 rclose: float 

27 is the residual criterion for convergence, in units of cubic length per 

28 time. The units for length and time are the same as established for all 

29 model data. When the maximum absolute value of the residual at all nodes 

30 during an iteration is less than or equal to RCLOSE, and the criterion 

31 for HCLOSE is also satisfied (see below), iteration stops. 

32 

33 Default value: 100.0. **Nota bene**: this is aimed at regional modelling. 

34 For detailed studies (e.g. lab experiments) much smaller values can be 

35 required. 

36 Very general rule of thumb: should be less than 10% of smallest cell volume. 

37 hclose: float 

38 is the head change criterion for convergence, in units of length. When 

39 the maximum absolute value of head change from all nodes during an 

40 iteration is less than or equal to HCLOSE, and the criterion for RCLOSE 

41 is also satisfied (see above), iteration stops. 

42 Default value: 1.0e-4. **Nota bene**: This is aimed at regional modelling, ` 

43 for detailed studies (e.g. lab experiments) much smaller values can be 

44 required. 

45 relax: float, optional 

46 is the relaxation parameter used. Usually, RELAX = 1.0, but for some 

47 problems a value of 0.99, 0.98, or 0.97 will reduce the number of 

48 iterations required for convergence. 

49 Default value: 0.98. 

50 damp: float, optional 

51 is the damping factor. It is typically set equal to one, which indicates 

52 no damping. A value less than 1 and greater than 0 causes damping. DAMP 

53 does not affect inner iterations; instead, it affects outer iterations. 

54 Default value: 1.0. 

55 """ 

56 

57 _pkg_id = "pcg" 

58 _template = ( 

59 "[pcg]\n" 

60 " mxiter = {max_iter}\n" 

61 " iter1 = {inner_iter}\n" 

62 " npcond = 1\n" 

63 " hclose = {hclose}\n" 

64 " rclose = {rclose}\n" 

65 " relax = {relax}\n" 

66 " iprpcg = 1\n" 

67 " mutpcg = 0\n" 

68 " damp = {damp}" 

69 ) 

70 

71 def __init__( 

72 self, 

73 max_iter=150, 

74 inner_iter=100, 

75 rclose=1000.0, 

76 hclose=1.0e-4, 

77 relax=0.98, 

78 damp=1.0, 

79 ): 

80 super().__init__() 

81 self["max_iter"] = max_iter 

82 self["inner_iter"] = inner_iter 

83 self["rclose"] = rclose 

84 self["hclose"] = hclose 

85 self["relax"] = relax 

86 self["damp"] = damp 

87 

88 def _pkgcheck(self, ibound=None): 

89 to_check = ["max_iter", "inner_iter", "rclose", "hclose", "damp"] 

90 self._check_positive(to_check) 

91 

92 

93class GeneralizedConjugateGradientSolver(Package): 

94 """ 

95 The Generalized Conjugate Gradient Solver solves the matrix equations 

96 resulting from the implicit solution of the transport equation. 

97 

98 Parameters 

99 ---------- 

100 max_iter: int 

101 is the maximum number of outer iterations (MXITER); it should be set to an 

102 integer greater than one (1) only when a nonlinear sorption isotherm is 

103 included in simulation. 

104 iter1: int 

105 is the maximum number of inner iterations (iter1); a value of 30-50 

106 should be adequate for most problems. 

107 isolve: int 

108 is the type of preconditioners to be used with the Lanczos/ORTHOMIN 

109 acceleration scheme: 

110 isolve = 1: Jacobi 

111 isolve = 2: SSOR 

112 isolve = 3: Modified Incomplete Cholesky (MIC) 

113 (MIC usually converges faster, but it needs significantly more memory) 

114 lump_dispersion: bool 

115 is an integer flag for treatment of dispersion tensor cross terms: 

116 ncrs = 0: lump all dispersion cross terms to the right-hand-side 

117 (approximate but highly efficient). 

118 ncrs = 1: include full dispersion tensor (memory intensive). 

119 cclose: float 

120 is the convergence criterion in terms of relative concentration; a real 

121 value between 10-4 and 10-6 is generally adequate. 

122 """ 

123 

124 _pkg_id = "gcg" 

125 _template = ( 

126 "[gcg]\n" 

127 " mxiter = {max_iter}\n" 

128 " iter1 = {inner_iter}\n" 

129 " isolve = {preconditioner}\n" 

130 " ncrs = {lump_dispersion}\n" 

131 " cclose = {cclose}\n" 

132 " iprgcg = 0" 

133 ) 

134 

135 _keywords = { 

136 "preconditioner": {"jacobi": 1, "ssor": 2, "mic": 3}, 

137 "lump_dispersion": {True: 0, False: 1}, 

138 } 

139 

140 def __init__( 

141 self, 

142 max_iter=1, 

143 inner_iter=50, 

144 cclose=1.0e-6, 

145 preconditioner="mic", 

146 lump_dispersion=True, 

147 ): 

148 super().__init__() 

149 self["max_iter"] = max_iter 

150 self["inner_iter"] = inner_iter 

151 self["cclose"] = cclose 

152 self["preconditioner"] = preconditioner 

153 self["lump_dispersion"] = lump_dispersion 

154 

155 def _pkgcheck(self, ibound=None): 

156 to_check = ["max_iter", "inner_iter", "cclose"] 

157 self._check_positive(to_check) 

158 

159 

160class ParallelSolver(Package): 

161 """ 

162 Base package for the parallel solvers. 

163 """ 

164 

165 def _compute_load_balance_weight(self, ibound): 

166 if self["partition"] == "rcb": 

167 if self["load_balance_weight"].values[()] is None: 

168 self["load_balance_weight"] = (ibound != 0).sum("layer").astype(float) 

169 

170 def _render(self, directory): 

171 d = {k: v.values for k, v in self.dataset.data_vars.items()} 

172 if hasattr(self, "_keywords"): 

173 for key in self._keywords.keys(): 

174 self._replace_keyword(d, key) 

175 

176 if self["partition"] == "rcb": 

177 d["load_balance_weight"] = self._compose_path( 

178 { 

179 "directory": directory, 

180 "name": "load_balance_weight", 

181 "extension": ".asc", 

182 } 

183 ) 

184 

185 return self._template.render(**d) 

186 

187 def save(self, directory): 

188 """ 

189 Overloaded method to write .asc instead of .idf. 

190 (This is an idiosyncracy of the parallel iMODwq code.) 

191 """ 

192 # TODO: remove this when iMOD_wq accepts idfs for the load grid. 

193 da = self["load_balance_weight"] 

194 if not da.dims == ("y", "x"): 

195 raise ValueError( 

196 "load_balance_weight dimensions must be exactly ('y', 'x')." 

197 ) 

198 path = pathlib.Path(directory).joinpath("load_balance_weight.asc") 

199 path.parent.mkdir(exist_ok=True, parents=True) 

200 pd.DataFrame(da.values).to_csv( 

201 path, sep="\t", header=False, index=False, float_format="%8.2f" 

202 ) 

203 

204 

205class ParallelKrylovFlowSolver(ParallelSolver): 

206 """ 

207 The Parallel Krylov Flow Solver is used for parallel solving of the flow 

208 model. 

209 

210 Parameters 

211 ---------- 

212 max_iter: int 

213 is the maximum number of outer iterations (MXITER); it should be set to 

214 an integer greater than one (1) only when a nonlinear sorption isotherm 

215 is included in simulation. 

216 inner_iter: int 

217 is the maximum number of inner iterations (INNERIT); a value of 30-50 

218 should be adequate for most problems. 

219 hclose: float 

220 is the head change criterion for convergence (HCLOSEPKS), in units of 

221 length. When the maximum absolute value of head change from all nodes 

222 during an iteration is less than or equal to HCLOSE, and the criterion 

223 for RCLOSE is also satisfied (see below), iteration stops. 

224 rclose: float 

225 is the residual criterion for convergence (RCLOSEPKS), in units of cubic 

226 length per time. The units for length and time are the same as 

227 established for all model data. When the maximum absolute value of the 

228 residual at all nodes during an iteration is less than or equal to 

229 RCLOSE, and the criterion for HCLOSE is also satisfied (see above), 

230 iteration stops. 

231 relax: float 

232 is the relaxation parameter used. Usually, RELAX = 1.0, but for some 

233 problems a value of 0.99, 0.98, or 0.97 will reduce the number of 

234 iterations required for convergence. 

235 h_fstrict: float, optional 

236 is a factor to apply to HCLOSE to set a stricter hclose for the linear 

237 inner iterations (H_FSTRICTPKS). HCLOSE_inner is calculated as follows: 

238 HCLOSEPKS * H_FSTRICTPKS. 

239 r_fstrict: float, optional 

240 is a factor to apply to RCLOSE to set a stricter rclose for the linear 

241 inner iterations (R_FSTRICTPKS). RCLOSE_inner is calculated as follows: 

242 RCLOSEPKS * R_FSTRICTPKS. 

243 partition: {"uniform", "rcb"}, optional 

244 Partitioning option (PARTOPT). "uniform" partitions the model domain 

245 into equally sized subdomains. "rcb" (Recursive Coordinate Bisection) 

246 uses a 2D pointer grid with weights to partition the model domain. 

247 Default value: "uniform" 

248 solver: {"pcg"}, optional 

249 Flag indicating the linear solver to be used (ISOLVER). 

250 Default value: "pcg" 

251 preconditioner: {"ilu"}, optional 

252 Flag inicating the preconditioner to be used (NPC). 

253 Devault value: "ilu" 

254 deflate: {True, False}, optional 

255 Flag for deflation preconditioner. 

256 Default value: False 

257 debug: {True, False}, optional 

258 Debug option. 

259 Default value: False 

260 load_balance_weight: xarray.DataArray, optional 

261 2D grid with load balance weights, used when partition = "rcb" 

262 (Recursive Coordinate Bisection). If None (default), then the module 

263 will create a load balance grid by summing active cells over layers: 

264 `(ibound != 0).sum("layer")` 

265 

266 Note that even though the iMOD-SEAWAT helpfile states .idf is 

267 accepted, it is not. This load balance grid should be a .asc file 

268 (without a header). Formatting is done as follows: 

269 `pd.DataFrame(load_balance_weight.values).to_csv(path, sep='\\t', 

270 header=False, index=False, float_format = "%8.2f")` 

271 """ 

272 

273 _pkg_id = "pksf" 

274 _template = jinja2.Template( 

275 "[pksf]\n" 

276 " mxiter = {{max_iter}}\n" 

277 " innerit = {{inner_iter}}\n" 

278 " hclosepks = {{hclose}}\n" 

279 " rclosepks = {{rclose}}\n" 

280 " relax = {{relax}}\n" 

281 " h_fstrictpks = {{h_fstrict}}\n" 

282 " r_fstrictpks = {{r_fstrict}}\n" 

283 " partopt = {{partition}}\n" 

284 " isolver = {{solver}}\n" 

285 " npc = {{preconditioner}}\n" 

286 " npcdef = {{deflate}}\n" 

287 "{% if load_balance_weight %} loadptr = {{load_balance_weight}}\n{% endif %}" 

288 " pressakey = {{debug}}\n" 

289 ) 

290 _keywords = { 

291 "partition": {"uniform": 0, "rcb": 5}, 

292 "solver": {"pcg": 1}, 

293 "preconditioner": {"ilu": 2}, 

294 "deflate": {False: 0, True: 1}, 

295 } 

296 

297 def __init__( 

298 self, 

299 max_iter=150, 

300 inner_iter=100, 

301 hclose=1.0e-4, 

302 rclose=1000.0, 

303 relax=0.98, 

304 h_fstrict=1.0, 

305 r_fstrict=1.0, 

306 partition="uniform", 

307 solver="pcg", 

308 preconditioner="ilu", 

309 deflate=False, 

310 debug=False, 

311 load_balance_weight=None, 

312 ): 

313 super().__init__() 

314 self["max_iter"] = max_iter 

315 self["inner_iter"] = inner_iter 

316 self["hclose"] = hclose 

317 self["rclose"] = rclose 

318 self["relax"] = relax 

319 self["h_fstrict"] = h_fstrict 

320 self["r_fstrict"] = r_fstrict 

321 self["partition"] = partition 

322 self["solver"] = solver 

323 self["preconditioner"] = preconditioner 

324 self["deflate"] = deflate 

325 self["debug"] = debug 

326 self["load_balance_weight"] = load_balance_weight 

327 

328 def _pkgcheck(self, ibound=None): 

329 to_check = [ 

330 "hclose", 

331 "rclose", 

332 "h_fstrict", 

333 "r_fstrict", 

334 "max_iter", 

335 "inner_iter", 

336 "relax", 

337 ] 

338 self._check_positive(to_check) 

339 # TODO: fix 

340 # Check whether option is actually an available option 

341 # for opt_arg in self._keywords.keys(): 

342 # if self[opt_arg] not in self._keywords[opt_arg].keys(): 

343 # raise ValueError( 

344 # "Argument for {} not in {}, instead got {}".format( 

345 # opt_arg, self._keywords[opt_arg].keys(), self[opt_arg] 

346 # ) 

347 # ) 

348 # 

349 

350 

351class ParallelKrylovTransportSolver(ParallelSolver): 

352 """ 

353 The Parallel Krylov Transport Solver is used for parallel solving of the 

354 transport model. 

355 

356 Parameters 

357 ---------- 

358 max_iter: int 

359 is the maximum number of outer iterations (MXITER); it should be set to 

360 an integer greater than one (1) only when a nonlinear sorption isotherm 

361 is included in simulation. 

362 inner_iter: int 

363 is the maximum number of inner iterations (INNERIT); a value of 30-50 

364 should be adequate for most problems. 

365 cclose: float, optional 

366 is the convergence criterion in terms of relative concentration; a real 

367 value between 10-4 and 10-6 is generally adequate. 

368 Default value: 1.0e-6. 

369 relax: float, optional 

370 is the relaxation parameter used. Usually, RELAX = 1.0, but for some 

371 problems a value of 0.99, 0.98, or 0.97 will reduce the number of 

372 iterations required for convergence. 

373 Default value: 0.98. 

374 partition: {"uniform", "rcb"}, optional 

375 Partitioning option (PARTOPT). "uniform" partitions the model domain 

376 into equally sized subdomains. "rcb" (Recursive Coordinate Bisection) 

377 uses a 2D pointer grid with weights to partition the model domain. 

378 Default value: "uniform". 

379 solver: {"bicgstab", "gmres", "gcr"}, optional 

380 Flag indicating the linear solver to be used (ISOLVER). 

381 Default value: "bicgstab" 

382 preconditioner: {"ilu"}, optional 

383 Flag inicating the preconditioner to be used (NPC). 

384 Devault value: "ilu". 

385 debug: {True, False}, optional 

386 Debug option. 

387 Default value: False 

388 load_balance_weight: xarray.DataArray, optional 

389 2D grid with load balance weights, used when partition = "rcb" 

390 (Recursive Coordinate Bisection). If None (default), then the module 

391 will create a load balance grid by summing active cells over layers: 

392 `(ibound != 0).sum("layer")` 

393 

394 Note that even though the iMOD-SEAWAT helpfile states .idf is 

395 accepted, it is not. This load balance grid should be a .asc file 

396 (without a header). Formatting is done as follows: 

397 `pd.DataFrame(load_balance_weight.values).to_csv(path, sep='\\t', 

398 header=False, index=False, float_format = "%8.2f")` 

399 """ 

400 

401 _pkg_id = "pkst" 

402 _template = jinja2.Template( 

403 "[pkst]\n" 

404 " mxiter = {{max_iter}}\n" 

405 " innerit = {{inner_iter}}\n" 

406 " cclosepks = {{cclose}}\n" 

407 " relax = {{relax}}\n" 

408 " partopt = {{partition}}\n" 

409 " isolver = {{solver}}\n" 

410 " npc = {{preconditioner}}\n" 

411 "{% if load_balance_weight %} loadptr = {{load_balance_weight}}\n{% endif %}" 

412 " pressakey = {{debug}}\n" 

413 ) 

414 

415 _keywords = { 

416 "partition": {"uniform": 0, "rcb": 5}, 

417 "solver": {"bicgstab": 2, "gmres": 3, "gcr": 4}, 

418 "preconditioner": {"ilu": 2}, 

419 } 

420 

421 def __init__( 

422 self, 

423 max_iter=1, 

424 inner_iter=50, 

425 cclose=1.0e-6, 

426 relax=0.98, 

427 partition="uniform", 

428 solver="bicgstab", 

429 preconditioner="ilu", 

430 debug=False, 

431 load_balance_weight=None, 

432 ): 

433 super().__init__() 

434 self["max_iter"] = max_iter 

435 self["inner_iter"] = inner_iter 

436 self["cclose"] = cclose 

437 self["relax"] = relax 

438 self["partition"] = partition 

439 self["solver"] = solver 

440 self["preconditioner"] = preconditioner 

441 self["debug"] = debug 

442 self["load_balance_weight"] = load_balance_weight 

443 

444 def _pkgcheck(self, ibound=None): 

445 to_check = ["cclose", "max_iter", "inner_iter", "relax"] 

446 self._check_positive(to_check) 

447 # TODO: fix 

448 # Check whether option is actually an available option 

449 # for opt_arg in self._keywords.keys(): 

450 # if self[opt_arg] not in self._keywords[opt_arg].keys(): 

451 # raise ValueError( 

452 # "Argument for {} not in {}, instead got {}".format( 

453 # opt_arg, self._keywords[opt_arg].keys(), self[opt_arg] 

454 # ) 

455 # )