• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

p-ortega / mf6rtm / 21449732138

28 Jan 2026 04:20AM UTC coverage: 84.527% (+0.01%) from 84.513%
21449732138

push

github

web-flow
Merge pull request #51 from p-ortega/develop

Patch mup3d into main

30 of 39 new or added lines in 5 files covered. (76.92%)

2 existing lines in 2 files now uncovered.

1464 of 1732 relevant lines covered (84.53%)

0.85 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

88.82
/mf6rtm/simulation/solver.py
1
"""The solver module provides the Mf6RTM class that couples modflowapi and
2
phreeqcrm, along with functions to run the coupled simulations.
3
"""
4
import os
1✔
5
# import warnings
6
import numpy as np
1✔
7

8
from datetime import datetime
1✔
9
from typing import Any, Union
1✔
10
from pathlib import Path
1✔
11

12
from PIL import Image
1✔
13
from mf6rtm.simulation.mf6api import Mf6API
1✔
14
from mf6rtm.simulation.phreeqcbmi import PhreeqcBMI
1✔
15
from mf6rtm.simulation.discretization import total_cells_in_grid
1✔
16
from mf6rtm.config.config import MF6RTMConfig
1✔
17
from mf6rtm.io.externalio import SelectedOutput
1✔
18
from mf6rtm.utils import utils
1✔
19

20
# warnings.filterwarnings("ignore")
21
# warnings.filterwarnings("ignore", category=DeprecationWarning)
22

23
# global variables
24
DT_FMT = "%Y-%m-%d %H:%M:%S"
1✔
25

26
time_units_dict = {
1✔
27
    "seconds": 1,
28
    "minutes": 60,
29
    "hours": 3600,
30
    "days": 86400,
31
    "years": 31536000,
32
    "unknown": 1,  # if unknown assume seconds
33
}
34

35
def check_config_file(wd: os.PathLike) -> tuple[os.PathLike, os.PathLike]:
1✔
36
    assert os.path.exists(
1✔
37
        os.path.join(wd, "mf6rtm.toml")
38
        ), "mf6rtm.toml not found in model directory"
39
    config_file= os.path.join(wd, "mf6rtm.toml")
1✔
40
    config = MF6RTMConfig.from_toml_file(config_file)
1✔
41

42
    # validate config values like timing and tsteps
43
    config._validate_config()
1✔
44
    return config
1✔
45

46
def check_nam_files(wd:os.PathLike) -> tuple[os.PathLike,os.PathLike]:
1✔
47
    """Check if the nam files are present in the model directory"""
48
    nam = [f for f in os.listdir(wd) if f.endswith(".nam")]
1✔
49
    assert "mfsim.nam" in nam, "mfsim.nam file not found in model directory"
1✔
50
    # assert "gwf.nam" in nam, "gwf.nam file not found in model directory"
51
    return os.path.join(wd, "mfsim.nam") #, os.path.join(wd, "gwf.nam")
1✔
52

53
def prep_to_run(wd:os.PathLike) -> tuple[os.PathLike,os.PathLike]:
1✔
54
    """
55
    Prepares the model to run by checking if the model directory (wd) contains the necessary files
56
    and returns the path to the yaml file (phreeqcrm) and the dll file (mf6 api)
57

58
    Parameters
59
    ----------
60
    wd :os.PathLike
61
        The path to the working directory of model directory
62
    Returns
63
    -------
64
    tuple[PathLike,os.PathLike]
65
        The path to the phreeqcrm model file (yaml) and the path to the MODFLOW 6 dll (associated with mf6api).
66
    """
67
    # check if wd exists
68
    assert os.path.exists(wd), f"Path {wd} not found"
1✔
69
    # check if file starting with libmf6 exists
70
    dll = [f for f in os.listdir(wd) if f.startswith("libmf6")]
1✔
71
    assert len(dll) == 1, "libmf6 dll not found in model directory"
1✔
72
    dll = os.path.join(wd, "libmf6")
1✔
73

74
    config = check_config_file(wd)
1✔
75
    check_nam_files(wd)
1✔
76
    if config.reactive['externalio']:
1✔
77
        from mf6rtm.io.externalio import Regenerator
1✔
78
        print("WARNING: Flag for external IO mode is active")
1✔
79
        regcls = Regenerator.regenerate_from_external_files(wd=wd,
1✔
80
                                                phinpfile='phinp.dat',
81
                                                yamlfile='mf6rtm.yaml',
82
                                                dllfile=dll
83
                                                )
84
        yamlfile = regcls.yamlfile
1✔
85
    else:
86
        yamlfile = os.path.join(wd, "mf6rtm.yaml")
1✔
87
    assert os.path.exists(
1✔
88
        os.path.join(wd, yamlfile)
89
    ), f"{yamlfile} not found in model directory {wd}"
90
    return yamlfile, dll
1✔
91

92
def solve(wd:os.PathLike, reactive: Union[bool, None] = None, nthread: int = 1) -> bool:
1✔
93
    """Wrapper to prepare and call solve functions"""
94

95
    mf6rtm = initialize_interfaces(wd, nthread=nthread)
1✔
96
    if reactive is not None and isinstance(reactive, bool) and reactive != mf6rtm.reactive:
1✔
97
        print(
1✔
98
                f"Mode changed from "
99
                f"{'reactive' if mf6rtm.reactive else 'non-reactive'} to "
100
                f"{'reactive' if reactive else 'non-reactive'}\n"
101
            )
102
        mf6rtm._set_reactive(reactive)
1✔
103
        # let mf6 manage this for conservative runs
104
        mf6rtm.selected_output.get_selected_output_on = False
1✔
105
    mf6rtm.print_warning_user_active()
1✔
106
    success = mf6rtm._solve()
1✔
107
    return success
1✔
108

109

110
# TODO: we should maybe move this into the Mf6API as an alternative constructor
111
def initialize_interfaces(wd:os.PathLike, nthread: int = 1) -> Mf6API:
1✔
112
    """Function to initialize the interfaces for modflowapi and phreeqcrm and returns the mf6rtm object"""
113

114
    yamlfile, dll = prep_to_run(wd)
1✔
115

116
    if nthread > 1:
1✔
117
        # set nthreds to nthread
118
        set_nthread_yaml(yamlfile, nthread=nthread)
1✔
119

120
    # initialize the interfaces
121
    mf6api = Mf6API(wd, dll)
1✔
122
    phreeqcrm = PhreeqcBMI(yamlfile) #FIXME: Does not work with path like
1✔
123
    mf6rtm = Mf6RTM(wd, mf6api, phreeqcrm)
1✔
124
    return mf6rtm
1✔
125

126

127
def set_nthread_yaml(yamlfile:os.PathLike, nthread: int = 1) -> None:
1✔
128
    """Function to set the number of threads in the yaml file"""
129
    with open(yamlfile, "r") as f:
1✔
130
        lines = f.readlines()
1✔
131
    for i, line in enumerate(lines):
1✔
132
        if "nthreads" in line:
1✔
133
            lines[i] = f"  nthreads: {nthread}\n"
1✔
134
    with open(yamlfile, "w") as f:
1✔
135
        f.writelines(lines)
1✔
136
    return
1✔
137

138

139
class Mf6RTM(object):
1✔
140
    def __init__(
1✔
141
        self,
142
        wd:os.PathLike,
143
        mf6api: Mf6API,
144
        phreeqcbmi: PhreeqcBMI,
145
    ) -> None:
146
        """
147
        Initialize the Mf6RTM instance with specified working directory, MF6API,
148
        and PhreeqcBMI instances.
149

150
        Parameters
151
        ----------
152
        wd :os.PathLike
153
            The working directory path for the model.
154
        mf6api : Mf6API
155
            An instance of the Mf6API class, representing the Modflow 6 API.
156
        phreeqcbmi : PhreeqcBMI
157
            An instance of the PhreeqcBMI class, representing the PHREEQC BMI.
158

159
        Attributes
160
        ----------
161
        mf6api : Mf6API
162
            The Modflow 6 API instance.
163
        phreeqcbmi : PhreeqcBMI
164
            The PHREEQC BMI instance.
165
        charge_offset : float
166
            Offset for charge, initialized to 0.0.
167
        wd :os.PathLike
168
            The working directory path.
169
        sout_fname : str
170
            Filename for the output, default is "sout.csv".
171
        reactive : bool
172
            Flag indicating if the model is reactive, default is True.
173
        epsaqu : float
174
            ??Epsaqueous value??, initialized to 0.0.
175
        fixed_components : Any
176
            Fixed components, default is None.
177
        get_selected_output_on : bool
178
            Flag indicating if selected output is on, default is True.
179
        component_model_dict : dict[str, str]
180
            Dictionary mapping PHREEQC aqueous chemical components to their
181
            corresponding Modflow 6 groundwater transport (gwt6) model names.
182
        conservative_transport_models: list[str]
183
            List of Modflow 6 groundwater transport (gwt6) model not coupled
184
            with PhreeqcRM
185
        nxyz : int
186
            Total number of cells in the grid.
187
        """
188
        assert isinstance(mf6api, Mf6API), "MF6API must be an instance of Mf6API"
1✔
189
        assert isinstance(
1✔
190
            phreeqcbmi, PhreeqcBMI
191
        ), "PhreeqcBMI must be an instance of PhreeqcBMI"
192
        self.mf6api = mf6api
1✔
193
        self.phreeqcbmi = phreeqcbmi
1✔
194
        self.charge_offset = 0.0
1✔
195
        self.wd = Path(wd)
1✔
196
        self.epsaqu = 0.0
1✔
197
        self.fixed_components = None
1✔
198
        self.selected_output = SelectedOutput(self)
1✔
199

200
        # set component model dictionary & list of conservative_transport_models
201
        self.component_model_dict, self.conservative_transport_models = self._create_component_model_dict()
1✔
202

203
        # set discretization
204
        self.nxyz = total_cells_in_grid(self.mf6api)
1✔
205
        # set time conversion factor
206
        self.set_time_conversion()
1✔
207

208
        self.config = MF6RTMConfig.from_toml_file(self.wd/"mf6rtm.toml")
1✔
209
        self.reactive = self.config.reactive['enabled']
1✔
210
        self.set_emulator_training()
1✔
211

212
    def set_emulator_training(self) -> None:
1✔
213
        """
214
        Configure emulator training output.
215

216
        Reads ``emulator_training_data`` from the configuration. If enabled,
217
        sets up emulator output variables; otherwise disables training data.
218

219
        Attributes
220
        ----------
221
        ml_output : bool
222
            Whether emulator training data output is enabled.
223

224
        Returns
225
        -------
226
        None
227
        """
228
        self.ml_output = bool(getattr(self.config, "emulator_training_data", False))
1✔
229

230
        if self.ml_output:
1✔
231
            self.set_emulator_output_add_variables()
×
232
            print("Saving emulator training data for surrogating")
×
233

234

235
    def set_emulator_output_add_variables(self) -> None:
1✔
236
        """
237
        Add emulator target and feature variables to the output.
238

239
        Updates ``selected_output`` with variables defined in the configuration.
240
        Defaults to empty lists if not provided.
241

242
        Attributes
243
        ----------
244
        selected_output.target_var : list of str
245
            Target variables for emulator training.
246
        selected_output.feat_var : list of str
247
            Feature variables for emulator training.
248

249
        Returns
250
        -------
251
        None
252
        """
253
        self.selected_output.target_var = getattr(
×
254
            self.config, "emulator_target_variables", []
255
        )
256
        self.selected_output.feat_var = getattr(
×
257
            self.config, "emulator_feature_variables", []
258
        )
259

260
    def print_warning_user_active(self):
1✔
261
        """
262
        Prints a warning if reaction timing is set to 'user'.
263
        """
264
        if self.config.reactive['timing'] == 'user':
1✔
265
            print(f"WARNING: Running reaction only in the following periods and time steps:")
×
NEW
266
            for period, timestep in self.config.reactive['tsteps']:
×
267
                print(f"  Period {period}, Time step {timestep}")
×
268
        else:
269
            return
1✔
270

271
    def get_saturation_from_mf6(self) -> dict[Any, np.ndarray]:
1✔
272
        """
273
        Get the saturation
274

275
        Parameters
276
        ----------
277
        mf6 (modflowapi): the modflow api object
278

279
        Returns
280
        -------
281
        array: the saturation
282
        """
283
        sat = {
1✔
284
            component: self.mf6api.get_value(
285
                self.mf6api.get_var_address(
286
                    "FMI/GWFSAT",
287
                    f"{self.component_model_dict[component]}"
288
                )
289
            )
290
            for component in self.phreeqcbmi.components
291
        }
292
        # select the first component to get the length of the array
293
        sat = sat[
1✔
294
            self.phreeqcbmi.components[0]
295
        ]  # saturation is the same for all components
296
        self.phreeqcbmi.sat_now = sat  # set phreeqcmbi saturation
1✔
297
        return sat
1✔
298

299
    def get_time_units_from_mf6(self) -> str:
1✔
300
        """Function to get the time units from mf6"""
301
        return self.mf6api.sim.tdis.time_units.get_data()
1✔
302

303
    def set_time_conversion(self) -> None:
1✔
304
        """Function to set the time conversion factor"""
305
        time_units = self.get_time_units_from_mf6()
1✔
306
        self.time_conversion = 1.0 / time_units_dict[time_units]
1✔
307
        self.phreeqcbmi.SetTimeConversion(self.time_conversion)
1✔
308

309
    def _create_component_model_dict(self)-> tuple[dict[str, str], list[str]]:
1✔
310
        """
311
        Create a dictionary of PHREEQC aqueous chemical component names and
312
        their corresponding Modflow 6 Groundwater Transport (GWT) model names.
313

314
        Returns
315
        -------
316
        component_model_dict : dict[str, str]
317
            A dictionary where the keys are the component names and the values are
318
            the corresponding transport model names.
319
        """
320
        components = self.phreeqcbmi.get_value_ptr("Components")
1✔
321
        # convert np.array to list of pure python strings
322
        components = [str(component) for component in components]
1✔
323

324
        gwt_model_names = [
1✔
325
            name for name in self.mf6api.sim.model_names
326
            if (self.mf6api.sim.get_model(name).model_type == 'gwt6')
327
        ]
328
        gwt_name_prefix = longest_common_substring(gwt_model_names)
1✔
329

330
        component_model_dict = dict(zip(components, [None]*len(components)))
1✔
331
        for component in components:
1✔
332
            for model_name in gwt_model_names:
1✔
333
                if model_name.replace(gwt_name_prefix, "").lower() == component.lower():
1✔
334
                    component_model_dict[component] = model_name
1✔
335
            if (component.lower() == 'charge') and (component_model_dict[component] == None):
1✔
336
                for model_name in gwt_model_names:
×
337
                    if model_name.replace(gwt_name_prefix, "").lower() == 'ch':
×
338
                        component_model_dict[component] = model_name
×
339
            assert (component_model_dict[component] != None,
1✔
340
                f"Component {component} is not matched with a transport model"
341
            )
342

343
        conservative_transport_models = list(
1✔
344
            set(gwt_model_names) - set(component_model_dict.values())
345
        )
346

347
        return component_model_dict, conservative_transport_models
1✔
348

349
    # TODO: remove or have raise not implemented error
350
    def _set_fixed_components(self, fixed_components): ...
351

352
    # TODO: make reactive a property
353
    def _set_reactive(self, reactive: bool) -> None:
1✔
354
        """Set the model to run only transport or transport and reactions"""
355
        self.reactive = reactive
1✔
356

357
    def _prepare_to_solve(self) -> None:
1✔
358
        """Prepare the model to solve"""
359
        # check if sout fname exists
360
        if self.selected_output._check_sout_exist():
1✔
361
            # if found remove it
362
            self.selected_output._rm_sout_file()
1✔
363

364
        self.mf6api._prepare_mf6()
1✔
365
        self.phreeqcbmi._prepare_phreeqcrm_bmi()
1✔
366

367
        # get and write sout headers
368
        self.selected_output._write_sout_headers()
1✔
369

370
    def _set_ctime(self) -> float:
1✔
371
        """Set the current time of the simulation from mf6api"""
372
        self.ctime = self.mf6api.get_current_time()
1✔
373
        self.phreeqcbmi._set_ctime(self.ctime)
1✔
374
        return self.ctime
1✔
375

376
    def _set_etime(self) -> float:
1✔
377
        """Set the end time of the simulation from mf6api"""
378
        self.etime = self.mf6api.get_end_time()
1✔
379
        return self.etime
1✔
380

381
    def _set_time_step(self) -> float:
1✔
382
        self.time_step = self.mf6api.get_time_step()
1✔
383
        return self.time_step
1✔
384

385
    def _finalize(self) -> None:
1✔
386
        """Finalize the APIs"""
387
        self._finalize_mf6api()
1✔
388
        self._finalize_phreeqcrm()
1✔
389

390
    def _finalize_mf6api(self) -> None:
1✔
391
        """Finalize the mf6api"""
392
        self.mf6api.finalize()
1✔
393

394
    def _finalize_phreeqcrm(self) -> None:
1✔
395
        """Finalize the phreeqcrm api"""
396
        self.phreeqcbmi.finalize()
1✔
397

398

399
    def _get_cdlbl_vect(self) -> np.ndarray[np.float64]:
1✔
400
        """Get the concentration array from phreeqc bmi reshape to (ncomps, nxyz)"""
401
        c_dbl_vect = self.phreeqcbmi.GetConcentrations()
1✔
402

403
        conc = [
1✔
404
            c_dbl_vect[i : i + self.nxyz] for i in range(0, len(c_dbl_vect), self.nxyz)
405
        ]  # reshape array
406
        # TODO: refactor to use np.reshape(), which is 2x faster
407
        return conc
1✔
408

409
    def _set_conc_at_current_kstep(self, c_dbl_vect: np.ndarray[np.float64]):
1✔
410
        """Saves the current concentration array to the object"""
411
        self.current_iteration_conc = np.reshape(
1✔
412
            c_dbl_vect, (self.phreeqcbmi.ncomps, self.nxyz)
413
        )
414

415
    def _set_conc_at_previous_kstep(self, c_dbl_vect: np.ndarray[np.float64]):
1✔
416
        """Saves the current concentration array to the object"""
417
        self.previous_iteration_conc = np.reshape(
1✔
418
            c_dbl_vect, (self.phreeqcbmi.ncomps, self.nxyz)
419
        )
420

421
    def _transfer_array_to_mf6(self) -> np.ndarray[np.float64]:
1✔
422
        """Transfer the concentration array to mf6"""
423
        c_dbl_vect = self._get_cdlbl_vect()
1✔
424

425
        # check if reactive cells were skipped due to small changes from transport and replace with previous conc
426
        if self._check_previous_conc_exists() and self._check_inactive_cells_exist(
1✔
427
            self.diffmask
428
        ):
429
            c_dbl_vect = self._replace_inactive_cells(c_dbl_vect, self.diffmask)
×
430
        else:
431
            pass
1✔
432

433
        conc_dict = {}
1✔
434
        for i, c in enumerate(self.phreeqcbmi.components):
1✔
435
            conc_dict[c] = c_dbl_vect[i]
1✔
436
            # Set concentrations in mf6
437
            gwt_model_name = self.component_model_dict[c]
1✔
438
            if gwt_model_name.lower() == "charge":
1✔
439
                self.mf6api.set_value(
1✔
440
                    f"{gwt_model_name.upper()}/X",
441
                    utils.concentration_l_to_m3(conc_dict[c]) + self.charge_offset,
442
                )
443
            else:
444
                self.mf6api.set_value(
1✔
445
                    f"{gwt_model_name.upper()}/X",
446
                    utils.concentration_l_to_m3(conc_dict[c]),
447
                )
448
        return c_dbl_vect
1✔
449

450
    def _check_previous_conc_exists(self) -> bool:
1✔
451
        """Function to replace inactive cells in the concentration array"""
452
        # check if self.previous_iteration_conc is a property
453
        return hasattr(self, "previous_iteration_conc")
1✔
454

455
    def _check_inactive_cells_exist(self, diffmask: np.ndarray[np.float64]) -> bool:
1✔
456
        """Function to check if inactive cells exist in the concentration array"""
457
        inact = utils.get_indices(0, diffmask)
1✔
458
        return len(inact) > 0
1✔
459

460
    def _replace_inactive_cells(
1✔
461
        self,
462
        c_dbl_vect: np.ndarray[np.float64],
463
        diffmask: np.ndarray[np.float64],
464
    ) -> np.ndarray[np.float64]:
465
        """Function to replace inactive cells in the concentration array"""
466
        c_dbl_vect = np.reshape(c_dbl_vect, (self.phreeqcbmi.ncomps, self.nxyz))
×
467
        # get inactive cells
468
        inactive_idx = [
×
469
            utils.get_indices(0, diffmask) for k in range(self.phreeqcbmi.ncomps)
470
        ]
471
        c_dbl_vect[:, inactive_idx] = self.previous_iteration_conc[:, inactive_idx]
×
472
        c_dbl_vect = c_dbl_vect.flatten()
×
473
        conc = [
×
474
            c_dbl_vect[i : i + self.nxyz] for i in range(0, len(c_dbl_vect), self.nxyz)
475
        ]
476
        return conc
×
477

478
    def _transfer_array_to_phreeqcrm(self) -> np.ndarray[np.float64]:
1✔
479
        """Transfer the concentration array to phreeqc bmi"""
480
        mf6_conc_array = []
1✔
481
        for c in self.phreeqcbmi.components:
1✔
482
            if c.lower() == "charge":
1✔
483
                mf6_conc_array.append(
1✔
484
                    utils.concentration_m3_to_l(
485
                        self.mf6api.get_value(
486
                            self.mf6api.get_var_address(
487
                                "X",
488
                                f"{self.component_model_dict[c].upper()}",
489
                            )
490
                        )
491
                        - self.charge_offset
492
                    )
493
                )
494

495
            else:
496
                mf6_conc_array.append(
1✔
497
                    utils.concentration_m3_to_l(
498
                        self.mf6api.get_value(
499
                            self.mf6api.get_var_address(
500
                                "X",
501
                                f"{self.component_model_dict[c].upper()}",
502
                            )
503
                        )
504
                    )
505
                )
506
        c_dbl_vect = np.reshape(mf6_conc_array, self.nxyz * self.phreeqcbmi.ncomps)
1✔
507
        self.phreeqcbmi.SetConcentrations(c_dbl_vect)
1✔
508

509
        # set the kper and kstp
510
        self.phreeqcbmi._get_kper_kstp_from_mf6api(
1✔
511
            self.mf6api
512
        )  # FIXME: calling this func here is not ideal
513

514
        return c_dbl_vect
1✔
515

516
    def _solve(self) -> bool:
1✔
517
        """Alias for the solve method to provide backward compatibility"""
518
        return self.solve()
1✔
519

520
    def is_reactive_tstep(self) -> bool:
1✔
521
        """
522
        Check if the current timestep should be reactive based on configuration.
523

524
        Returns:
525
            bool: True if current timestep should be reactive, False otherwise
526
        """
527
        # Early return if not in reactive mode
528

529
        if not self.reactive:
1✔
530
            return False
1✔
531

532
        # Get current timestep
533
        current_tstep = [self.mf6api.kper, self.mf6api.kstp]
1✔
534

535
        # Check strategy
536
        if self.config.reactive['timing'] == 'all':
1✔
537
            return True
1✔
NEW
538
        elif self.config.reactive['timing'] == 'user':
×
NEW
539
            return current_tstep in self.config.reactive['tsteps']
×
540
        else:
541
            # Handle unknown strategy
NEW
542
            print(f"Warning: Unknown strategy '{self.config.reactive['timing']}'. Defaulting to reactive.")
×
543
            return True
×
544

545
    def set_kiter(self) -> int:
1✔
546
        if hasattr(self, "kiter"):
1✔
547
            self.kiter += 1
1✔
548
        else:
549
            self.kiter = 0
1✔
550
        return self.kiter
1✔
551
    def solve(self) -> bool:
1✔
552
        """Solve the model"""
553
        success = False  # initialize success flag
1✔
554
        sim_start = datetime.now()
1✔
555
        self._prepare_to_solve()
1✔
556

557
        # check sout was created
558
        assert self.selected_output._check_sout_exist(), f"{self.selected_output.sout_fname} not found"
1✔
559

560
        print("Starting Solution at {0}".format(sim_start.strftime(DT_FMT)))
1✔
561
        ctime = self._set_ctime()
1✔
562
        etime = self._set_etime()
1✔
563
        while ctime < etime:
1✔
564
            # self iteration counter
565
            self.set_kiter()
1✔
566
            # length of the current solve time
567
            dt = self._set_time_step()
1✔
568
            self.mf6api.prepare_time_step(dt)
1✔
569
            self.mf6api._solve_gwt()
1✔
570

571
            # get saturation
572
            self.get_saturation_from_mf6()
1✔
573
            # check_reactive_kstp()
574
            if self.is_reactive_tstep():
1✔
575
                c_dbl_vect = self._transfer_array_to_phreeqcrm()
1✔
576
                self._set_conc_at_current_kstep(c_dbl_vect)
1✔
577

578
                # Export ML feature arrays if option is on
579
                if self.ml_output:
1✔
580
                    self.selected_output.write_ml_arrays(self.current_iteration_conc,
×
581
                                                    self.kiter,
582
                                                    add_var_names=self.selected_output.feat_var,
583
                                                    fname='_features.csv'
584
                                                )
585

586
                if ctime == 0.0:
1✔
587
                    self.diffmask = np.ones(self.nxyz)
1✔
588
                else:
589
                    diffmask = get_conc_change_mask(
1✔
590
                        self.current_iteration_conc,
591
                        self.previous_iteration_conc,
592
                        self.phreeqcbmi.ncomps,
593
                        self.nxyz,
594
                        treshold=self.epsaqu,
595
                    )
596
                    self.diffmask = diffmask
1✔
597
                # solve reactions
598
                self.phreeqcbmi._solve_phreeqcrm(dt, diffmask=self.diffmask)
1✔
599
                c_dbl_vect = self._transfer_array_to_mf6()
1✔
600

601
                self._set_conc_at_previous_kstep(c_dbl_vect)
1✔
602

603
            self.mf6api.finalize_time_step()
1✔
604
            ctime = self._set_ctime()  # update the current time tracking
1✔
605
            if self.selected_output.get_selected_output_on:
1✔
606
                # get sout and update df
607
                self.selected_output._update_selected_output()
1✔
608
                # append current sout rows to file
609
                self.selected_output._append_to_soutdf_file()
1✔
610
                # Export ML target arrays if option is on
611
                if self.ml_output:
1✔
612
                    self.selected_output.write_ml_arrays(self.previous_iteration_conc,
×
613
                                            self.kiter,
614
                                            add_var_names=self.selected_output.target_var,
615
                                            fname='_targets.csv'
616
                                                )
617

618
        sim_end = datetime.now()
1✔
619
        td = (sim_end - sim_start).total_seconds() / 60.0
1✔
620

621
        self.mf6api._check_num_fails()
1✔
622

623
        # Clean up and close api objs
624
        try:
1✔
625
            self._finalize()
1✔
626
            success = True
1✔
627
            print(mrbeaker())
1✔
628
            print(
1✔
629
                "\nMR BEAKER IMPORTANT MESSAGE: MODEL RUN FINISHED BUT CHECK THE RESULTS .. THEY ARE PROLY RUBBISH\n"
630
            )
631
        except:
×
632
            print("MR BEAKER IMPORTANT MESSAGE: SOMETHING WENT WRONG. BUMMER\n")
×
633
            pass
×
634
        print(
1✔
635
            "Solution finished at {0}. Running time: {1:10.5G} mins".format(
636
                sim_end.strftime(DT_FMT), td
637
            )
638
        )
639
        return success
1✔
640

641

642
def get_less_than_zero_idx(arr):
1✔
643
    """Function to get the index of all occurrences of <0 in an array"""
644
    idx = np.where(arr < 0)
×
645
    return idx
×
646

647

648
def get_inactive_idx(arr: np.ndarray, val: float = 1e30):
1✔
649
    """Function to get the index of all occurrences of <0 in an array"""
650
    idx = list(np.where(arr >= val)[0])
×
651
    return idx
×
652

653

654
def get_conc_change_mask(
1✔
655
    ci: np.ndarray[np.float64],
656
    ck: np.ndarray[np.float64],
657
    ncomp: int,
658
    nxyz: int,
659
    treshold: float = 1e-10,
660
) -> np.ndarray[np.float64]:
661
    """Function to get the active-inactive cell mask for concentration change to inform phreeqc which cells to update"""
662
    # reshape arrays to 2D (nxyz, ncomp)
663
    ci = ci.reshape(nxyz, ncomp)
1✔
664
    ck = ck.reshape(nxyz, ncomp) + 1e-30
1✔
665

666
    # get the difference between the two arrays and divide by ci
667
    diff = np.abs((ci - ck.reshape(-1 * nxyz, ncomp)) / ci) < treshold
1✔
668
    diff = np.where(diff, 0, 1)
1✔
669
    diff = diff.sum(axis=1)
1✔
670

671
    # where values <0 put -1 else 1
672
    diff = np.where(diff == 0, 0, 1)
1✔
673
    return diff
1✔
674

675

676
def longest_common_substring(strings):
1✔
677
    """Function to find the longest common substring of a list of strings
678
    Used here to find the common "stem" of the GWT model names for matching
679
    with PhreeqcRM components.
680
    """
681
    if not strings:
1✔
682
        return ""
×
683

684
    # Start with the first string as a reference
685
    reference_string = strings[0]
1✔
686
    longest_lcs = ""
1✔
687

688
    # Iterate through all possible substrings of the reference string
689
    for i in range(len(reference_string)):
1✔
690
        for j in range(i + 1, len(reference_string) + 1):
1✔
691
            current_substring = reference_string[i:j]
1✔
692

693
            # Check if this substring exists in all other strings
694
            is_common = True
1✔
695
            for other_string in strings[1:]:
1✔
696
                if current_substring not in other_string:
1✔
697
                    is_common = False
1✔
698
                    break
1✔
699

700
            # If it's common and longer than the current longest, update
701
            if is_common and len(current_substring) > len(longest_lcs):
1✔
702
                longest_lcs = current_substring
×
703

704
    return longest_lcs
1✔
705

706

707
def mrbeaker() -> str:
1✔
708
    """ASCII art of Mr. Beaker"""
709

710
    from mf6rtm.assets import mrbeaker_path
1✔
711

712
    mr_beaker_image = Image.open(mrbeaker_path())
1✔
713

714
    # Resize the image to fit the terminal width
715
    terminal_width = 70  # Adjust this based on your terminal width
1✔
716
    aspect_ratio = mr_beaker_image.width / mr_beaker_image.height
1✔
717
    terminal_height = int(terminal_width / aspect_ratio * 0.5)
1✔
718
    mr_beaker_image = mr_beaker_image.resize((terminal_width, terminal_height))
1✔
719

720
    # Convert the image to grayscale
721
    mr_beaker_image = mr_beaker_image.convert("L")
1✔
722

723
    # Convert the grayscale image to ASCII art
724
    ascii_chars = "%,.?>#*+=-:."
1✔
725

726
    mrbeaker = ""
1✔
727
    for y in range(int(mr_beaker_image.height)):
1✔
728
        mrbeaker += "\n"
1✔
729
        for x in range(int(mr_beaker_image.width)):
1✔
730
            pixel_value = mr_beaker_image.getpixel((x, y))
1✔
731
            mrbeaker += ascii_chars[pixel_value // 64]
1✔
732
        # mrbeaker += "\n"
733

734
    return mrbeaker
1✔
735

736
def run_cmd():
1✔
737
    # get the current directory
738
    cwd = os.getcwd()
×
739
    # run the solve function
740
    solve(cwd)
×
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc