Why is Flopy giving the same results?

Question

I tried to use flopy, so I could try running some optimization procedures with python and modflow. Modflow requires a number of data to work with and we provide that info using different files.

We provide the input and flopy runs modflow

My trouble is that, flopy seems to disregard the input files and gives the same result, no matter what input I give.

Here's the code:

nper = 10
class BASreader:
    def __init__(self):
        self.ibound = None
        self.head = None

        with open("inps/bas/ibound", "r") as f:
            data = f.read().replace("-", " -").split()
            data = [int(x) for x in data]
            data = np.array(data).reshape(1, 71, 24)
            self.ibound = data

        with open("inps/bas/head", "r") as f:
            data = f.read().split()
            data = [float(x) for x in data]
            self.head = np.array(data).reshape(1, 71, 24)


class Modflow:
    def __init__(self):
        self.modelname = 'outs/gen1'
        self.mf = flopy.modflow.Modflow(self.modelname, exe_name=r'G:\Program Files (x86)\Visual MODFLOW\mf2005.exe')

        dis = flopy.modflow.ModflowDis.load("inps/yo.dis", self.mf)

        basreader = BASreader()
        bas = flopy.modflow.ModflowBas(self.mf, ibound=basreader.ibound, strt=basreader.head)
        self.prev_headdata = basreader.head

        wel = flopy.modflow.ModflowWel(self.mf, stress_period_data=WellProvider(nper).wells)

        fname = 'inps/yo.evt'
        fhandle = open(fname, 'r')
        packages = []
        ext_unit_dict = {22: flopy.utils.mfreadnam.NamData('EVT', fname, fhandle, packages)}
        evt = flopy.modflow.ModflowEvt.load(fhandle, self.mf, ext_unit_dict=ext_unit_dict)
        fhandle.close()

        rech = {}
        for x in range(nper):
            rech[x] = 1.44e-5
        rch = flopy.modflow.ModflowRch(self.mf, rech=rech)

        stress_period_data = {}
        for kper in range(nper):
            for kstp in range(int(nper/10)):
                stress_period_data[(kper, kstp)] = ['save head',
                                                    'save drawdown',
                                                    'save budget',
                                                    'print head',
                                                    'print drawdown',
                                                    'print budget']
        oc = flopy.modflow.ModflowOc(self.mf, stress_period_data=stress_period_data, compact=True)
        chd = flopy.modflow.ModflowChd.load("inps/yo.chd", self.mf)
        lpf = flopy.modflow.ModflowLpf(self.mf, hk=14.44, vka=14.44, ipakcb=53, sy=0.22)

        self.mf.write_input()

    def run(self):
        success, buff = self.mf.run_model(silent=True)
        headobj = bf.HeadFile(self.modelname + '.hds')
        newheaddata = headobj.get_alldata()[-1][0]
        return newheaddata

    mdflw = Modflow()
    mdflw.run()

Now, even if I change the EVT, RCH or WEL info, the results are same. I even tried to not include the above files, still the results were same. Any pointers?

Answer 1

This revision of your script results in different head results (since different groundwater recharge values are used). I used the bcf2ss MODFLOW-2005 input files as base files in the revised example.

Not entirely sure why your script wasn't working since your input files weren't available.

import os
import numpy as np
import flopy


class BASreader:
    def __init__(self):
        self.ibound = None
        self.head = None

        with open("inps/bas/ibound", "r") as f:
            data = f.read().replace("-", " -").split()
            data = [int(x) for x in data]
            data = np.array(data).reshape(2, 10, 15)
            self.ibound = data

        with open("inps/bas/head", "r") as f:
            data = f.read().split()
            data = [float(x) for x in data]
            self.head = np.array(data).reshape(2, 10, 15)


class Modflow:
    def __init__(self, ws='outs', rech_val=.0040, wel_data=None):
        self.model_ws = ws
        self.modelname = 'gen1'
        self.mf = flopy.modflow.Modflow(self.modelname, exe_name='mf2005',
                                        model_ws=self.model_ws)

        dis = flopy.modflow.ModflowDis.load("inps/yo.dis", self.mf, check=False)
        nper = dis.nper
        nstp = dis.nstp.array

        basreader = BASreader()
        bas = flopy.modflow.ModflowBas(self.mf, ibound=basreader.ibound,
                                       strt=basreader.head)
        self.prev_headdata = basreader.head
        bcf = flopy.modflow.ModflowBcf.load("inps/yo.bc6", self.mf)

        if wel_data is None:
            wel_data = [[1, 2, 3,  -35000.],
                        [1, 7, 3,  -35000.]]
        wel_spd = {1: wel_data}

        wel = flopy.modflow.ModflowWel(self.mf,
                                       stress_period_data=wel_spd)

        riv = flopy.modflow.ModflowRiv.load('inps/yo.riv', self.mf, check=False)

        rech = {}
        for x in range(nper):
            rech[x] = rech_val
        rch = flopy.modflow.ModflowRch(self.mf, rech=rech)

        pcg = flopy.modflow.ModflowPcg.load('inps/yo.pcg', self.mf)

        stress_period_data = {}
        for kper in range(nper):
            for kstp in range(nstp[kper]):
                stress_period_data[(kper, kstp)] = ['save head',
                                                    'save drawdown',
                                                    'save budget',
                                                    'print head',
                                                    'print drawdown',
                                                    'print budget']
        oc = flopy.modflow.ModflowOc(self.mf,
                                     stress_period_data=stress_period_data,
                                     compact=True)

        self.mf.write_input()

    def run(self):
        success, buff = self.mf.run_model(silent=True)
        fpth = os.path.join(self.model_ws, self.modelname + '.hds')
        headobj = flopy.utils.HeadFile(fpth)
        newheaddata = headobj.get_data(idx=1)
        return newheaddata


if __name__ == "__main__":
    m1 = Modflow()
    h1 = m1.run()

    m2 = Modflow(rech_val=.0041)
    h2 = m2.run()

    assert np.array_equal(h1, h2), 'h1 does not equal h2'