diff --git a/.docs/conf.py b/.docs/conf.py index c0b275c6c1..4639f4ddb6 100644 --- a/.docs/conf.py +++ b/.docs/conf.py @@ -79,9 +79,9 @@ # -- Project information ----------------------------------------------------- -project = 'flopy' -copyright = '2020, Bakker, Mark, Post, Vincent, Langevin, C. D., Hughes, J. D., White, J. T., Leaf, A. T., Paulinski, S. R., Larsen, J. D., Toews, M. W., Morway, E. D., Bellino, J. C., Starn, J. J., and Fienen, M. N.' -author = 'Bakker, Mark, Post, Vincent, Langevin, C. D., Hughes, J. D., White, J. T., Leaf, A. T., Paulinski, S. R., Larsen, J. D., Toews, M. W., Morway, E. D., Bellino, J. C., Starn, J. J., and Fienen, M. N.' +project = "flopy" +copyright = "2020, Bakker, Mark, Post, Vincent, Langevin, C. D., Hughes, J. D., White, J. T., Leaf, A. T., Paulinski, S. R., Larsen, J. D., Toews, M. W., Morway, E. D., Bellino, J. C., Starn, J. J., and Fienen, M. N." +author = "Bakker, Mark, Post, Vincent, Langevin, C. D., Hughes, J. D., White, J. T., Leaf, A. T., Paulinski, S. R., Larsen, J. D., Toews, M. W., Morway, E. D., Bellino, J. C., Starn, J. J., and Fienen, M. N." # The version. version = __version__ @@ -112,7 +112,6 @@ ] - # Add any paths that contain templates here, relative to this directory. templates_path = ["_templates"] diff --git a/.docs/pysrc/tutorial1.py b/.docs/pysrc/tutorial1.py index 0186b7187c..c611220f9f 100644 --- a/.docs/pysrc/tutorial1.py +++ b/.docs/pysrc/tutorial1.py @@ -3,80 +3,82 @@ import numpy as np import flopy -#Assign name and create modflow model object -modelname = 'tutorial1' +# Assign name and create modflow model object +modelname = "tutorial1" exe_name = os.path.join("..", ".bin", "mf2005") mf = flopy.modflow.Modflow(modelname, exe_name=exe_name) -#model domain and grid definition -Lx = 1000. -Ly = 1000. -ztop = 0. -zbot = -50. +# model domain and grid definition +Lx = 1000.0 +Ly = 1000.0 +ztop = 0.0 +zbot = -50.0 nlay = 1 nrow = 10 ncol = 10 -delr = Lx/ncol -delc = Ly/nrow +delr = Lx / ncol +delc = Ly / nrow delv = (ztop - zbot) / nlay botm = np.linspace(ztop, zbot, nlay + 1) # Create the discretization object -dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, - top=ztop, botm=botm[1:]) +dis = flopy.modflow.ModflowDis( + mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm[1:] +) # Variables for the BAS package ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) ibound[:, :, 0] = -1 ibound[:, :, -1] = -1 strt = np.ones((nlay, nrow, ncol), dtype=np.float32) -strt[:, :, 0] = 10. -strt[:, :, -1] = 0. +strt[:, :, 0] = 10.0 +strt[:, :, -1] = 0.0 bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) -#Add LPF package to the MODFLOW model -lpf = flopy.modflow.ModflowLpf(mf, hk=10., vka=10., ipakcb=53) +# Add LPF package to the MODFLOW model +lpf = flopy.modflow.ModflowLpf(mf, hk=10.0, vka=10.0, ipakcb=53) # Add OC package to the MODFLOW model -spd = {(0, 0): ['print head', 'print budget', 'save head', 'save budget']} +spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]} oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True) -#Add PCG package to the MODFLOW model +# Add PCG package to the MODFLOW model pcg = flopy.modflow.ModflowPcg(mf) -#Write the MODFLOW model input files +# Write the MODFLOW model input files mf.write_input() -#Run the MODFLOW model +# Run the MODFLOW model mf.run_model() import matplotlib.pyplot as plt import flopy.utils.binaryfile as bf -plt.subplot(1,1,1,aspect='equal') -hds = bf.HeadFile(modelname+'.hds') + +plt.subplot(1, 1, 1, aspect="equal") +hds = bf.HeadFile(modelname + ".hds") head = hds.get_data(totim=1.0) -levels = np.arange(1,10,1) -extent = (delr/2., Lx - delc/2., Ly - delc/2., delc/2.) +levels = np.arange(1, 10, 1) +extent = (delr / 2.0, Lx - delc / 2.0, Ly - delc / 2.0, delc / 2.0) plt.contour(head[0, :, :], levels=levels, extent=extent) -plt.savefig('tutorial1a.png') +plt.savefig("tutorial1a.png") -fig = plt.figure(figsize=(10,10)) -ax = fig.add_subplot(1, 1, 1, aspect='equal') +fig = plt.figure(figsize=(10, 10)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") -hds = bf.HeadFile(modelname + '.hds') +hds = bf.HeadFile(modelname + ".hds") times = hds.get_times() head = hds.get_data(totim=times[-1]) levels = np.linspace(0, 10, 11) -cbb = bf.CellBudgetFile(modelname + '.cbc') +cbb = bf.CellBudgetFile(modelname + ".cbc") kstpkper_list = cbb.get_kstpkper() -frf = cbb.get_data(text='FLOW RIGHT FACE', totim=times[-1])[0] -fff = cbb.get_data(text='FLOW FRONT FACE', totim=times[-1])[0] +frf = cbb.get_data(text="FLOW RIGHT FACE", totim=times[-1])[0] +fff = cbb.get_data(text="FLOW FRONT FACE", totim=times[-1])[0] pmv = flopy.plot.PlotMapView(model=mf, layer=0) qm = pmv.plot_ibound() lc = pmv.plot_grid() cs = pmv.contour_array(head, levels=levels) quiver = pmv.plot_discharge(frf, fff, head=head) -plt.savefig('tutorial1b.png') +plt.savefig("tutorial1b.png") diff --git a/.docs/pysrc/tutorial2.py b/.docs/pysrc/tutorial2.py index 11b340bbd5..5fc9396457 100644 --- a/.docs/pysrc/tutorial2.py +++ b/.docs/pysrc/tutorial2.py @@ -3,10 +3,10 @@ import flopy # Model domain and grid definition -Lx = 1000. -Ly = 1000. -ztop = 10. -zbot = -50. +Lx = 1000.0 +Ly = 1000.0 +ztop = 10.0 +zbot = -50.0 nlay = 1 nrow = 10 ncol = 10 @@ -14,16 +14,16 @@ delc = Ly / nrow delv = (ztop - zbot) / nlay botm = np.linspace(ztop, zbot, nlay + 1) -hk = 1. -vka = 1. +hk = 1.0 +vka = 1.0 sy = 0.1 -ss = 1.e-4 +ss = 1.0e-4 laytyp = 1 # Variables for the BAS package # Note that changes from the previous tutorial! ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) -strt = 10. * np.ones((nlay, nrow, ncol), dtype=np.float32) +strt = 10.0 * np.ones((nlay, nrow, ncol), dtype=np.float32) # Time step parameters nper = 3 @@ -32,21 +32,32 @@ steady = [True, False, False] # Flopy objects -modelname = 'tutorial2' +modelname = "tutorial2" exe_name = os.path.join("..", ".bin", "mf2005") mf = flopy.modflow.Modflow(modelname, exe_name=exe_name) -dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, - top=ztop, botm=botm[1:], - nper=nper, perlen=perlen, nstp=nstp, - steady=steady) +dis = flopy.modflow.ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=ztop, + botm=botm[1:], + nper=nper, + perlen=perlen, + nstp=nstp, + steady=steady, +) bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) -lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp, - ipakcb=53) +lpf = flopy.modflow.ModflowLpf( + mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp, ipakcb=53 +) pcg = flopy.modflow.ModflowPcg(mf) # Make list for stress period 1 -stageleft = 10. -stageright = 10. +stageleft = 10.0 +stageright = 10.0 bound_sp1 = [] for il in range(nlay): condleft = hk * (stageleft - zbot) * delc @@ -54,11 +65,11 @@ for ir in range(nrow): bound_sp1.append([il, ir, 0, stageleft, condleft]) bound_sp1.append([il, ir, ncol - 1, stageright, condright]) -print('Adding ', len(bound_sp1), 'GHBs for stress period 1.') +print("Adding ", len(bound_sp1), "GHBs for stress period 1.") # Make list for stress period 2 -stageleft = 10. -stageright = 0. +stageleft = 10.0 +stageright = 0.0 condleft = hk * (stageleft - zbot) * delc condright = hk * (stageright - zbot) * delc bound_sp2 = [] @@ -66,7 +77,7 @@ for ir in range(nrow): bound_sp2.append([il, ir, 0, stageleft, condleft]) bound_sp2.append([il, ir, ncol - 1, stageright, condright]) -print('Adding ', len(bound_sp2), 'GHBs for stress period 2.') +print("Adding ", len(bound_sp2), "GHBs for stress period 2.") # We do not need to add a dictionary entry for stress period 3. # Flopy will automatically take the list from stress period 2 and apply it @@ -78,9 +89,9 @@ # Create the well package # Remember to use zero-based layer, row, column indices! -pumping_rate = -100. -wel_sp1 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.]] -wel_sp2 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.]] +pumping_rate = -100.0 +wel_sp1 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.0]] +wel_sp2 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.0]] wel_sp3 = [[0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate]] stress_period_data = {0: wel_sp1, 1: wel_sp2, 2: wel_sp3} wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data) @@ -88,13 +99,16 @@ 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 budget'] -oc = flopy.modflow.ModflowOc(mf, stress_period_data=stress_period_data, - compact=True) + stress_period_data[(kper, kstp)] = [ + "save head", + "save drawdown", + "save budget", + "print head", + "print budget", + ] +oc = flopy.modflow.ModflowOc( + mf, stress_period_data=stress_period_data, compact=True +) # Write the model input files mf.write_input() @@ -102,7 +116,7 @@ # Run the model success, mfoutput = mf.run_model(silent=True, pause=False, report=True) if not success: - raise Exception('MODFLOW did not terminate normally.') + raise Exception("MODFLOW did not terminate normally.") # Imports @@ -110,66 +124,72 @@ import flopy.utils.binaryfile as bf # Create the headfile and budget file objects -headobj = bf.HeadFile(modelname+'.hds') +headobj = bf.HeadFile(modelname + ".hds") times = headobj.get_times() -cbb = bf.CellBudgetFile(modelname+'.cbc') +cbb = bf.CellBudgetFile(modelname + ".cbc") # Setup contour parameters levels = np.linspace(0, 10, 11) -extent = (delr/2., Lx - delr/2., delc/2., Ly - delc/2.) -print('Levels: ', levels) -print('Extent: ', extent) +extent = (delr / 2.0, Lx - delr / 2.0, delc / 2.0, Ly - delc / 2.0) +print("Levels: ", levels) +print("Extent: ", extent) # Well point -wpt = ((float(ncol/2)-0.5)*delr, (float(nrow/2-1)+0.5)*delc) -wpt = (450., 550.) +wpt = ((float(ncol / 2) - 0.5) * delr, (float(nrow / 2 - 1) + 0.5) * delc) +wpt = (450.0, 550.0) # Make the plots mytimes = [1.0, 101.0, 201.0] for iplot, time in enumerate(mytimes): - print('*****Processing time: ', time) - head = headobj.get_data(totim=time) - #Print statistics - print('Head statistics') - print(' min: ', head.min()) - print(' max: ', head.max()) - print(' std: ', head.std()) - - # Extract flow right face and flow front face - frf = cbb.get_data(text='FLOW RIGHT FACE', totim=time)[0] - fff = cbb.get_data(text='FLOW FRONT FACE', totim=time)[0] - - #Create the plot - f = plt.figure() - plt.subplot(1, 1, 1, aspect='equal') - plt.title('stress period ' + str(iplot + 1)) - - - modelmap = flopy.plot.PlotMapView(model=mf, layer=0) - qm = modelmap.plot_ibound() - lc = modelmap.plot_grid() - qm = modelmap.plot_bc('GHB', alpha=0.5) - cs = modelmap.contour_array(head, levels=levels) - plt.clabel(cs, inline=1, fontsize=10, fmt='%1.1f') - quiver = modelmap.plot_discharge(frf, fff, head=head) - - - mfc = 'None' - if (iplot+1) == len(mytimes): - mfc='black' - plt.plot(wpt[0], wpt[1], lw=0, marker='o', markersize=8, - markeredgewidth=0.5, - markeredgecolor='black', markerfacecolor=mfc, zorder=9) - plt.text(wpt[0]+25, wpt[1]-25, 'well', size=12, zorder=12) - plt.savefig('tutorial2-{}.png'.format(iplot)) + print("*****Processing time: ", time) + head = headobj.get_data(totim=time) + # Print statistics + print("Head statistics") + print(" min: ", head.min()) + print(" max: ", head.max()) + print(" std: ", head.std()) + + # Extract flow right face and flow front face + frf = cbb.get_data(text="FLOW RIGHT FACE", totim=time)[0] + fff = cbb.get_data(text="FLOW FRONT FACE", totim=time)[0] + + # Create the plot + f = plt.figure() + plt.subplot(1, 1, 1, aspect="equal") + plt.title("stress period " + str(iplot + 1)) + + modelmap = flopy.plot.PlotMapView(model=mf, layer=0) + qm = modelmap.plot_ibound() + lc = modelmap.plot_grid() + qm = modelmap.plot_bc("GHB", alpha=0.5) + cs = modelmap.contour_array(head, levels=levels) + plt.clabel(cs, inline=1, fontsize=10, fmt="%1.1f") + quiver = modelmap.plot_discharge(frf, fff, head=head) + + mfc = "None" + if (iplot + 1) == len(mytimes): + mfc = "black" + plt.plot( + wpt[0], + wpt[1], + lw=0, + marker="o", + markersize=8, + markeredgewidth=0.5, + markeredgecolor="black", + markerfacecolor=mfc, + zorder=9, + ) + plt.text(wpt[0] + 25, wpt[1] - 25, "well", size=12, zorder=12) + plt.savefig("tutorial2-{}.png".format(iplot)) # Plot the head versus time idx = (0, int(nrow / 2) - 1, int(ncol / 2) - 1) ts = headobj.get_ts(idx) plt.subplot(1, 1, 1) -ttl = 'Head at cell ({0},{1},{2})'.format(idx[0] + 1, idx[1] + 1, idx[2] + 1) +ttl = "Head at cell ({0},{1},{2})".format(idx[0] + 1, idx[1] + 1, idx[2] + 1) plt.title(ttl) -plt.xlabel('time') -plt.ylabel('head') -plt.plot(ts[:, 0], ts[:, 1], 'bo-') -plt.savefig('tutorial2-ts.png') +plt.xlabel("time") +plt.ylabel("head") +plt.plot(ts[:, 0], ts[:, 1], "bo-") +plt.savefig("tutorial2-ts.png") diff --git a/examples/Notebooks/setup_pmv_demo.py b/examples/Notebooks/setup_pmv_demo.py index 268ae18015..09d48caec1 100644 --- a/examples/Notebooks/setup_pmv_demo.py +++ b/examples/Notebooks/setup_pmv_demo.py @@ -11,21 +11,22 @@ def run(): try: import flopy except: - fpth = os.path.abspath(os.path.join('..', '..')) + fpth = os.path.abspath(os.path.join("..", "..")) sys.path.append(fpth) import flopy print(sys.version) - print('numpy version: {}'.format(np.__version__)) - print('matplotlib version: {}'.format(mpl.__version__)) - print('flopy version: {}'.format(flopy.__version__)) + print("numpy version: {}".format(np.__version__)) + print("matplotlib version: {}".format(mpl.__version__)) + print("flopy version: {}".format(flopy.__version__)) if not os.path.exists("data"): os.mkdir("data") from flopy.utils.gridgen import Gridgen - Lx = 10000. - Ly = 10500. + + Lx = 10000.0 + Ly = 10500.0 nlay = 3 nrow = 21 ncol = 20 @@ -35,180 +36,267 @@ def run(): botm = [220, 200, 0] ms = flopy.modflow.Modflow() - dis5 = flopy.modflow.ModflowDis(ms, nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, - delc=delc, top=top, botm=botm) - - model_name = 'mp7p2' - model_ws = os.path.join('data', 'mp7_ex2', 'mf6') - gridgen_ws = os.path.join(model_ws, 'gridgen') + dis5 = flopy.modflow.ModflowDis( + ms, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + model_name = "mp7p2" + model_ws = os.path.join("data", "mp7_ex2", "mf6") + gridgen_ws = os.path.join(model_ws, "gridgen") g = Gridgen(dis5, model_ws=gridgen_ws) - rf0shp = os.path.join(gridgen_ws, 'rf0') + rf0shp = os.path.join(gridgen_ws, "rf0") xmin = 7 * delr xmax = 12 * delr ymin = 8 * delc ymax = 13 * delc - rfpoly = [[[(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax), (xmin, ymin)]]] - g.add_refinement_features(rfpoly, 'polygon', 1, range(nlay)) - - rf1shp = os.path.join(gridgen_ws, 'rf1') + rfpoly = [ + [ + [ + (xmin, ymin), + (xmax, ymin), + (xmax, ymax), + (xmin, ymax), + (xmin, ymin), + ] + ] + ] + g.add_refinement_features(rfpoly, "polygon", 1, range(nlay)) + + rf1shp = os.path.join(gridgen_ws, "rf1") xmin = 8 * delr xmax = 11 * delr ymin = 9 * delc ymax = 12 * delc - rfpoly = [[[(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax), (xmin, ymin)]]] - g.add_refinement_features(rfpoly, 'polygon', 2, range(nlay)) - - rf2shp = os.path.join(gridgen_ws, 'rf2') + rfpoly = [ + [ + [ + (xmin, ymin), + (xmax, ymin), + (xmax, ymax), + (xmin, ymax), + (xmin, ymin), + ] + ] + ] + g.add_refinement_features(rfpoly, "polygon", 2, range(nlay)) + + rf2shp = os.path.join(gridgen_ws, "rf2") xmin = 9 * delr xmax = 10 * delr ymin = 10 * delc ymax = 11 * delc - rfpoly = [[[(xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax), (xmin, ymin)]]] - g.add_refinement_features(rfpoly, 'polygon', 3, range(nlay)) + rfpoly = [ + [ + [ + (xmin, ymin), + (xmax, ymin), + (xmax, ymax), + (xmin, ymax), + (xmin, ymin), + ] + ] + ] + g.add_refinement_features(rfpoly, "polygon", 3, range(nlay)) g.build(verbose=False) gridprops = g.get_gridprops_disv() - ncpl = gridprops['ncpl'] - top = gridprops['top'] - botm = gridprops['botm'] - nvert = gridprops['nvert'] - vertices = gridprops['vertices'] - cell2d = gridprops['cell2d'] + ncpl = gridprops["ncpl"] + top = gridprops["top"] + botm = gridprops["botm"] + nvert = gridprops["nvert"] + vertices = gridprops["vertices"] + cell2d = gridprops["cell2d"] # cellxy = gridprops['cellxy'] # create simulation - sim = flopy.mf6.MFSimulation(sim_name=model_name, version='mf6', exe_name='mf6', - sim_ws=model_ws) + sim = flopy.mf6.MFSimulation( + sim_name=model_name, version="mf6", exe_name="mf6", sim_ws=model_ws + ) # create tdis package tdis_rc = [(1000.0, 1, 1.0)] - tdis = flopy.mf6.ModflowTdis(sim, pname='tdis', time_units='DAYS', - perioddata=tdis_rc) + tdis = flopy.mf6.ModflowTdis( + sim, pname="tdis", time_units="DAYS", perioddata=tdis_rc + ) # create gwf model - gwf = flopy.mf6.ModflowGwf(sim, modelname=model_name, - model_nam_file='{}.nam'.format(model_name)) + gwf = flopy.mf6.ModflowGwf( + sim, modelname=model_name, model_nam_file="{}.nam".format(model_name) + ) gwf.name_file.save_flows = True # create iterative model solution and register the gwf model with it - ims = flopy.mf6.ModflowIms(sim, pname='ims', print_option='SUMMARY', - complexity='SIMPLE', outer_hclose=1.e-5, - outer_maximum=100, under_relaxation='NONE', - inner_maximum=100, inner_hclose=1.e-6, - rcloserecord=0.1, linear_acceleration='BICGSTAB', - scaling_method='NONE', reordering_method='NONE', - relaxation_factor=0.99) + ims = flopy.mf6.ModflowIms( + sim, + pname="ims", + print_option="SUMMARY", + complexity="SIMPLE", + outer_hclose=1.0e-5, + outer_maximum=100, + under_relaxation="NONE", + inner_maximum=100, + inner_hclose=1.0e-6, + rcloserecord=0.1, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=0.99, + ) sim.register_ims_package(ims, [gwf.name]) # disv - disv = flopy.mf6.ModflowGwfdisv(gwf, nlay=nlay, ncpl=ncpl, - top=top, botm=botm, - nvert=nvert, vertices=vertices, - cell2d=cell2d) + disv = flopy.mf6.ModflowGwfdisv( + gwf, + nlay=nlay, + ncpl=ncpl, + top=top, + botm=botm, + nvert=nvert, + vertices=vertices, + cell2d=cell2d, + ) # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, pname='ic', strt=320.) + ic = flopy.mf6.ModflowGwfic(gwf, pname="ic", strt=320.0) # node property flow - npf = flopy.mf6.ModflowGwfnpf(gwf, xt3doptions=[('xt3d')], - save_specific_discharge=True, - icelltype=[1, 0, 0], - k=[50.0, 0.01, 200.0], - k33=[10., 0.01, 20.]) + npf = flopy.mf6.ModflowGwfnpf( + gwf, + xt3doptions=[("xt3d")], + save_specific_discharge=True, + icelltype=[1, 0, 0], + k=[50.0, 0.01, 200.0], + k33=[10.0, 0.01, 20.0], + ) # wel - wellpoints = [(4750., 5250.)] - welcells = g.intersect(wellpoints, 'point', 0) + wellpoints = [(4750.0, 5250.0)] + welcells = g.intersect(wellpoints, "point", 0) # welspd = flopy.mf6.ModflowGwfwel.stress_period_data.empty(gwf, maxbound=1, aux_vars=['iface']) - welspd = [[(2, icpl), -150000, 0] for icpl in welcells['nodenumber']] - wel = flopy.mf6.ModflowGwfwel(gwf, print_input=True, - auxiliary=[('iface',)], - stress_period_data=welspd) + welspd = [[(2, icpl), -150000, 0] for icpl in welcells["nodenumber"]] + wel = flopy.mf6.ModflowGwfwel( + gwf, + print_input=True, + auxiliary=[("iface",)], + stress_period_data=welspd, + ) # rch aux = [np.ones(ncpl, dtype=np.int) * 6] - rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=0.005, - auxiliary=[('iface',)], - aux={0: [6]}) + rch = flopy.mf6.ModflowGwfrcha( + gwf, recharge=0.005, auxiliary=[("iface",)], aux={0: [6]} + ) # riv - riverline = [[[(Lx - 1., Ly), (Lx - 1., 0.)]]] - rivcells = g.intersect(riverline, 'line', 0) - rivspd = [[(0, icpl), 320., 100000., 318] for icpl in rivcells['nodenumber']] + riverline = [[[(Lx - 1.0, Ly), (Lx - 1.0, 0.0)]]] + rivcells = g.intersect(riverline, "line", 0) + rivspd = [ + [(0, icpl), 320.0, 100000.0, 318] for icpl in rivcells["nodenumber"] + ] riv = flopy.mf6.ModflowGwfriv(gwf, stress_period_data=rivspd) # output control - oc = flopy.mf6.ModflowGwfoc(gwf, pname='oc', budget_filerecord='{}.cbb'.format(model_name), - head_filerecord='{}.hds'.format(model_name), - headprintrecord=[('COLUMNS', 10, 'WIDTH', 15, - 'DIGITS', 6, 'GENERAL')], - saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')], - printrecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')]) + oc = flopy.mf6.ModflowGwfoc( + gwf, + pname="oc", + budget_filerecord="{}.cbb".format(model_name), + head_filerecord="{}.hds".format(model_name), + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) sim.write_simulation() sim.run_simulation() - mp_namea = model_name + 'a_mp' - mp_nameb = model_name + 'b_mp' - - pcoord = np.array([[0.000, 0.125, 0.500], - [0.000, 0.375, 0.500], - [0.000, 0.625, 0.500], - [0.000, 0.875, 0.500], - [1.000, 0.125, 0.500], - [1.000, 0.375, 0.500], - [1.000, 0.625, 0.500], - [1.000, 0.875, 0.500], - [0.125, 0.000, 0.500], - [0.375, 0.000, 0.500], - [0.625, 0.000, 0.500], - [0.875, 0.000, 0.500], - [0.125, 1.000, 0.500], - [0.375, 1.000, 0.500], - [0.625, 1.000, 0.500], - [0.875, 1.000, 0.500]]) - nodew = gwf.disv.ncpl.array * 2 + welcells['nodenumber'][0] + mp_namea = model_name + "a_mp" + mp_nameb = model_name + "b_mp" + + pcoord = np.array( + [ + [0.000, 0.125, 0.500], + [0.000, 0.375, 0.500], + [0.000, 0.625, 0.500], + [0.000, 0.875, 0.500], + [1.000, 0.125, 0.500], + [1.000, 0.375, 0.500], + [1.000, 0.625, 0.500], + [1.000, 0.875, 0.500], + [0.125, 0.000, 0.500], + [0.375, 0.000, 0.500], + [0.625, 0.000, 0.500], + [0.875, 0.000, 0.500], + [0.125, 1.000, 0.500], + [0.375, 1.000, 0.500], + [0.625, 1.000, 0.500], + [0.875, 1.000, 0.500], + ] + ) + nodew = gwf.disv.ncpl.array * 2 + welcells["nodenumber"][0] plocs = [nodew for i in range(pcoord.shape[0])] # create particle data - pa = flopy.modpath.ParticleData(plocs, structured=False, - localx=pcoord[:, 0], - localy=pcoord[:, 1], - localz=pcoord[:, 2], - drape=0) + pa = flopy.modpath.ParticleData( + plocs, + structured=False, + localx=pcoord[:, 0], + localy=pcoord[:, 1], + localz=pcoord[:, 2], + drape=0, + ) # create backward particle group - fpth = mp_namea + '.sloc' - pga = flopy.modpath.ParticleGroup(particlegroupname='BACKWARD1', particledata=pa, - filename=fpth) - - facedata = flopy.modpath.FaceDataType(drape=0, - verticaldivisions1=10, horizontaldivisions1=10, - verticaldivisions2=10, horizontaldivisions2=10, - verticaldivisions3=10, horizontaldivisions3=10, - verticaldivisions4=10, horizontaldivisions4=10, - rowdivisions5=0, columndivisions5=0, - rowdivisions6=4, columndivisions6=4) + fpth = mp_namea + ".sloc" + pga = flopy.modpath.ParticleGroup( + particlegroupname="BACKWARD1", particledata=pa, filename=fpth + ) + + facedata = flopy.modpath.FaceDataType( + drape=0, + verticaldivisions1=10, + horizontaldivisions1=10, + verticaldivisions2=10, + horizontaldivisions2=10, + verticaldivisions3=10, + horizontaldivisions3=10, + verticaldivisions4=10, + horizontaldivisions4=10, + rowdivisions5=0, + columndivisions5=0, + rowdivisions6=4, + columndivisions6=4, + ) pb = flopy.modpath.NodeParticleData(subdivisiondata=facedata, nodes=nodew) # create forward particle group - fpth = mp_nameb + '.sloc' - pgb = flopy.modpath.ParticleGroupNodeTemplate(particlegroupname='BACKWARD2', - particledata=pb, - filename=fpth) + fpth = mp_nameb + ".sloc" + pgb = flopy.modpath.ParticleGroupNodeTemplate( + particlegroupname="BACKWARD2", particledata=pb, filename=fpth + ) # create modpath files - mp = flopy.modpath.Modpath7(modelname=mp_namea, flowmodel=gwf, - exe_name='mp7', model_ws=model_ws) + mp = flopy.modpath.Modpath7( + modelname=mp_namea, flowmodel=gwf, exe_name="mp7", model_ws=model_ws + ) flopy.modpath.Modpath7Bas(mp, porosity=0.1) - flopy.modpath.Modpath7Sim(mp, simulationtype='combined', - trackingdirection='backward', - weaksinkoption='pass_through', - weaksourceoption='pass_through', - referencetime=0., - stoptimeoption='extend', - timepointdata=[500, 1000.], - particlegroups=pga) + flopy.modpath.Modpath7Sim( + mp, + simulationtype="combined", + trackingdirection="backward", + weaksinkoption="pass_through", + weaksourceoption="pass_through", + referencetime=0.0, + stoptimeoption="extend", + timepointdata=[500, 1000.0], + particlegroups=pga, + ) # write modpath datasets mp.write_input() @@ -217,16 +305,20 @@ def run(): mp.run_model() # create modpath files - mp = flopy.modpath.Modpath7(modelname=mp_nameb, flowmodel=gwf, - exe_name='mp7', model_ws=model_ws) + mp = flopy.modpath.Modpath7( + modelname=mp_nameb, flowmodel=gwf, exe_name="mp7", model_ws=model_ws + ) flopy.modpath.Modpath7Bas(mp, porosity=0.1) - flopy.modpath.Modpath7Sim(mp, simulationtype='endpoint', - trackingdirection='backward', - weaksinkoption='pass_through', - weaksourceoption='pass_through', - referencetime=0., - stoptimeoption='extend', - particlegroups=pgb) + flopy.modpath.Modpath7Sim( + mp, + simulationtype="endpoint", + trackingdirection="backward", + weaksinkoption="pass_through", + weaksourceoption="pass_through", + referencetime=0.0, + stoptimeoption="extend", + particlegroups=pgb, + ) # write modpath datasets mp.write_input() @@ -235,5 +327,6 @@ def run(): mp.run_model() return + if __name__ == "__main__": - run() \ No newline at end of file + run() diff --git a/examples/Tutorials/Tutorial01/tutorial01.py b/examples/Tutorials/Tutorial01/tutorial01.py index 940ed41da7..f08ab8b637 100644 --- a/examples/Tutorials/Tutorial01/tutorial01.py +++ b/examples/Tutorials/Tutorial01/tutorial01.py @@ -1,42 +1,42 @@ - import numpy as np import flopy # Assign name and create modflow model object -modelname = 'tutorial1' -mf = flopy.modflow.Modflow(modelname, exe_name='mf2005') +modelname = "tutorial1" +mf = flopy.modflow.Modflow(modelname, exe_name="mf2005") # Model domain and grid definition -Lx = 1000. -Ly = 1000. -ztop = 0. -zbot = -50. +Lx = 1000.0 +Ly = 1000.0 +ztop = 0.0 +zbot = -50.0 nlay = 1 nrow = 10 ncol = 10 -delr = Lx/ncol -delc = Ly/nrow +delr = Lx / ncol +delc = Ly / nrow delv = (ztop - zbot) / nlay botm = np.linspace(ztop, zbot, nlay + 1) # Create the discretization object -dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, - top=ztop, botm=botm[1:]) +dis = flopy.modflow.ModflowDis( + mf, nlay, nrow, ncol, delr=delr, delc=delc, top=ztop, botm=botm[1:] +) # Variables for the BAS package ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) ibound[:, :, 0] = -1 ibound[:, :, -1] = -1 strt = np.ones((nlay, nrow, ncol), dtype=np.float32) -strt[:, :, 0] = 10. -strt[:, :, -1] = 0. +strt[:, :, 0] = 10.0 +strt[:, :, -1] = 0.0 bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) # Add LPF package to the MODFLOW model -lpf = flopy.modflow.ModflowLpf(mf, hk=10., vka=10., ipakcb=53) +lpf = flopy.modflow.ModflowLpf(mf, hk=10.0, vka=10.0, ipakcb=53) # Add OC package to the MODFLOW model -spd = {(0, 0): ['print head', 'print budget', 'save head', 'save budget']} +spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]} oc = flopy.modflow.ModflowOc(mf, stress_period_data=spd, compact=True) # Add PCG package to the MODFLOW model @@ -52,33 +52,30 @@ import matplotlib.pyplot as plt import flopy.utils.binaryfile as bf -plt.subplot(1, 1, 1, aspect='equal') -hds = bf.HeadFile(modelname + '.hds') +plt.subplot(1, 1, 1, aspect="equal") +hds = bf.HeadFile(modelname + ".hds") head = hds.get_data(totim=1.0) levels = np.arange(1, 10, 1) -extent = (delr / 2., Lx - delr / 2., Ly - delc / 2., delc / 2.) +extent = (delr / 2.0, Lx - delr / 2.0, Ly - delc / 2.0, delc / 2.0) plt.contour(head[0, :, :], levels=levels, extent=extent) -plt.savefig('tutorial1a.png') +plt.savefig("tutorial1a.png") -fig = plt.figure(figsize=(10,10)) -ax = fig.add_subplot(1, 1, 1, aspect='equal') +fig = plt.figure(figsize=(10, 10)) +ax = fig.add_subplot(1, 1, 1, aspect="equal") -hds = bf.HeadFile(modelname+'.hds') +hds = bf.HeadFile(modelname + ".hds") times = hds.get_times() head = hds.get_data(totim=times[-1]) levels = np.linspace(0, 10, 11) -cbb = bf.CellBudgetFile(modelname+'.cbc') +cbb = bf.CellBudgetFile(modelname + ".cbc") kstpkper_list = cbb.get_kstpkper() -frf = cbb.get_data(text='FLOW RIGHT FACE', totim=times[-1])[0] -fff = cbb.get_data(text='FLOW FRONT FACE', totim=times[-1])[0] +frf = cbb.get_data(text="FLOW RIGHT FACE", totim=times[-1])[0] +fff = cbb.get_data(text="FLOW FRONT FACE", totim=times[-1])[0] modelmap = flopy.plot.ModelMap(model=mf, layer=0) qm = modelmap.plot_ibound() lc = modelmap.plot_grid() cs = modelmap.contour_array(head, levels=levels) quiver = modelmap.plot_discharge(frf, fff, head=head) -plt.savefig('tutorial1b.png') - - - +plt.savefig("tutorial1b.png") diff --git a/examples/Tutorials/Tutorial02/tutorial02.py b/examples/Tutorials/Tutorial02/tutorial02.py index 588d9ec4d3..af0c837819 100644 --- a/examples/Tutorials/Tutorial02/tutorial02.py +++ b/examples/Tutorials/Tutorial02/tutorial02.py @@ -1,12 +1,11 @@ - import numpy as np import flopy # Model domain and grid definition -Lx = 1000. -Ly = 1000. -ztop = 10. -zbot = -50. +Lx = 1000.0 +Ly = 1000.0 +ztop = 10.0 +zbot = -50.0 nlay = 1 nrow = 10 ncol = 10 @@ -14,16 +13,16 @@ delc = Ly / nrow delv = (ztop - zbot) / nlay botm = np.linspace(ztop, zbot, nlay + 1) -hk = 1. -vka = 1. +hk = 1.0 +vka = 1.0 sy = 0.1 -ss = 1.e-4 +ss = 1.0e-4 laytyp = 1 # Variables for the BAS package # Note that changes from the previous tutorial! ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) -strt = 10. * np.ones((nlay, nrow, ncol), dtype=np.float32) +strt = 10.0 * np.ones((nlay, nrow, ncol), dtype=np.float32) # Time step parameters nper = 3 @@ -32,20 +31,31 @@ steady = [True, False, False] # Flopy objects -modelname = 'tutorial2' -mf = flopy.modflow.Modflow(modelname, exe_name='mf2005') -dis = flopy.modflow.ModflowDis(mf, nlay, nrow, ncol, delr=delr, delc=delc, - top=ztop, botm=botm[1:], - nper=nper, perlen=perlen, nstp=nstp, - steady=steady) +modelname = "tutorial2" +mf = flopy.modflow.Modflow(modelname, exe_name="mf2005") +dis = flopy.modflow.ModflowDis( + mf, + nlay, + nrow, + ncol, + delr=delr, + delc=delc, + top=ztop, + botm=botm[1:], + nper=nper, + perlen=perlen, + nstp=nstp, + steady=steady, +) bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt) -lpf = flopy.modflow.ModflowLpf(mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp, - ipakcb=53) +lpf = flopy.modflow.ModflowLpf( + mf, hk=hk, vka=vka, sy=sy, ss=ss, laytyp=laytyp, ipakcb=53 +) pcg = flopy.modflow.ModflowPcg(mf) # Make list for stress period 1 -stageleft = 10. -stageright = 10. +stageleft = 10.0 +stageright = 10.0 bound_sp1 = [] for il in range(nlay): condleft = hk * (stageleft - zbot) * delc @@ -53,11 +63,11 @@ for ir in range(nrow): bound_sp1.append([il, ir, 0, stageleft, condleft]) bound_sp1.append([il, ir, ncol - 1, stageright, condright]) -print('Adding ', len(bound_sp1), 'GHBs for stress period 1.') +print("Adding ", len(bound_sp1), "GHBs for stress period 1.") # Make list for stress period 2 -stageleft = 10. -stageright = 0. +stageleft = 10.0 +stageright = 0.0 condleft = hk * (stageleft - zbot) * delc condright = hk * (stageright - zbot) * delc bound_sp2 = [] @@ -65,7 +75,7 @@ for ir in range(nrow): bound_sp2.append([il, ir, 0, stageleft, condleft]) bound_sp2.append([il, ir, ncol - 1, stageright, condright]) -print('Adding ', len(bound_sp2), 'GHBs for stress period 2.') +print("Adding ", len(bound_sp2), "GHBs for stress period 2.") # We do not need to add a dictionary entry for stress period 3. # Flopy will automatically take the list from stress period 2 and apply it @@ -77,10 +87,10 @@ # Create the well package # Remember to use zero-based layer, row, column indices! -pumping_rate = -500. -wel_sp1 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]] -wel_sp2 = [[0, nrow/2 - 1, ncol/2 - 1, 0.]] -wel_sp3 = [[0, nrow/2 - 1, ncol/2 - 1, pumping_rate]] +pumping_rate = -500.0 +wel_sp1 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.0]] +wel_sp2 = [[0, nrow / 2 - 1, ncol / 2 - 1, 0.0]] +wel_sp3 = [[0, nrow / 2 - 1, ncol / 2 - 1, pumping_rate]] stress_period_data = {0: wel_sp1, 1: wel_sp2, 2: wel_sp3} wel = flopy.modflow.ModflowWel(mf, stress_period_data=stress_period_data) @@ -88,13 +98,16 @@ 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 budget'] -oc = flopy.modflow.ModflowOc(mf, stress_period_data=stress_period_data, - compact=True) + stress_period_data[(kper, kstp)] = [ + "save head", + "save drawdown", + "save budget", + "print head", + "print budget", + ] +oc = flopy.modflow.ModflowOc( + mf, stress_period_data=stress_period_data, compact=True +) # Write the model input files mf.write_input() @@ -102,7 +115,7 @@ # Run the model success, mfoutput = mf.run_model(silent=False, pause=False) if not success: - raise Exception('MODFLOW did not terminate normally.') + raise Exception("MODFLOW did not terminate normally.") # Imports @@ -110,67 +123,73 @@ import flopy.utils.binaryfile as bf # Create the headfile and budget file objects -headobj = bf.HeadFile(modelname+'.hds') +headobj = bf.HeadFile(modelname + ".hds") times = headobj.get_times() -cbb = bf.CellBudgetFile(modelname+'.cbc') +cbb = bf.CellBudgetFile(modelname + ".cbc") # Setup contour parameters levels = np.linspace(0, 10, 11) -extent = (delr/2., Lx - delr/2., delc/2., Ly - delc/2.) -print('Levels: ', levels) -print('Extent: ', extent) +extent = (delr / 2.0, Lx - delr / 2.0, delc / 2.0, Ly - delc / 2.0) +print("Levels: ", levels) +print("Extent: ", extent) # Well point -wpt = ((float(ncol/2)-0.5)*delr, (float(nrow/2-1)+0.5)*delc) -wpt = (450., 550.) +wpt = ((float(ncol / 2) - 0.5) * delr, (float(nrow / 2 - 1) + 0.5) * delc) +wpt = (450.0, 550.0) # Make the plots mytimes = [1.0, 101.0, 201.0] for iplot, time in enumerate(mytimes): - print('*****Processing time: ', time) + print("*****Processing time: ", time) head = headobj.get_data(totim=time) - #Print statistics - print('Head statistics') - print(' min: ', head.min()) - print(' max: ', head.max()) - print(' std: ', head.std()) + # Print statistics + print("Head statistics") + print(" min: ", head.min()) + print(" max: ", head.max()) + print(" std: ", head.std()) # Extract flow right face and flow front face - frf = cbb.get_data(text='FLOW RIGHT FACE', totim=time)[0] - fff = cbb.get_data(text='FLOW FRONT FACE', totim=time)[0] - - #Create the plot - #plt.subplot(1, len(mytimes), iplot + 1, aspect='equal') - plt.subplot(1, 1, 1, aspect='equal') - plt.title('stress period ' + str(iplot + 1)) + frf = cbb.get_data(text="FLOW RIGHT FACE", totim=time)[0] + fff = cbb.get_data(text="FLOW FRONT FACE", totim=time)[0] + # Create the plot + # plt.subplot(1, len(mytimes), iplot + 1, aspect='equal') + plt.subplot(1, 1, 1, aspect="equal") + plt.title("stress period " + str(iplot + 1)) modelmap = flopy.plot.ModelMap(model=mf, layer=0) qm = modelmap.plot_ibound() lc = modelmap.plot_grid() - qm = modelmap.plot_bc('GHB', alpha=0.5) + qm = modelmap.plot_bc("GHB", alpha=0.5) cs = modelmap.contour_array(head, levels=levels) - plt.clabel(cs, inline=1, fontsize=10, fmt='%1.1f') + plt.clabel(cs, inline=1, fontsize=10, fmt="%1.1f") quiver = modelmap.plot_discharge(frf, fff, head=head) - - mfc = 'None' - if (iplot+1) == len(mytimes): - mfc='black' - plt.plot(wpt[0], wpt[1], lw=0, marker='o', markersize=8, - markeredgewidth=0.5, - markeredgecolor='black', markerfacecolor=mfc, zorder=9) - plt.text(wpt[0]+25, wpt[1]-25, 'well', size=12, zorder=12) - plt.savefig('tutorial2-{}.png'.format(iplot)) + mfc = "None" + if (iplot + 1) == len(mytimes): + mfc = "black" + plt.plot( + wpt[0], + wpt[1], + lw=0, + marker="o", + markersize=8, + markeredgewidth=0.5, + markeredgecolor="black", + markerfacecolor=mfc, + zorder=9, + ) + plt.text(wpt[0] + 25, wpt[1] - 25, "well", size=12, zorder=12) + plt.savefig("tutorial2-{}.png".format(iplot)) # Plot the head versus time -idx = (0, int(nrow/2) - 1, int(ncol/2) - 1) +idx = (0, int(nrow / 2) - 1, int(ncol / 2) - 1) ts = headobj.get_ts(idx) plt.subplot(1, 1, 1) -ttl = 'Head at cell ({0},{1},{2})'.format(idx[0] + 1, idx[1] + 1, idx[2] + 1) +ttl = "Head at cell ({0},{1},{2})".format(idx[0] + 1, idx[1] + 1, idx[2] + 1) plt.title(ttl) -plt.xlabel('time') -plt.ylabel('head') -plt.plot(ts[:, 0], ts[:, 1], 'bo-') -plt.savefig('tutorial2-ts.png') +plt.xlabel("time") +plt.ylabel("head") +plt.plot(ts[:, 0], ts[:, 1], "bo-") +plt.savefig("tutorial2-ts.png") diff --git a/examples/groundwater_paper/scripts/ex01-script.py b/examples/groundwater_paper/scripts/ex01-script.py index 0a0fb9ade3..29587bd8de 100644 --- a/examples/groundwater_paper/scripts/ex01-script.py +++ b/examples/groundwater_paper/scripts/ex01-script.py @@ -1,4 +1,5 @@ import flopy + model = flopy.modflow.Modflow() dis = flopy.modflow.ModflowDis() bas = flopy.modflow.ModflowBas() diff --git a/examples/groundwater_paper/scripts/uspb_capture.py b/examples/groundwater_paper/scripts/uspb_capture.py index 42ea1a7b5d..046ce380d9 100644 --- a/examples/groundwater_paper/scripts/uspb_capture.py +++ b/examples/groundwater_paper/scripts/uspb_capture.py @@ -5,19 +5,24 @@ import flopy -base_pth = os.path.join('data', 'uspb', 'flopy') -cf_pth = os.path.join('data', 'uspb', 'cf') -res_pth = os.path.join('data', 'uspb', 'results') +base_pth = os.path.join("data", "uspb", "flopy") +cf_pth = os.path.join("data", "uspb", "cf") +res_pth = os.path.join("data", "uspb", "results") + def cf_model(model, k, i, j, base, Q=-100): wd = {1: [[k, i, j, Q]]} - model.remove_package('WEL') + model.remove_package("WEL") wel = flopy.modflow.ModflowWel(model=model, stress_period_data=wd) wel.write_file() model.run_model(silent=True) # get the results - hedObj = flopy.utils.HeadFile(os.path.join(cf_pth, 'DG.hds'), precision='double') - cbcObj = flopy.utils.CellBudgetFile(os.path.join(cf_pth, 'DG.cbc'), precision='double') + hedObj = flopy.utils.HeadFile( + os.path.join(cf_pth, "DG.hds"), precision="double" + ) + cbcObj = flopy.utils.CellBudgetFile( + os.path.join(cf_pth, "DG.cbc"), precision="double" + ) kk = hedObj.get_kstpkper() v = np.zeros((len(kk)), dtype=np.float) h = hedObj.get_ts((k, i, j)) @@ -25,15 +30,22 @@ def cf_model(model, k, i, j, base, Q=-100): if h[idx, 1] == model.lpf.hdry: v[idx] = np.nan else: - v1 = cbcObj.get_data(kstpkper=kon, text='DRAINS', full3D=True)[0] - v2 = cbcObj.get_data(kstpkper=kon, text='STREAM LEAKAGE', full3D=True)[0] - v3 = cbcObj.get_data(kstpkper=kon, text='ET', full3D=True)[0] + v1 = cbcObj.get_data(kstpkper=kon, text="DRAINS", full3D=True)[0] + v2 = cbcObj.get_data( + kstpkper=kon, text="STREAM LEAKAGE", full3D=True + )[0] + v3 = cbcObj.get_data(kstpkper=kon, text="ET", full3D=True)[0] v[idx] = ((v1.sum() + v2.sum() + v3.sum()) - base) / (-Q) return v -ml = flopy.modflow.Modflow.load('DG.nam', version='mf2005', exe_name='mf2005dbl', - verbose=True, model_ws=base_pth) +ml = flopy.modflow.Modflow.load( + "DG.nam", + version="mf2005", + exe_name="mf2005dbl", + verbose=True, + model_ws=base_pth, +) # set a few variables from the model nrow, ncol = ml.dis.nrow, ml.dis.ncol @@ -45,34 +57,38 @@ def cf_model(model, k, i, j, base, Q=-100): ml.run_model() # get base model results -cbcObj = flopy.utils.CellBudgetFile(os.path.join(cf_pth, 'DG.cbc'), precision='double') -v1 = cbcObj.get_data(kstpkper=(0, 0), text='DRAINS', full3D=True)[0] -v2 = cbcObj.get_data(kstpkper=(0, 0), text='STREAM LEAKAGE', full3D=True)[0] -v3 = cbcObj.get_data(kstpkper=(0, 0), text='ET', full3D=True)[0] +cbcObj = flopy.utils.CellBudgetFile( + os.path.join(cf_pth, "DG.cbc"), precision="double" +) +v1 = cbcObj.get_data(kstpkper=(0, 0), text="DRAINS", full3D=True)[0] +v2 = cbcObj.get_data(kstpkper=(0, 0), text="STREAM LEAKAGE", full3D=True)[0] +v3 = cbcObj.get_data(kstpkper=(0, 0), text="ET", full3D=True)[0] baseQ = v1.sum() + v2.sum() + v3.sum() # modify OC -ml.remove_package('OC') -stress_period_data = {(1,9): ['save head', 'save budget', 'print budget'], - (1,10): [], - (1,19): ['save head', 'save budget', 'print budget'], - (1,20): [], - (1,29): ['save head', 'save budget', 'print budget'], - (1,30): [], - (1,39): ['save head', 'save budget', 'print budget'], - (1,40): [], - (1,49): ['save head', 'save budget', 'print budget'], - (1,50): [], - (1,59): ['save head', 'save budget', 'print budget'], - (1,60): [], - (1,69): ['save head', 'save budget', 'print budget'], - (1,70): [], - (1,79): ['save head', 'save budget', 'print budget'], - (1,80): [], - (1,89): ['save head', 'save budget', 'print budget'], - (1,90): [], - (1,99): ['save head', 'save budget', 'print budget'], - (1,100): []} +ml.remove_package("OC") +stress_period_data = { + (1, 9): ["save head", "save budget", "print budget"], + (1, 10): [], + (1, 19): ["save head", "save budget", "print budget"], + (1, 20): [], + (1, 29): ["save head", "save budget", "print budget"], + (1, 30): [], + (1, 39): ["save head", "save budget", "print budget"], + (1, 40): [], + (1, 49): ["save head", "save budget", "print budget"], + (1, 50): [], + (1, 59): ["save head", "save budget", "print budget"], + (1, 60): [], + (1, 69): ["save head", "save budget", "print budget"], + (1, 70): [], + (1, 79): ["save head", "save budget", "print budget"], + (1, 80): [], + (1, 89): ["save head", "save budget", "print budget"], + (1, 90): [], + (1, 99): ["save head", "save budget", "print budget"], + (1, 100): [], +} oc = flopy.modflow.ModflowOc(ml, stress_period_data=stress_period_data) oc.write_file() @@ -82,12 +98,22 @@ def cf_model(model, k, i, j, base, Q=-100): ncol2 = ncol // nstep # open summary file -fs = open(os.path.join('data', 'uspb', 'uspb_capture_{}.out'.format(nstep)), 'w', 0) +fs = open( + os.path.join("data", "uspb", "uspb_capture_{}.out".format(nstep)), "w", 0 +) # write some summary information -fs.write('Problem size: {} rows and {} columns.\n'.format(nrow, ncol)) -fs.write('Capture fraction analysis performed every {} rows and columns.\n'.format(nstep)) -fs.write('Maximum number of analyses: {} rows and {} columns.\n'.format(nrow2, ncol2)) +fs.write("Problem size: {} rows and {} columns.\n".format(nrow, ncol)) +fs.write( + "Capture fraction analysis performed every {} rows and columns.\n".format( + nstep + ) +) +fs.write( + "Maximum number of analyses: {} rows and {} columns.\n".format( + nrow2, ncol2 + ) +) # create array to store capture fraction data (subset of model) cf_array = np.empty((10, nrow2, ncol2), dtype=np.float) @@ -104,15 +130,17 @@ def cf_model(model, k, i, j, base, Q=-100): jcnt = 0 for j in range(0, ncol, nstep): if ibound[i, j] < 1: - sys.stdout.write('.') + sys.stdout.write(".") else: - line = '\nrow {} of {} - col {} of {}\n'.format(icnt+1, nrow2, jcnt+1, ncol2) + line = "\nrow {} of {} - col {} of {}\n".format( + icnt + 1, nrow2, jcnt + 1, ncol2 + ) fs.write(line) sys.stdout.write(line) s0 = time.time() cf = cf_model(ml, 3, i, j, baseQ) s1 = time.time() - line = ' model {} run time: {} seconds\n'.format(idx, s1-s0) + line = " model {} run time: {} seconds\n".format(idx, s1 - s0) fs.write(line) sys.stdout.write(line) idx += 1 @@ -127,14 +155,16 @@ def cf_model(model, k, i, j, base, Q=-100): # end timer for capture fraction analysis end = time.time() ets = end - start -line = '\n' + \ - 'streamflow capture analysis took {} seconds.\n'.format(ets) + \ - 'streamflow capture analysis took {} minutes.\n'.format(ets/60.) + \ - 'streamflow capture analysis took {} hours.\n'.format(ets/3600.) +line = ( + "\n" + + "streamflow capture analysis took {} seconds.\n".format(ets) + + "streamflow capture analysis took {} minutes.\n".format(ets / 60.0) + + "streamflow capture analysis took {} hours.\n".format(ets / 3600.0) +) fs.write(line) sys.stdout.write(line) -#close summary file +# close summary file fs.close() # clean up working directory @@ -146,7 +176,9 @@ def cf_model(model, k, i, j, base, Q=-100): if not os.path.exists(res_pth): os.makedirs(res_pth) for idx in range(10): - fn = os.path.join(res_pth, 'USPB_capture_fraction_{:02d}_{:02d}.dat'.format(nstep, idx+1)) - print('saving capture fraction data to...{}'.format(os.path.basename(fn))) - np.savetxt(fn, cf_array[idx, :, :], delimiter=' ') - + fn = os.path.join( + res_pth, + "USPB_capture_fraction_{:02d}_{:02d}.dat".format(nstep, idx + 1), + ) + print("saving capture fraction data to...{}".format(os.path.basename(fn))) + np.savetxt(fn, cf_array[idx, :, :], delimiter=" ") diff --git a/examples/groundwater_paper/scripts/uspb_capture_par.py b/examples/groundwater_paper/scripts/uspb_capture_par.py index 32c84303f7..5bb06e6879 100644 --- a/examples/groundwater_paper/scripts/uspb_capture_par.py +++ b/examples/groundwater_paper/scripts/uspb_capture_par.py @@ -15,83 +15,102 @@ import flopy # global executable name and output precision -precision = 'double' -exe_name = 'mf2005dbl' -if platform.system() == 'Windows': - exe_name = 'mf2005dbl_x64.exe' +precision = "double" +exe_name = "mf2005dbl" +if platform.system() == "Windows": + exe_name = "mf2005dbl_x64.exe" # functions that do all of the work def load_base_model(klay): # paths - base_pth = os.path.join('data', 'uspb', 'flopy') + base_pth = os.path.join("data", "uspb", "flopy") + + sys.stdout.write("loading base model\n") + ml = flopy.modflow.Modflow.load( + "DG.nam", + version="mf2005", + exe_name=exe_name, + verbose=True, + model_ws=base_pth, + ) - sys.stdout.write('loading base model\n') - ml = flopy.modflow.Modflow.load('DG.nam', version='mf2005', exe_name=exe_name, - verbose=True, model_ws=base_pth) - # set a few variables from the model nrow, ncol = ml.dis.nrow, ml.dis.ncol ibound = ml.bas6.ibound[klay, :, :] - + return ml, nrow, ncol, ibound - - + + def get_baseQ(model): - sys.stdout.write('\nrunning base model to get base head-dependent flow\n\n') + sys.stdout.write( + "\nrunning base model to get base head-dependent flow\n\n" + ) success, report = model.run_model(silent=True, report=True) - sys.stdout.write('Base model run: {}\n'.format(report[-3])) - + sys.stdout.write("Base model run: {}\n".format(report[-3])) + # get base model results - cbcObj = flopy.utils.CellBudgetFile(os.path.join(model.model_ws, 'DG.cbc'), precision=precision) - v1 = cbcObj.get_data(kstpkper=(0, 0), text='DRAINS', full3D=True)[0] - v2 = cbcObj.get_data(kstpkper=(0, 0), text='STREAM LEAKAGE', full3D=True)[0] - v3 = cbcObj.get_data(kstpkper=(0, 0), text='ET', full3D=True)[0] + cbcObj = flopy.utils.CellBudgetFile( + os.path.join(model.model_ws, "DG.cbc"), precision=precision + ) + v1 = cbcObj.get_data(kstpkper=(0, 0), text="DRAINS", full3D=True)[0] + v2 = cbcObj.get_data(kstpkper=(0, 0), text="STREAM LEAKAGE", full3D=True)[ + 0 + ] + v3 = cbcObj.get_data(kstpkper=(0, 0), text="ET", full3D=True)[0] return v1.sum() + v2.sum() + v3.sum() - - -def copy_files(ml, nproc): + + +def copy_files(ml, nproc): # path - cf_base = os.path.join('data', 'uspb') + cf_base = os.path.join("data", "uspb") - exclude = ['hds', 'cbc', 'list', 'ddn'] + exclude = ["hds", "cbc", "list", "ddn"] cf_pths = [] for idx in range(nproc): - cf_pths.append(os.path.join(cf_base, 'cf{:02d}'.format(idx))) + cf_pths.append(os.path.join(cf_base, "cf{:02d}".format(idx))) # create base model in each directory if idx == 0: ml.model_ws = cf_pths[idx] # modify the oc file - ml.remove_package('OC') - stress_period_data = {(1,9): ['save head', 'save budget', 'print budget'], - (1,10): [], - (1,19): ['save head', 'save budget', 'print budget'], - (1,20): [], - (1,29): ['save head', 'save budget', 'print budget'], - (1,30): [], - (1,39): ['save head', 'save budget', 'print budget'], - (1,40): [], - (1,49): ['save head', 'save budget', 'print budget'], - (1,50): [], - (1,59): ['save head', 'save budget', 'print budget'], - (1,60): [], - (1,69): ['save head', 'save budget', 'print budget'], - (1,70): [], - (1,79): ['save head', 'save budget', 'print budget'], - (1,80): [], - (1,89): ['save head', 'save budget', 'print budget'], - (1,90): [], - (1,99): ['save head', 'save budget', 'print budget'], - (1,100): []} - oc = flopy.modflow.ModflowOc(ml, stress_period_data=stress_period_data) + ml.remove_package("OC") + stress_period_data = { + (1, 9): ["save head", "save budget", "print budget"], + (1, 10): [], + (1, 19): ["save head", "save budget", "print budget"], + (1, 20): [], + (1, 29): ["save head", "save budget", "print budget"], + (1, 30): [], + (1, 39): ["save head", "save budget", "print budget"], + (1, 40): [], + (1, 49): ["save head", "save budget", "print budget"], + (1, 50): [], + (1, 59): ["save head", "save budget", "print budget"], + (1, 60): [], + (1, 69): ["save head", "save budget", "print budget"], + (1, 70): [], + (1, 79): ["save head", "save budget", "print budget"], + (1, 80): [], + (1, 89): ["save head", "save budget", "print budget"], + (1, 90): [], + (1, 99): ["save head", "save budget", "print budget"], + (1, 100): [], + } + oc = flopy.modflow.ModflowOc( + ml, stress_period_data=stress_period_data + ) # write the input files ml.write_input() else: if not os.path.exists(cf_pths[idx]): os.makedirs(cf_pths[idx]) filelist = [f for f in os.listdir(cf_pths[0])] - sys.stdout.write('copying files from {} to {}\n'.format(cf_pths[0], cf_pths[idx])) + sys.stdout.write( + "copying files from {} to {}\n".format( + cf_pths[0], cf_pths[idx] + ) + ) for f in filelist: if os.path.splitext(f)[1].lower() in exclude: continue @@ -100,46 +119,48 @@ def copy_files(ml, nproc): shutil.copyfile(src, dst) return ml, cf_pths - + + # functions to run the models in parallel def unpack_args(args): try: f, idx, nmax, k, i, j, Qt, base, hdry = args except: - sys.stdout.write('could not unpack args\n') + sys.stdout.write("could not unpack args\n") raise try: current = mp.current_process() - imod = current._identity[0]-1 + imod = current._identity[0] - 1 except: - sys.stdout.write('could not get current process\n') + sys.stdout.write("could not get current process\n") raise return f(imod, idx, nmax, k, i, j, Qt, base, hdry) + def make_well(pth, k, i, j, Qt): - fn = os.path.join(pth, 'DG.wel') - f = open(fn, 'w') - f.write('# Well file for MODFLOW, generated by Flopy.\n') - f.write(' 1 0\n') - f.write(' 0 0 # stress period 0\n') - f.write(' 1 0 # stress period 1\n') - f.write('{:10d}{:10d}{:10d}{:10.2f}\n'.format(k+1, i+1, j+1, Qt)) + fn = os.path.join(pth, "DG.wel") + f = open(fn, "w") + f.write("# Well file for MODFLOW, generated by Flopy.\n") + f.write(" 1 0\n") + f.write(" 0 0 # stress period 0\n") + f.write(" 1 0 # stress period 1\n") + f.write("{:10d}{:10d}{:10d}{:10.2f}\n".format(k + 1, i + 1, j + 1, Qt)) f.close() + def run_model(pth): - proc = sp.Popen([exe_name, 'DG.nam'], - stdout=sp.PIPE, cwd=pth) - sys.stdout.write(' running {} in {}\n'.format(exe_name, pth)) + proc = sp.Popen([exe_name, "DG.nam"], stdout=sp.PIPE, cwd=pth) + sys.stdout.write(" running {} in {}\n".format(exe_name, pth)) success = False buff = [] - elt = 'Normal model termination did not occur' + elt = "Normal model termination did not occur" while True: line = proc.stdout.readline() - c = line.decode('utf-8') - if c != '': - if 'normal termination of simulation' in c.lower(): + c = line.decode("utf-8") + if c != "": + if "normal termination of simulation" in c.lower(): success = True - c = c.rstrip('\r\n') + c = c.rstrip("\r\n") buff.append(c) else: break @@ -148,96 +169,136 @@ def run_model(pth): elt = buff[-3].strip() except: pass - return ([success, elt]) + return [success, elt] + # function to create well file, run model, and extract results def cf_model(imod, ion, nmax, k, i, j, Qt, base, hdry): - pth = os.path.join('data', 'uspb', 'cf{:02d}'.format(imod)) - sys.stdout.write('\nRunning model number: {}\n'.format(imod)) - sys.stdout.write(' model run: {} of {}\n'.format(ion+1, nmax)) - sys.stdout.write(' model number {} working directory: {}\n'.format(imod, pth)) + pth = os.path.join("data", "uspb", "cf{:02d}".format(imod)) + sys.stdout.write("\nRunning model number: {}\n".format(imod)) + sys.stdout.write(" model run: {} of {}\n".format(ion + 1, nmax)) + sys.stdout.write( + " model number {} working directory: {}\n".format(imod, pth) + ) make_well(pth, k, i, j, Qt) success, elt = run_model(pth) - line = '\nModel run: {} of {} (model number {})\n'.format(ion+1, nmax, imod) - line += ' row {} - col {}\n'.format(i+1, j+1) - line += ' {}\n'.format(elt) + line = "\nModel run: {} of {} (model number {})\n".format( + ion + 1, nmax, imod + ) + line += " row {} - col {}\n".format(i + 1, j + 1) + line += " {}\n".format(elt) # get the results v = np.zeros((10), dtype=np.float) if success: try: - hedObj = flopy.utils.HeadFile(os.path.join(pth, 'DG.hds'), precision=precision) - cbcObj = flopy.utils.CellBudgetFile(os.path.join(pth, 'DG.cbc'), precision=precision) + hedObj = flopy.utils.HeadFile( + os.path.join(pth, "DG.hds"), precision=precision + ) + cbcObj = flopy.utils.CellBudgetFile( + os.path.join(pth, "DG.cbc"), precision=precision + ) kk = hedObj.get_kstpkper() h = hedObj.get_ts((k, i, j)) for idx, kon in enumerate(kk): if h[idx, 1] == hdry: v[idx] = np.nan else: - v1 = cbcObj.get_data(kstpkper=kon, text='DRAINS', full3D=True)[0] - v2 = cbcObj.get_data(kstpkper=kon, text='STREAM LEAKAGE', full3D=True)[0] - v3 = cbcObj.get_data(kstpkper=kon, text='ET', full3D=True)[0] + v1 = cbcObj.get_data( + kstpkper=kon, text="DRAINS", full3D=True + )[0] + v2 = cbcObj.get_data( + kstpkper=kon, text="STREAM LEAKAGE", full3D=True + )[0] + v3 = cbcObj.get_data(kstpkper=kon, text="ET", full3D=True)[ + 0 + ] v[idx] = ((v1.sum() + v2.sum() + v3.sum()) - base) / (-Qt) except: - line += ' Error: Model run: {} of {} (model number {}) - '.format(ion+1, nmax, imod) - line += 'could not process model results.\n' + line += " Error: Model run: {} of {} (model number {}) - ".format( + ion + 1, nmax, imod + ) + line += "could not process model results.\n" v[:] = np.nan else: - line += ' Error: Model run: {} of {} (model number {}) '.format(ion+1, nmax, imod) - line += 'did not execute successfully\n' + line += " Error: Model run: {} of {} (model number {}) ".format( + ion + 1, nmax, imod + ) + line += "did not execute successfully\n" v[:] = np.nan sys.stdout.write(line) return (v, line) + def doit(): # multi processing information nproc = 3 ncores = mp.cpu_count() if nproc > ncores: - sys.stdout.write('Requested {} cores but only {} cores are available.\n\n\n'.format(nproc, ncores)) + sys.stdout.write( + "Requested {} cores but only {} cores are available.\n\n\n".format( + nproc, ncores + ) + ) else: - sys.stdout.write('Requested {} cores and {} cores are available.\n\n\n'.format(nproc, ncores)) - - #paths - res_pth = os.path.join('data', 'uspb', 'results') - - #model data + sys.stdout.write( + "Requested {} cores and {} cores are available.\n\n\n".format( + nproc, ncores + ) + ) + + # paths + res_pth = os.path.join("data", "uspb", "results") + + # model data klay = 3 - Qcf = -100. + Qcf = -100.0 nstep = 4 # load base model ml, nrow, ncol, ibound = load_base_model(klay) - + # run first model created to get base model results baseQ = get_baseQ(ml) - sys.stdout.write('Base head-dependent flux = {}'.format(baseQ)) - + sys.stdout.write("Base head-dependent flux = {}".format(baseQ)) + # modify oc file copy model files ml, cf_pths = copy_files(ml, nproc) - + # calculate subset of model to run nrow2 = nrow // nstep ncol2 = ncol // nstep # open summary file - fs = open(os.path.join('data', 'uspb', 'uspb_capture_{}.out'.format(nstep)), 'w', 0) - + fs = open( + os.path.join("data", "uspb", "uspb_capture_{}.out".format(nstep)), + "w", + 0, + ) + # write some summary information - fs.write('Problem size: {} rows and {} columns.\n'.format(nrow, ncol)) - fs.write('Capture fraction analysis performed every {} rows and columns.\n'.format(nstep)) - fs.write('Maximum number of analyses: {} rows and {} columns.\n'.format(nrow2, ncol2)) - + fs.write("Problem size: {} rows and {} columns.\n".format(nrow, ncol)) + fs.write( + "Capture fraction analysis performed every {} rows and columns.\n".format( + nstep + ) + ) + fs.write( + "Maximum number of analyses: {} rows and {} columns.\n".format( + nrow2, ncol2 + ) + ) + # create array to store capture fraction data (subset of model) cf_array = np.empty((10, nrow2, ncol2), dtype=np.float) cf_array.fill(np.nan) - + # timer for capture fraction analysis start = time.time() - + # capture fraction analysis icnt = 0 jcnt = 0 - + # build tuple with list of cells cells = [] cellmap = [] @@ -252,57 +313,67 @@ def doit(): jcnt += 1 # increment icnt icnt += 1 - + ## test cg_model function - #t = cf_model(models[0], klay, cells[0][0], cells[0][1], Qcf, baseQ) - #sys.stdout.write(t) - + # t = cf_model(models[0], klay, cells[0][0], cells[0][1], Qcf, baseQ) + # sys.stdout.write(t) + # create multiprocessing pool pool = mp.Pool(processes=nproc) - args = [(cf_model, idx, len(cells), klay, i, j, Qcf, baseQ, ml.lpf.hdry) for idx, (i, j) in enumerate(cells)] - #sys.stdout.write(args) + args = [ + (cf_model, idx, len(cells), klay, i, j, Qcf, baseQ, ml.lpf.hdry) + for idx, (i, j) in enumerate(cells) + ] + # sys.stdout.write(args) output = pool.map(unpack_args, args, nproc) pool.close() pool.join() - + for v in output: fs.write(v[1]) - + for idx, (icnt, jcnt) in enumerate(cellmap): # add values to the array if icnt < nrow2 and jcnt < ncol2: cf_array[:, icnt, jcnt] = output[idx][0].copy() - - + # end timer for capture fraction analysis end = time.time() ets = end - start - line = '\n' + \ - 'streamflow capture analysis took {} seconds.\n'.format(ets) + \ - 'streamflow capture analysis took {} minutes.\n'.format(ets/60.) + \ - 'streamflow capture analysis took {} hours.\n'.format(ets/3600.) + line = ( + "\n" + + "streamflow capture analysis took {} seconds.\n".format(ets) + + "streamflow capture analysis took {} minutes.\n".format(ets / 60.0) + + "streamflow capture analysis took {} hours.\n".format(ets / 3600.0) + ) fs.write(line) sys.stdout.write(line) - - #close summary file + + # close summary file fs.close() - + # clean up working directories for idx in range(nproc): filelist = [f for f in os.listdir(cf_pths[idx])] for f in filelist: os.remove(os.path.join(cf_pths[idx], f)) - + # create res_pth (if it doesn't exist) and save data if not os.path.exists(res_pth): os.makedirs(res_pth) for idx in range(10): - fn = os.path.join(res_pth, 'USPB_capture_fraction_{:02d}_{:02d}.dat'.format(nstep, idx+1)) - sys.stdout.write('saving capture fraction data to...{}\n'.format(os.path.basename(fn))) - np.savetxt(fn, cf_array[idx, :, :], delimiter=' ') - + fn = os.path.join( + res_pth, + "USPB_capture_fraction_{:02d}_{:02d}.dat".format(nstep, idx + 1), + ) + sys.stdout.write( + "saving capture fraction data to...{}\n".format( + os.path.basename(fn) + ) + ) + np.savetxt(fn, cf_array[idx, :, :], delimiter=" ") -if __name__ == '__main__': - - doit() +if __name__ == "__main__": + + doit() diff --git a/examples/scripts/flopy_henry.py b/examples/scripts/flopy_henry.py index 3b3aa48e0f..ab1ef65aac 100644 --- a/examples/scripts/flopy_henry.py +++ b/examples/scripts/flopy_henry.py @@ -6,44 +6,55 @@ def run(): - workspace = os.path.join('henry') + workspace = os.path.join("henry") # make sure workspace directory exists if not os.path.exists(workspace): os.makedirs(workspace) - fext = 'png' + fext = "png" narg = len(sys.argv) iarg = 0 if narg > 1: while iarg < narg - 1: iarg += 1 basearg = sys.argv[iarg].lower() - if basearg == '--pdf': - fext = 'pdf' + if basearg == "--pdf": + fext = "pdf" # Input variables for the Henry Problem - Lx = 2. - Lz = 1. + Lx = 2.0 + Lz = 1.0 nlay = 50 nrow = 1 ncol = 100 delr = Lx / ncol delc = 1.0 delv = Lz / nlay - henry_top = 1. - henry_botm = np.linspace(henry_top - delv, 0., nlay) + henry_top = 1.0 + henry_botm = np.linspace(henry_top - delv, 0.0, nlay) qinflow = 5.702 # m3/day dmcoef = 0.57024 # m2/day Could also try 1.62925 as another case of the Henry problem - hk = 864. # m/day + hk = 864.0 # m/day # Create the basic MODFLOW model data - modelname = 'henry' + modelname = "henry" m = flopy.seawat.Seawat(modelname, exe_name="swtv4", model_ws=workspace) # Add DIS package to the MODFLOW model - dis = flopy.modflow.ModflowDis(m, nlay, nrow, ncol, nper=1, delr=delr, - delc=delc, laycbd=0, top=henry_top, - botm=henry_botm, perlen=1.5, nstp=15) + dis = flopy.modflow.ModflowDis( + m, + nlay, + nrow, + ncol, + nper=1, + delr=delr, + delc=delc, + laycbd=0, + top=henry_top, + botm=henry_botm, + perlen=1.5, + nstp=15, + ) # Variables for the BAS package ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) @@ -54,13 +65,14 @@ def run(): lpf = flopy.modflow.ModflowLpf(m, hk=hk, vka=hk, ipakcb=53) # Add PCG Package to the MODFLOW model - pcg = flopy.modflow.ModflowPcg(m, hclose=1.e-8) + pcg = flopy.modflow.ModflowPcg(m, hclose=1.0e-8) # Add OC package to the MODFLOW model - oc = flopy.modflow.ModflowOc(m, - stress_period_data={ - (0, 0): ['save head', 'save budget']}, - compact=True) + oc = flopy.modflow.ModflowOc( + m, + stress_period_data={(0, 0): ["save head", "save budget"]}, + compact=True, + ) # Create WEL and SSM data itype = flopy.mt3d.Mt3dSsm.itype_dict() @@ -70,32 +82,48 @@ def run(): ssm_sp1 = [] for k in range(nlay): wel_sp1.append([k, 0, 0, qinflow / nlay]) - ssm_sp1.append([k, 0, 0, 0., itype['WEL']]) - ssm_sp1.append([k, 0, ncol - 1, 35., itype['BAS6']]) + ssm_sp1.append([k, 0, 0, 0.0, itype["WEL"]]) + ssm_sp1.append([k, 0, ncol - 1, 35.0, itype["BAS6"]]) wel_data[0] = wel_sp1 ssm_data[0] = ssm_sp1 wel = flopy.modflow.ModflowWel(m, stress_period_data=wel_data) # Create the basic MT3DMS model data - btn = flopy.mt3d.Mt3dBtn(m, nprs=-5, prsity=0.35, sconc=35., ifmtcn=0, - chkmas=False, nprobs=10, nprmas=10, dt0=0.001) + btn = flopy.mt3d.Mt3dBtn( + m, + nprs=-5, + prsity=0.35, + sconc=35.0, + ifmtcn=0, + chkmas=False, + nprobs=10, + nprmas=10, + dt0=0.001, + ) adv = flopy.mt3d.Mt3dAdv(m, mixelm=0) - dsp = flopy.mt3d.Mt3dDsp(m, al=0., trpt=1., trpv=1., dmcoef=dmcoef) + dsp = flopy.mt3d.Mt3dDsp(m, al=0.0, trpt=1.0, trpv=1.0, dmcoef=dmcoef) gcg = flopy.mt3d.Mt3dGcg(m, iter1=500, mxiter=1, isolve=1, cclose=1e-7) ssm = flopy.mt3d.Mt3dSsm(m, stress_period_data=ssm_data) # Create the SEAWAT model data - vdf = flopy.seawat.SeawatVdf(m, iwtable=0, densemin=0, densemax=0, - denseref=1000., denseslp=0.7143, firstdt=1e-3) + vdf = flopy.seawat.SeawatVdf( + m, + iwtable=0, + densemin=0, + densemax=0, + denseref=1000.0, + denseslp=0.7143, + firstdt=1e-3, + ) # Write the input files m.write_input() # Try to delete the output files, to prevent accidental use of older files try: - os.remove(os.path.join(workspace, 'MT3D001.UCN')) - os.remove(os.path.join(workspace, modelname + '.hds')) - os.remove(os.path.join(workspace, modelname + '.cbc')) + os.remove(os.path.join(workspace, "MT3D001.UCN")) + os.remove(os.path.join(workspace, modelname + ".hds")) + os.remove(os.path.join(workspace, modelname + ".cbc")) except: pass @@ -106,61 +134,72 @@ def run(): # Load data ucnobj = flopy.utils.binaryfile.UcnFile( - os.path.join(workspace, 'MT3D001.UCN'), - model=m) + os.path.join(workspace, "MT3D001.UCN"), model=m + ) times = ucnobj.get_times() concentration = ucnobj.get_data(totim=times[-1]) cbbobj = flopy.utils.binaryfile.CellBudgetFile( - os.path.join(workspace, 'henry.cbc')) + os.path.join(workspace, "henry.cbc") + ) times = cbbobj.get_times() - qx = cbbobj.get_data(text='flow right face', totim=times[-1])[0] - qz = cbbobj.get_data(text='flow lower face', totim=times[-1])[0] + qx = cbbobj.get_data(text="flow right face", totim=times[-1])[0] + qz = cbbobj.get_data(text="flow lower face", totim=times[-1])[0] # Average flows to cell centers qx_avg = np.empty(qx.shape, dtype=qx.dtype) - qx_avg[:, :, 1:] = 0.5 * (qx[:, :, 0:ncol - 1] + qx[:, :, 1:ncol]) + qx_avg[:, :, 1:] = 0.5 * (qx[:, :, 0 : ncol - 1] + qx[:, :, 1:ncol]) qx_avg[:, :, 0] = 0.5 * qx[:, :, 0] qz_avg = np.empty(qz.shape, dtype=qz.dtype) - qz_avg[1:, :, :] = 0.5 * (qz[0:nlay - 1, :, :] + qz[1:nlay, :, :]) + qz_avg[1:, :, :] = 0.5 * (qz[0 : nlay - 1, :, :] + qz[1:nlay, :, :]) qz_avg[0, :, :] = 0.5 * qz[0, :, :] # Make the plot # import matplotlib.pyplot as plt fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot(1, 1, 1, aspect='equal') - ax.imshow(concentration[:, 0, :], interpolation='nearest', - extent=(0, Lx, 0, Lz)) + ax = fig.add_subplot(1, 1, 1, aspect="equal") + ax.imshow( + concentration[:, 0, :], interpolation="nearest", extent=(0, Lx, 0, Lz) + ) y, x, z = dis.get_node_coordinates() X, Z = np.meshgrid(x, z[:, 0, 0]) iskip = 3 - ax.quiver(X[::iskip, ::iskip], Z[::iskip, ::iskip], - qx_avg[::iskip, 0, ::iskip], -qz_avg[::iskip, 0, ::iskip], - color='w', scale=5, headwidth=3, headlength=2, - headaxislength=2, width=0.0025) - - outfig = os.path.join(workspace, 'henry_flows.{0}'.format(fext)) + ax.quiver( + X[::iskip, ::iskip], + Z[::iskip, ::iskip], + qx_avg[::iskip, 0, ::iskip], + -qz_avg[::iskip, 0, ::iskip], + color="w", + scale=5, + headwidth=3, + headlength=2, + headaxislength=2, + width=0.0025, + ) + + outfig = os.path.join(workspace, "henry_flows.{0}".format(fext)) fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) # Extract the heads - fname = os.path.join(workspace, 'henry.hds') + fname = os.path.join(workspace, "henry.hds") headobj = flopy.utils.binaryfile.HeadFile(fname) times = headobj.get_times() head = headobj.get_data(totim=times[-1]) # Make a simple head plot fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot(1, 1, 1, aspect='equal') - im = ax.imshow(head[:, 0, :], interpolation='nearest', - extent=(0, Lx, 0, Lz)) - ax.set_title('Simulated Heads') + ax = fig.add_subplot(1, 1, 1, aspect="equal") + im = ax.imshow( + head[:, 0, :], interpolation="nearest", extent=(0, Lx, 0, Lz) + ) + ax.set_title("Simulated Heads") - outfig = os.path.join(workspace, 'henry_heads.{0}'.format(fext)) + outfig = os.path.join(workspace, "henry_heads.{0}".format(fext)) fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) return 0 -if __name__ == '__main__': +if __name__ == "__main__": success = run() diff --git a/examples/scripts/flopy_lake_example.py b/examples/scripts/flopy_lake_example.py index 0defd5be02..cb6720aadf 100644 --- a/examples/scripts/flopy_lake_example.py +++ b/examples/scripts/flopy_lake_example.py @@ -7,20 +7,20 @@ def run(): - workspace = os.path.join('lake') + workspace = os.path.join("lake") # make sure workspace directory exists if not os.path.exists(workspace): os.makedirs(workspace) - fext = 'png' + fext = "png" narg = len(sys.argv) iarg = 0 if narg > 1: while iarg < narg - 1: iarg += 1 basearg = sys.argv[iarg].lower() - if basearg == '--pdf': - fext = 'pdf' + if basearg == "--pdf": + fext = "pdf" # save the starting path cwdpth = os.getcwd() @@ -33,7 +33,7 @@ def run(): # of the model and the parameters of the model: the number of layers `Nlay`, the number of rows # and columns `N`, lengths of the sides of the model `L`, aquifer thickness `H`, hydraulic # conductivity `Kh` - name = 'lake_example' + name = "lake_example" h1 = 100 h2 = 90 Nlay = 10 @@ -46,8 +46,9 @@ def run(): # whatever you want). The modelname will be the name given to all MODFLOW files (input and output). # The exe_name should be the full path to your MODFLOW executable. The version is either 'mf2k' # for MODFLOW2000 or 'mf2005'for MODFLOW2005. - ml = flopy.modflow.Modflow(modelname=name, exe_name='mf2005', - version='mf2005') + ml = flopy.modflow.Modflow( + modelname=name, exe_name="mf2005", version="mf2005" + ) # Define the discretization of the model. All layers are given equal thickness. The `bot` array # is build from the `Hlay` values to indicate top and bottom of each layer, and `delrow` and @@ -55,8 +56,17 @@ def run(): # the Discretization file is built. bot = np.linspace(-H / Nlay, -H, Nlay) delrow = delcol = L / (N - 1) - dis = flopy.modflow.ModflowDis(ml, nlay=Nlay, nrow=N, ncol=N, delr=delrow, - delc=delcol, top=0.0, botm=bot, laycbd=0) + dis = flopy.modflow.ModflowDis( + ml, + nlay=Nlay, + nrow=N, + ncol=N, + delr=delrow, + delc=delcol, + top=0.0, + botm=bot, + laycbd=0, + ) # Next we specify the boundary conditions and starting heads with the Basic package. The `ibound` # array will be `1` in all cells in all layers, except for along the boundary and in the cell at @@ -77,14 +87,14 @@ def run(): # create external ibound array and starting head files files = [] - hfile = '{}_strt.ref'.format(name) + hfile = "{}_strt.ref".format(name) np.savetxt(hfile, start) hfiles = [] for kdx in range(Nlay): - file = '{}_ib{:02d}.ref'.format(name, kdx + 1) + file = "{}_ib{:02d}.ref".format(name, kdx + 1) files.append(file) hfiles.append(hfile) - np.savetxt(file, ibound[kdx, :, :], fmt='%5d') + np.savetxt(file, ibound[kdx, :, :], fmt="%5d") bas = flopy.modflow.ModflowBas(ml, ibound=files, strt=hfiles) @@ -108,39 +118,39 @@ def run(): # specifying, in this case, the step number and period number for which we want to retrieve data. # A three-dimensional array is returned of size `nlay, nrow, ncol`. Matplotlib contouring functions # are used to make contours of the layers or a cross-section. - hds = flopy.utils.HeadFile(os.path.join(workspace, name + '.hds')) + hds = flopy.utils.HeadFile(os.path.join(workspace, name + ".hds")) h = hds.get_data(kstpkper=(0, 0)) x = y = np.linspace(0, L, N) c = plt.contour(x, y, h[0], np.arange(90, 100.1, 0.2)) - plt.clabel(c, fmt='%2.1f') - plt.axis('scaled') + plt.clabel(c, fmt="%2.1f") + plt.axis("scaled") - outfig = os.path.join(workspace, 'lake1.{0}'.format(fext)) + outfig = os.path.join(workspace, "lake1.{0}".format(fext)) fig = plt.gcf() fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) x = y = np.linspace(0, L, N) c = plt.contour(x, y, h[-1], np.arange(90, 100.1, 0.2)) - plt.clabel(c, fmt='%1.1f') - plt.axis('scaled') + plt.clabel(c, fmt="%1.1f") + plt.axis("scaled") - outfig = os.path.join(workspace, 'lake2.{0}'.format(fext)) + outfig = os.path.join(workspace, "lake2.{0}".format(fext)) fig = plt.gcf() fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) z = np.linspace(-H / Nlay / 2, -H + H / Nlay / 2, Nlay) - c = plt.contour(x, z, h[:, 50, :], np.arange(90, 100.1, .2)) - plt.axis('scaled') + c = plt.contour(x, z, h[:, 50, :], np.arange(90, 100.1, 0.2)) + plt.axis("scaled") - outfig = os.path.join(workspace, 'lake3.{0}'.format(fext)) + outfig = os.path.join(workspace, "lake3.{0}".format(fext)) fig = plt.gcf() fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) return 0 -if __name__ == '__main__': +if __name__ == "__main__": success = run() diff --git a/examples/scripts/flopy_swi2_ex1.py b/examples/scripts/flopy_swi2_ex1.py index 170e8d29c7..5fb36d003a 100755 --- a/examples/scripts/flopy_swi2_ex1.py +++ b/examples/scripts/flopy_swi2_ex1.py @@ -11,54 +11,56 @@ import matplotlib.pyplot as plt # --modify default matplotlib settings -updates = {'font.family': ['Univers 57 Condensed', 'Arial'], - 'mathtext.default': 'regular', - 'pdf.compression': 0, - 'pdf.fonttype': 42, - 'legend.fontsize': 7, - 'axes.labelsize': 8, - 'xtick.labelsize': 7, - 'ytick.labelsize': 7} +updates = { + "font.family": ["Univers 57 Condensed", "Arial"], + "mathtext.default": "regular", + "pdf.compression": 0, + "pdf.fonttype": 42, + "legend.fontsize": 7, + "axes.labelsize": 8, + "xtick.labelsize": 7, + "ytick.labelsize": 7, +} plt.rcParams.update(updates) def run(): - workspace = 'swiex1' + workspace = "swiex1" cleanFiles = False - fext = 'png' + fext = "png" narg = len(sys.argv) iarg = 0 if narg > 1: while iarg < narg - 1: iarg += 1 basearg = sys.argv[iarg].lower() - if basearg == '--clean': + if basearg == "--clean": cleanFiles = True - elif basearg == '--pdf': - fext = 'pdf' + elif basearg == "--pdf": + fext = "pdf" if cleanFiles: - print('cleaning all files') - print('excluding *.py files') + print("cleaning all files") + print("excluding *.py files") files = os.listdir(workspace) for f in files: if os.path.isdir(f): continue - if '.py' != os.path.splitext(f)[1].lower(): - print(' removing...{}'.format(os.path.basename(f))) + if ".py" != os.path.splitext(f)[1].lower(): + print(" removing...{}".format(os.path.basename(f))) os.remove(os.path.join(workspace, f)) return 1 - modelname = 'swiex1' - exe_name = 'mf2005' + modelname = "swiex1" + exe_name = "mf2005" nlay = 1 nrow = 1 ncol = 50 - delr = 5. - delc = 1. + delr = 5.0 + delc = 1.0 ibound = np.ones((nrow, ncol), np.int) ibound[0, -1] = -1 @@ -75,26 +77,45 @@ def run(): ocdict = {} for idx in range(49, 200, 50): key = (0, idx) - ocdict[key] = ['save head', 'save budget'] + ocdict[key] = ["save head", "save budget"] key = (0, idx + 1) ocdict[key] = [] # create flopy modflow object - ml = flopy.modflow.Modflow(modelname, version='mf2005', exe_name=exe_name, - model_ws=workspace) + ml = flopy.modflow.Modflow( + modelname, version="mf2005", exe_name=exe_name, model_ws=workspace + ) # create flopy modflow package objects - discret = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol, - delr=delr, delc=delc, - top=0, botm=[-40.0], - perlen=400, nstp=200) + discret = flopy.modflow.ModflowDis( + ml, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=0, + botm=[-40.0], + perlen=400, + nstp=200, + ) bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=0.0) - lpf = flopy.modflow.ModflowLpf(ml, hk=2., vka=2.0, vkcb=0, laytyp=0, - layavg=0) + lpf = flopy.modflow.ModflowLpf( + ml, hk=2.0, vka=2.0, vkcb=0, laytyp=0, layavg=0 + ) wel = flopy.modflow.ModflowWel(ml, stress_period_data={0: [(0, 0, 0, 1)]}) - swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, npln=1, istrat=1, - toeslope=0.2, tipslope=0.2, nu=[0, 0.025], - zeta=z, ssz=0.2, isource=isource, - nsolver=1) + swi = flopy.modflow.ModflowSwi2( + ml, + iswizt=55, + npln=1, + istrat=1, + toeslope=0.2, + tipslope=0.2, + nu=[0, 0.025], + zeta=z, + ssz=0.2, + isource=isource, + nsolver=1, + ) oc = flopy.modflow.ModflowOc(ml, stress_period_data=ocdict) pcg = flopy.modflow.ModflowPcg(ml) # create model files @@ -102,17 +123,17 @@ def run(): # run the model m = ml.run_model(silent=False) # read model heads - headfile = os.path.join(workspace, '{}.hds'.format(modelname)) + headfile = os.path.join(workspace, "{}.hds".format(modelname)) hdobj = flopy.utils.HeadFile(headfile) head = hdobj.get_alldata() head = np.array(head) # read model zeta - zetafile = os.path.join(workspace, '{}.zta'.format(modelname)) + zetafile = os.path.join(workspace, "{}.zta".format(modelname)) zobj = flopy.utils.CellBudgetFile(zetafile) zkstpkper = zobj.get_kstpkper() zeta = [] for kk in zkstpkper: - zeta.append(zobj.get_data(kstpkper=kk, text=' ZETASRF 1')[0]) + zeta.append(zobj.get_data(kstpkper=kk, text=" ZETASRF 1")[0]) zeta = np.array(zeta) x = np.arange(0.5 * delr, ncol * delr, delr) @@ -134,70 +155,116 @@ def run(): fbot = 0.125 ftop = 0.925 - fig = plt.figure(figsize=(fwid, fhgt), facecolor='w') - fig.subplots_adjust(wspace=0.25, hspace=0.25, left=flft, right=frgt, - bottom=fbot, top=ftop) + fig = plt.figure(figsize=(fwid, fhgt), facecolor="w") + fig.subplots_adjust( + wspace=0.25, hspace=0.25, left=flft, right=frgt, bottom=fbot, top=ftop + ) ax = fig.add_subplot(211) - ax.text(-0.075, 1.05, 'A', transform=ax.transAxes, va='center', - ha='center', - size='8') - ax.plot([80, 120], [0, -40], 'k') + ax.text( + -0.075, + 1.05, + "A", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) + ax.plot([80, 120], [0, -40], "k") ax.set_xlim(0, 250) ax.set_ylim(-40, 0) ax.set_yticks(np.arange(-40, 1, 10)) - ax.text(50, -10, 'salt') - ax.text(130, -10, 'fresh') - a = ax.annotate("", xy=(50, -25), xytext=(30, -25), - arrowprops=dict(arrowstyle='->', fc='k')) - ax.text(40, -22, 'groundwater flow velocity=0.125 m/d', ha='center', - size=7) - ax.set_ylabel('Elevation, in meters') + ax.text(50, -10, "salt") + ax.text(130, -10, "fresh") + a = ax.annotate( + "", + xy=(50, -25), + xytext=(30, -25), + arrowprops=dict(arrowstyle="->", fc="k"), + ) + ax.text( + 40, -22, "groundwater flow velocity=0.125 m/d", ha="center", size=7 + ) + ax.set_ylabel("Elevation, in meters") ax = fig.add_subplot(212) - ax.text(-0.075, 1.05, 'B', transform=ax.transAxes, va='center', - ha='center', - size='8') + ax.text( + -0.075, + 1.05, + "B", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) for i in range(4): Ltoe[i] = H * math.sqrt(k * nu * (t[i] + tzero) / n / H) - ax.plot([100 - Ltoe[i] + v * t[i], 100 + Ltoe[i] + v * t[i]], [0, -40], - 'k', label='_None') + ax.plot( + [100 - Ltoe[i] + v * t[i], 100 + Ltoe[i] + v * t[i]], + [0, -40], + "k", + label="_None", + ) for i in range(4): zi = zeta[i, 0, 0, :] p = (zi < 0) & (zi > -39.9) - ax.plot(x[p], zeta[i, 0, 0, p], 'bo', - markersize=3, markeredgecolor='blue', markerfacecolor='None', - label='_None') + ax.plot( + x[p], + zeta[i, 0, 0, p], + "bo", + markersize=3, + markeredgecolor="blue", + markerfacecolor="None", + label="_None", + ) ipos = 0 for jdx, t in enumerate(zeta[i, 0, 0, :]): if t > -39.9: ipos = jdx - ax.text(x[ipos], -37.75, '{0} days'.format(((i + 1) * 100)), size=5, - ha='left', va='center') + ax.text( + x[ipos], + -37.75, + "{0} days".format(((i + 1) * 100)), + size=5, + ha="left", + va="center", + ) # fake items for labels - ax.plot([-100., -100], [-100., -100], 'k', label='Analytical solution') - ax.plot([-100., -100], [-100., -100], 'bo', markersize=3, - markeredgecolor='blue', markerfacecolor='None', label='SWI2') + ax.plot([-100.0, -100], [-100.0, -100], "k", label="Analytical solution") + ax.plot( + [-100.0, -100], + [-100.0, -100], + "bo", + markersize=3, + markeredgecolor="blue", + markerfacecolor="None", + label="SWI2", + ) # legend - leg = ax.legend(loc='upper right', numpoints=1) + leg = ax.legend(loc="upper right", numpoints=1) leg._drawFrame = False # axes ax.set_xlim(0, 250) ax.set_ylim(-40, 0) ax.set_yticks(np.arange(-40, 1, 10)) - a = ax.annotate("", xy=(50, -25), xytext=(30, -25), - arrowprops=dict(arrowstyle='->', fc='k')) - ax.text(40, -22, 'groundwater flow velocity=0.125 m/d', ha='center', - size=7) - ax.set_ylabel('Elevation, in meters') - ax.set_xlabel('Horizontal distance, in meters') - - outfig = os.path.join(workspace, 'Figure06_swi2ex1.{0}'.format(fext)) + a = ax.annotate( + "", + xy=(50, -25), + xytext=(30, -25), + arrowprops=dict(arrowstyle="->", fc="k"), + ) + ax.text( + 40, -22, "groundwater flow velocity=0.125 m/d", ha="center", size=7 + ) + ax.set_ylabel("Elevation, in meters") + ax.set_xlabel("Horizontal distance, in meters") + + outfig = os.path.join(workspace, "Figure06_swi2ex1.{0}".format(fext)) fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) return 0 diff --git a/examples/scripts/flopy_swi2_ex2.py b/examples/scripts/flopy_swi2_ex2.py index 8ae80b1d26..ec301d8fcc 100755 --- a/examples/scripts/flopy_swi2_ex2.py +++ b/examples/scripts/flopy_swi2_ex2.py @@ -11,57 +11,59 @@ import matplotlib.pyplot as plt # --modify default matplotlib settings -updates = {'font.family': ['Univers 57 Condensed', 'Arial'], - 'mathtext.default': 'regular', - 'pdf.compression': 0, - 'pdf.fonttype': 42, - 'legend.fontsize': 7, - 'axes.labelsize': 8, - 'xtick.labelsize': 7, - 'ytick.labelsize': 7} +updates = { + "font.family": ["Univers 57 Condensed", "Arial"], + "mathtext.default": "regular", + "pdf.compression": 0, + "pdf.fonttype": 42, + "legend.fontsize": 7, + "axes.labelsize": 8, + "xtick.labelsize": 7, + "ytick.labelsize": 7, +} plt.rcParams.update(updates) def run(): - workspace = 'swiex2' + workspace = "swiex2" if not os.path.isdir(workspace): os.mkdir(workspace) cleanFiles = False skipRuns = False - fext = 'png' + fext = "png" narg = len(sys.argv) iarg = 0 if narg > 1: while iarg < narg - 1: iarg += 1 basearg = sys.argv[iarg].lower() - if basearg == '--clean': + if basearg == "--clean": cleanFiles = True - elif basearg == '--skipruns': + elif basearg == "--skipruns": skipRuns = True - elif basearg == '--pdf': - fext = 'pdf' + elif basearg == "--pdf": + fext = "pdf" - dirs = [os.path.join(workspace, 'SWI2'), os.path.join(workspace, 'SEAWAT')] + dirs = [os.path.join(workspace, "SWI2"), os.path.join(workspace, "SEAWAT")] if cleanFiles: - print('cleaning all files') - print('excluding *.py files') + print("cleaning all files") + print("excluding *.py files") file_dict = collections.OrderedDict() file_dict[0] = os.listdir(dirs[0]) file_dict[1] = os.listdir(dirs[1]) file_dict[-1] = os.listdir(workspace) for key, files in list(file_dict.items()): - pth = '.' + pth = "." if key >= 0: pth = dirs[key] for f in files: fpth = os.path.join(pth, f) if os.path.isdir(fpth): continue - if '.py' != os.path.splitext(f)[1].lower(): - print(' removing...{}'.format(os.path.basename(f))) + if ".py" != os.path.splitext(f)[1].lower(): + print(" removing...{}".format(os.path.basename(f))) try: os.remove(fpth) except: @@ -76,15 +78,15 @@ def run(): if not os.path.exists(d): os.mkdir(d) - modelname = 'swiex2' - mf_name = 'mf2005' + modelname = "swiex2" + mf_name = "mf2005" # problem data nper = 1 perlen = 2000 nstp = 1000 nlay, nrow, ncol = 1, 1, 60 - delr = 5. + delr = 5.0 nsurf = 2 x = np.arange(0.5 * delr, ncol * delr, delr) xedge = np.linspace(0, float(ncol) * delr, len(x) + 1) @@ -101,75 +103,111 @@ def run(): z.append(z0) z.append(z1) ssz = 0.2 - isource = np.ones((nrow, ncol), 'int') + isource = np.ones((nrow, ncol), "int") isource[0, 0] = 2 # stratified model - modelname = 'swiex2_strat' - print('creating...', modelname) - ml = flopy.modflow.Modflow(modelname, version='mf2005', exe_name=mf_name, - model_ws=dirs[0]) - discret = flopy.modflow.ModflowDis(ml, nlay=1, ncol=ncol, nrow=nrow, - delr=delr, - delc=1, - top=0, botm=[-40.0], - nper=nper, perlen=perlen, nstp=nstp) + modelname = "swiex2_strat" + print("creating...", modelname) + ml = flopy.modflow.Modflow( + modelname, version="mf2005", exe_name=mf_name, model_ws=dirs[0] + ) + discret = flopy.modflow.ModflowDis( + ml, + nlay=1, + ncol=ncol, + nrow=nrow, + delr=delr, + delc=1, + top=0, + botm=[-40.0], + nper=nper, + perlen=perlen, + nstp=nstp, + ) bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=0.05) bcf = flopy.modflow.ModflowBcf(ml, laycon=0, tran=2 * 40) - swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, nsrf=nsurf, istrat=1, - toeslope=0.2, tipslope=0.2, - nu=[0, 0.0125, 0.025], - zeta=z, ssz=ssz, isource=isource, - nsolver=1) - oc = flopy.modflow.ModflowOc(ml, - stress_period_data={(0, 999): ['save head']}) + swi = flopy.modflow.ModflowSwi2( + ml, + iswizt=55, + nsrf=nsurf, + istrat=1, + toeslope=0.2, + tipslope=0.2, + nu=[0, 0.0125, 0.025], + zeta=z, + ssz=ssz, + isource=isource, + nsolver=1, + ) + oc = flopy.modflow.ModflowOc( + ml, stress_period_data={(0, 999): ["save head"]} + ) pcg = flopy.modflow.ModflowPcg(ml) ml.write_input() # run stratified model if not skipRuns: m = ml.run_model(silent=False) # read stratified results - zetafile = os.path.join(dirs[0], '{}.zta'.format(modelname)) + zetafile = os.path.join(dirs[0], "{}.zta".format(modelname)) zobj = flopy.utils.CellBudgetFile(zetafile) zkstpkper = zobj.get_kstpkper() - zeta = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF 1')[0] - zeta2 = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF 2')[0] + zeta = zobj.get_data(kstpkper=zkstpkper[-1], text="ZETASRF 1")[0] + zeta2 = zobj.get_data(kstpkper=zkstpkper[-1], text="ZETASRF 2")[0] # # vd model - modelname = 'swiex2_vd' - print('creating...', modelname) - ml = flopy.modflow.Modflow(modelname, version='mf2005', exe_name=mf_name, - model_ws=dirs[0]) - discret = flopy.modflow.ModflowDis(ml, nlay=1, ncol=ncol, nrow=nrow, - delr=delr, - delc=1, - top=0, botm=[-40.0], - nper=nper, perlen=perlen, nstp=nstp) + modelname = "swiex2_vd" + print("creating...", modelname) + ml = flopy.modflow.Modflow( + modelname, version="mf2005", exe_name=mf_name, model_ws=dirs[0] + ) + discret = flopy.modflow.ModflowDis( + ml, + nlay=1, + ncol=ncol, + nrow=nrow, + delr=delr, + delc=1, + top=0, + botm=[-40.0], + nper=nper, + perlen=perlen, + nstp=nstp, + ) bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=0.05) bcf = flopy.modflow.ModflowBcf(ml, laycon=0, tran=2 * 40) - swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, nsrf=nsurf, istrat=0, - toeslope=0.2, tipslope=0.2, - nu=[0, 0, 0.025, 0.025], - zeta=z, ssz=ssz, isource=isource, - nsolver=1) - oc = flopy.modflow.ModflowOc(ml, - stress_period_data={(0, 999): ['save head']}) + swi = flopy.modflow.ModflowSwi2( + ml, + iswizt=55, + nsrf=nsurf, + istrat=0, + toeslope=0.2, + tipslope=0.2, + nu=[0, 0, 0.025, 0.025], + zeta=z, + ssz=ssz, + isource=isource, + nsolver=1, + ) + oc = flopy.modflow.ModflowOc( + ml, stress_period_data={(0, 999): ["save head"]} + ) pcg = flopy.modflow.ModflowPcg(ml) ml.write_input() # run vd model if not skipRuns: m = ml.run_model(silent=False) # read vd model data - zetafile = os.path.join(dirs[0], '{}.zta'.format(modelname)) + zetafile = os.path.join(dirs[0], "{}.zta".format(modelname)) zobj = flopy.utils.CellBudgetFile(zetafile) zkstpkper = zobj.get_kstpkper() - zetavd = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF 1')[0] - zetavd2 = zobj.get_data(kstpkper=zkstpkper[-1], text='ZETASRF 2')[0] + zetavd = zobj.get_data(kstpkper=zkstpkper[-1], text="ZETASRF 1")[0] + zetavd2 = zobj.get_data(kstpkper=zkstpkper[-1], text="ZETASRF 2")[0] # # seawat model - swtexe_name = 'swtv4' - modelname = 'swiex2_swt' - print('creating...', modelname) + swtexe_name = "swtv4" + modelname = "swiex2_swt" + print("creating...", modelname) swt_xmax = 300.0 swt_zmax = 40.0 swt_delr = 1.0 @@ -184,7 +222,7 @@ def run(): swt_ibound[0, 0, 0] = -1 swt_x = np.arange(0.5 * swt_delr, swt_ncol * swt_delr, swt_delr) swt_xedge = np.linspace(0, float(ncol) * delr, len(swt_x) + 1) - swt_top = 0. + swt_top = 0.0 z0 = swt_top swt_botm = np.zeros((swt_nlay), np.float) swt_z = np.zeros((swt_nlay), np.float) @@ -217,60 +255,91 @@ def run(): # ssm data itype = flopy.mt3d.Mt3dSsm.itype_dict() - ssm_data = {0: [0, 0, 0, 35., itype['BAS6']]} + ssm_data = {0: [0, 0, 0, 35.0, itype["BAS6"]]} # print sconcp # mt3d print times - timprs = (np.arange(5) + 1) * 2000. + timprs = (np.arange(5) + 1) * 2000.0 nprs = len(timprs) # create the MODFLOW files m = flopy.seawat.Seawat(modelname, exe_name=swtexe_name, model_ws=dirs[1]) - discret = flopy.modflow.ModflowDis(m, nrow=swt_nrow, ncol=swt_ncol, - nlay=swt_nlay, - delr=swt_delr, delc=swt_delc, laycbd=0, - top=swt_top, - botm=swt_botm, - nper=nper, perlen=perlen, nstp=1, - steady=False) + discret = flopy.modflow.ModflowDis( + m, + nrow=swt_nrow, + ncol=swt_ncol, + nlay=swt_nlay, + delr=swt_delr, + delc=swt_delc, + laycbd=0, + top=swt_top, + botm=swt_botm, + nper=nper, + perlen=perlen, + nstp=1, + steady=False, + ) bas = flopy.modflow.ModflowBas(m, ibound=swt_ibound, strt=0.05) - lpf = flopy.modflow.ModflowLpf(m, hk=2.0, vka=2.0, ss=0.0, sy=0.0, - laytyp=0, - layavg=0) - oc = flopy.modflow.ModflowOc(m, save_every=1, save_types=['save head']) + lpf = flopy.modflow.ModflowLpf( + m, hk=2.0, vka=2.0, ss=0.0, sy=0.0, laytyp=0, layavg=0 + ) + oc = flopy.modflow.ModflowOc(m, save_every=1, save_types=["save head"]) pcg = flopy.modflow.ModflowPcg(m) # Create the MT3DMS model files - adv = flopy.mt3d.Mt3dAdv(m, mixelm=-1, # -1 is TVD - percel=0.05, - nadvfd=0, - # 0 or 1 is upstream; 2 is central in space - # particle based methods - nplane=4, - mxpart=1e7, - itrack=2, - dceps=1e-4, - npl=16, - nph=16, - npmin=8, - npmax=256) - btn = flopy.mt3d.Mt3dBtn(m, icbund=1, prsity=ssz, sconc=sconc, ifmtcn=-1, - chkmas=False, nprobs=10, nprmas=10, dt0=0.0, - ttsmult=1.2, - ttsmax=100.0, - ncomp=1, nprs=nprs, timprs=timprs, mxstrn=1e8) - dsp = flopy.mt3d.Mt3dDsp(m, al=0., trpt=1., trpv=1., dmcoef=0.) - gcg = flopy.mt3d.Mt3dGcg(m, mxiter=1, iter1=50, isolve=3, cclose=1e-6, - iprgcg=5) + adv = flopy.mt3d.Mt3dAdv( + m, + mixelm=-1, # -1 is TVD + percel=0.05, + nadvfd=0, + # 0 or 1 is upstream; 2 is central in space + # particle based methods + nplane=4, + mxpart=1e7, + itrack=2, + dceps=1e-4, + npl=16, + nph=16, + npmin=8, + npmax=256, + ) + btn = flopy.mt3d.Mt3dBtn( + m, + icbund=1, + prsity=ssz, + sconc=sconc, + ifmtcn=-1, + chkmas=False, + nprobs=10, + nprmas=10, + dt0=0.0, + ttsmult=1.2, + ttsmax=100.0, + ncomp=1, + nprs=nprs, + timprs=timprs, + mxstrn=1e8, + ) + dsp = flopy.mt3d.Mt3dDsp(m, al=0.0, trpt=1.0, trpv=1.0, dmcoef=0.0) + gcg = flopy.mt3d.Mt3dGcg( + m, mxiter=1, iter1=50, isolve=3, cclose=1e-6, iprgcg=5 + ) ssm = flopy.mt3d.Mt3dSsm(m, stress_period_data=ssm_data) # Create the SEAWAT model files - vdf = flopy.seawat.SeawatVdf(m, nswtcpl=1, iwtable=0, densemin=0, - densemax=0, - denseref=1000., denseslp=25., firstdt=1.0e-03) + vdf = flopy.seawat.SeawatVdf( + m, + nswtcpl=1, + iwtable=0, + densemin=0, + densemax=0, + denseref=1000.0, + denseslp=25.0, + firstdt=1.0e-03, + ) m.write_input() # run seawat model if not skipRuns: m = m.run_model(silent=False) # read seawat model data - ucnfile = os.path.join(dirs[1], 'MT3D001.UCN') + ucnfile = os.path.join(dirs[1], "MT3D001.UCN") uobj = flopy.utils.UcnFile(ucnfile) times = uobj.get_times() print(times) @@ -290,85 +359,125 @@ def run(): fbot = 0.125 ftop = 0.925 - print('creating cross-section figure...') - xsf = plt.figure(figsize=(fwid, fhgt), facecolor='w') - xsf.subplots_adjust(wspace=0.25, hspace=0.25, left=flft, right=frgt, - bottom=fbot, top=ftop) + print("creating cross-section figure...") + xsf = plt.figure(figsize=(fwid, fhgt), facecolor="w") + xsf.subplots_adjust( + wspace=0.25, hspace=0.25, left=flft, right=frgt, bottom=fbot, top=ftop + ) # plot initial conditions ax = xsf.add_subplot(3, 1, 1) - ax.text(-0.075, 1.05, 'A', transform=ax.transAxes, va='center', - ha='center', - size='8') + ax.text( + -0.075, + 1.05, + "A", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) # text(.975, .1, '(a)', transform = ax.transAxes, va = 'center', ha = 'center') - ax.plot([110, 150], [0, -40], 'k') - ax.plot([150, 190], [0, -40], 'k') + ax.plot([110, 150], [0, -40], "k") + ax.plot([150, 190], [0, -40], "k") ax.set_xlim(0, 300) ax.set_ylim(-40, 0) ax.set_yticks(np.arange(-40, 1, 10)) - ax.text(50, -20, 'salt', va='center', ha='center') - ax.text(150, -20, 'brackish', va='center', ha='center') - ax.text(250, -20, 'fresh', va='center', ha='center') - ax.set_ylabel('Elevation, in meters') + ax.text(50, -20, "salt", va="center", ha="center") + ax.text(150, -20, "brackish", va="center", ha="center") + ax.text(250, -20, "fresh", va="center", ha="center") + ax.set_ylabel("Elevation, in meters") # plot stratified swi2 and seawat results ax = xsf.add_subplot(3, 1, 2) - ax.text(-0.075, 1.05, 'B', transform=ax.transAxes, va='center', - ha='center', - size='8') + ax.text( + -0.075, + 1.05, + "B", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) # zp = zeta[0, 0, :] - p = (zp < 0.0) & (zp > -40.) - ax.plot(x[p], zp[p], 'b', linewidth=1.5, drawstyle='steps-mid') + p = (zp < 0.0) & (zp > -40.0) + ax.plot(x[p], zp[p], "b", linewidth=1.5, drawstyle="steps-mid") zp = zeta2[0, 0, :] - p = (zp < 0.0) & (zp > -40.) - ax.plot(x[p], zp[p], 'b', linewidth=1.5, drawstyle='steps-mid') + p = (zp < 0.0) & (zp > -40.0) + ax.plot(x[p], zp[p], "b", linewidth=1.5, drawstyle="steps-mid") # seawat data - cc = ax.contour(swt_X, swt_Z, conc, levels=[0.25, 0.75], colors='k', - linestyles='solid', linewidths=0.75, zorder=101) + cc = ax.contour( + swt_X, + swt_Z, + conc, + levels=[0.25, 0.75], + colors="k", + linestyles="solid", + linewidths=0.75, + zorder=101, + ) # fake figures - ax.plot([-100., -100], [-100., -100], 'b', linewidth=1.5, label='SWI2') - ax.plot([-100., -100], [-100., -100], 'k', linewidth=0.75, label='SEAWAT') + ax.plot([-100.0, -100], [-100.0, -100], "b", linewidth=1.5, label="SWI2") + ax.plot( + [-100.0, -100], [-100.0, -100], "k", linewidth=0.75, label="SEAWAT" + ) # legend - leg = ax.legend(loc='lower left', numpoints=1) + leg = ax.legend(loc="lower left", numpoints=1) leg._drawFrame = False # axes ax.set_xlim(0, 300) ax.set_ylim(-40, 0) ax.set_yticks(np.arange(-40, 1, 10)) - ax.set_ylabel('Elevation, in meters') + ax.set_ylabel("Elevation, in meters") # plot vd model ax = xsf.add_subplot(3, 1, 3) - ax.text(-0.075, 1.05, 'C', transform=ax.transAxes, va='center', - ha='center', size='8') + ax.text( + -0.075, + 1.05, + "C", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) dr = zeta[0, 0, :] - ax.plot(x, dr, 'b', linewidth=1.5, drawstyle='steps-mid') + ax.plot(x, dr, "b", linewidth=1.5, drawstyle="steps-mid") dr = zeta2[0, 0, :] - ax.plot(x, dr, 'b', linewidth=1.5, drawstyle='steps-mid') + ax.plot(x, dr, "b", linewidth=1.5, drawstyle="steps-mid") dr = zetavd[0, 0, :] - ax.plot(x, dr, 'r', linewidth=0.75, drawstyle='steps-mid') + ax.plot(x, dr, "r", linewidth=0.75, drawstyle="steps-mid") dr = zetavd2[0, 0, :] - ax.plot(x, dr, 'r', linewidth=0.75, drawstyle='steps-mid') + ax.plot(x, dr, "r", linewidth=0.75, drawstyle="steps-mid") # fake figures - ax.plot([-100., -100], [-100., -100], 'b', linewidth=1.5, - label='SWI2 stratified option') - ax.plot([-100., -100], [-100., -100], 'r', linewidth=0.75, - label='SWI2 continuous option') + ax.plot( + [-100.0, -100], + [-100.0, -100], + "b", + linewidth=1.5, + label="SWI2 stratified option", + ) + ax.plot( + [-100.0, -100], + [-100.0, -100], + "r", + linewidth=0.75, + label="SWI2 continuous option", + ) # legend - leg = ax.legend(loc='lower left', numpoints=1) + leg = ax.legend(loc="lower left", numpoints=1) leg._drawFrame = False # axes ax.set_xlim(0, 300) ax.set_ylim(-40, 0) ax.set_yticks(np.arange(-40, 1, 10)) - ax.set_xlabel('Horizontal distance, in meters') - ax.set_ylabel('Elevation, in meters') + ax.set_xlabel("Horizontal distance, in meters") + ax.set_ylabel("Elevation, in meters") - outfig = os.path.join(workspace, 'Figure07_swi2ex2.{0}'.format(fext)) + outfig = os.path.join(workspace, "Figure07_swi2ex2.{0}".format(fext)) xsf.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) return 0 -if __name__ == '__main__': +if __name__ == "__main__": success = run() sys.exit(success) diff --git a/examples/scripts/flopy_swi2_ex3.py b/examples/scripts/flopy_swi2_ex3.py index 358321b153..ea28b3bf21 100755 --- a/examples/scripts/flopy_swi2_ex3.py +++ b/examples/scripts/flopy_swi2_ex3.py @@ -10,14 +10,16 @@ import matplotlib.pyplot as plt # --modify default matplotlib settings -updates = {'font.family': ['Univers 57 Condensed', 'Arial'], - 'mathtext.default': 'regular', - 'pdf.compression': 0, - 'pdf.fonttype': 42, - 'legend.fontsize': 7, - 'axes.labelsize': 8, - 'xtick.labelsize': 7, - 'ytick.labelsize': 7} +updates = { + "font.family": ["Univers 57 Condensed", "Arial"], + "mathtext.default": "regular", + "pdf.compression": 0, + "pdf.fonttype": 42, + "legend.fontsize": 7, + "axes.labelsize": 8, + "xtick.labelsize": 7, + "ytick.labelsize": 7, +} plt.rcParams.update(updates) @@ -42,104 +44,120 @@ def MergeData(ndim, zdata, tb): def LegBar(ax, x0, y0, t0, dx, dy, dt, cc): for c in cc: ax.plot([x0, x0 + dx], [y0, y0], color=c, linewidth=4) - ctxt = '{0:=3d} years'.format(t0) - ax.text(x0 + 2. * dx, y0 + dy / 2., ctxt, size=5) + ctxt = "{0:=3d} years".format(t0) + ax.text(x0 + 2.0 * dx, y0 + dy / 2.0, ctxt, size=5) y0 += dy t0 += dt return def run(): - workspace = 'swiex3' + workspace = "swiex3" cleanFiles = False - fext = 'png' + fext = "png" narg = len(sys.argv) iarg = 0 if narg > 1: while iarg < narg - 1: iarg += 1 basearg = sys.argv[iarg].lower() - if basearg == '--clean': + if basearg == "--clean": cleanFiles = True - elif basearg == '--pdf': - fext = 'pdf' + elif basearg == "--pdf": + fext = "pdf" if cleanFiles: - print('cleaning all files') - print('excluding *.py files') + print("cleaning all files") + print("excluding *.py files") files = os.listdir(workspace) for f in files: fpth = os.path.join(workspace, f) if os.path.isdir(fpth): continue - if '.py' != os.path.splitext(f)[1].lower(): - print(' removing...{}'.format(os.path.basename(f))) + if ".py" != os.path.splitext(f)[1].lower(): + print(" removing...{}".format(os.path.basename(f))) try: os.remove(fpth) except: pass return 0 - modelname = 'swiex3' - exe_name = 'mf2005' + modelname = "swiex3" + exe_name = "mf2005" nlay = 3 nrow = 1 ncol = 200 delr = 20.0 - delc = 1. + delc = 1.0 # well data lrcQ1 = np.array([(0, 0, 199, 0.01), (2, 0, 199, 0.02)]) lrcQ2 = np.array([(0, 0, 199, 0.01 * 0.5), (2, 0, 199, 0.02 * 0.5)]) # ghb data lrchc = np.zeros((30, 5)) - lrchc[:, [0, 1, 3, 4]] = [0, 0, 0., 0.8 / 2.0] + lrchc[:, [0, 1, 3, 4]] = [0, 0, 0.0, 0.8 / 2.0] lrchc[:, 2] = np.arange(0, 30) # swi2 data zini = np.hstack( - (-9 * np.ones(24), np.arange(-9, -50, -0.5), -50 * np.ones(94)))[ - np.newaxis, :] + (-9 * np.ones(24), np.arange(-9, -50, -0.5), -50 * np.ones(94)) + )[np.newaxis, :] iso = np.zeros((1, 200), dtype=np.int) iso[:, :30] = -2 # model objects - ml = flopy.modflow.Modflow(modelname, version='mf2005', exe_name=exe_name, - model_ws=workspace) - discret = flopy.modflow.ModflowDis(ml, nrow=nrow, ncol=ncol, nlay=3, - delr=delr, - delc=delc, - laycbd=[0, 0, 0], top=-9.0, - botm=[-29, -30, -50], - nper=2, perlen=[365 * 1000, 1000 * 365], - nstp=[500, 500]) + ml = flopy.modflow.Modflow( + modelname, version="mf2005", exe_name=exe_name, model_ws=workspace + ) + discret = flopy.modflow.ModflowDis( + ml, + nrow=nrow, + ncol=ncol, + nlay=3, + delr=delr, + delc=delc, + laycbd=[0, 0, 0], + top=-9.0, + botm=[-29, -30, -50], + nper=2, + perlen=[365 * 1000, 1000 * 365], + nstp=[500, 500], + ) bas = flopy.modflow.ModflowBas(ml, ibound=1, strt=1.0) - bcf = flopy.modflow.ModflowBcf(ml, laycon=[0, 0, 0], tran=[40.0, 1, 80.0], - vcont=[0.005, 0.005]) + bcf = flopy.modflow.ModflowBcf( + ml, laycon=[0, 0, 0], tran=[40.0, 1, 80.0], vcont=[0.005, 0.005] + ) wel = flopy.modflow.ModflowWel(ml, stress_period_data={0: lrcQ1, 1: lrcQ2}) ghb = flopy.modflow.ModflowGhb(ml, stress_period_data={0: lrchc}) - swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, nsrf=1, istrat=1, - toeslope=0.01, tipslope=0.04, - nu=[0, 0.025], - zeta=[zini, zini, zini], ssz=0.2, - isource=iso, - nsolver=1) - oc = flopy.modflow.ModflowOc(ml, save_every=100, save_types=['save head']) + swi = flopy.modflow.ModflowSwi2( + ml, + iswizt=55, + nsrf=1, + istrat=1, + toeslope=0.01, + tipslope=0.04, + nu=[0, 0.025], + zeta=[zini, zini, zini], + ssz=0.2, + isource=iso, + nsolver=1, + ) + oc = flopy.modflow.ModflowOc(ml, save_every=100, save_types=["save head"]) pcg = flopy.modflow.ModflowPcg(ml) # write the model files ml.write_input() # run the model m = ml.run_model(silent=True) - headfile = os.path.join(workspace, '{}.hds'.format(modelname)) + headfile = os.path.join(workspace, "{}.hds".format(modelname)) hdobj = flopy.utils.HeadFile(headfile) - head = hdobj.get_data(totim=3.65000E+05) + head = hdobj.get_data(totim=3.65000e05) - zetafile = os.path.join(workspace, '{}.zta'.format(modelname)) + zetafile = os.path.join(workspace, "{}.zta".format(modelname)) zobj = flopy.utils.CellBudgetFile(zetafile) zkstpkper = zobj.get_kstpkper() zeta = [] for kk in zkstpkper: - zeta.append(zobj.get_data(kstpkper=kk, text='ZETASRF 1')[0]) + zeta.append(zobj.get_data(kstpkper=kk, text="ZETASRF 1")[0]) zeta = np.array(zeta) fwid, fhgt = 7.00, 4.50 @@ -154,106 +172,150 @@ def run(): lw = 0.5 x = np.arange(-30 * delr + 0.5 * delr, (ncol - 30) * delr, delr) - xedge = np.linspace(-30. * delr, (ncol - 30.) * delr, len(x) + 1) - zedge = [[-9., -29.], [-29., -30.], [-30., -50.]] + xedge = np.linspace(-30.0 * delr, (ncol - 30.0) * delr, len(x) + 1) + zedge = [[-9.0, -29.0], [-29.0, -30.0], [-30.0, -50.0]] - fig = plt.figure(figsize=(fwid, fhgt), facecolor='w') - fig.subplots_adjust(wspace=0.25, hspace=0.25, left=flft, right=frgt, - bottom=fbot, top=ftop) + fig = plt.figure(figsize=(fwid, fhgt), facecolor="w") + fig.subplots_adjust( + wspace=0.25, hspace=0.25, left=flft, right=frgt, bottom=fbot, top=ftop + ) ax = fig.add_subplot(311) - ax.text(-0.075, 1.05, 'A', transform=ax.transAxes, va='center', - ha='center', - size='8') + ax.text( + -0.075, + 1.05, + "A", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) # confining unit - ax.fill([-600, 3400, 3400, -600], [-29, -29, -30, -30], fc=[.8, .8, .8], - ec=[.8, .8, .8]) + ax.fill( + [-600, 3400, 3400, -600], + [-29, -29, -30, -30], + fc=[0.8, 0.8, 0.8], + ec=[0.8, 0.8, 0.8], + ) # z = np.copy(zini[0, :]) zr = z.copy() - p = (zr < -9.) & (zr > -50.0) - ax.plot(x[p], zr[p], color=cc[0], linewidth=lw, drawstyle='steps-mid') + p = (zr < -9.0) & (zr > -50.0) + ax.plot(x[p], zr[p], color=cc[0], linewidth=lw, drawstyle="steps-mid") # for i in range(5): - zt = MergeData(ncol, - [zeta[i, 0, 0, :], zeta[i, 1, 0, :], zeta[i, 2, 0, :]], - zedge) + zt = MergeData( + ncol, [zeta[i, 0, 0, :], zeta[i, 1, 0, :], zeta[i, 2, 0, :]], zedge + ) dr = zt.copy() - ax.plot(x, dr, color=cc[i + 1], linewidth=lw, drawstyle='steps-mid') + ax.plot(x, dr, color=cc[i + 1], linewidth=lw, drawstyle="steps-mid") # Manufacture a legend bar - LegBar(ax, -200., -33.75, 0, 25, -2.5, 200, cc[0:6]) + LegBar(ax, -200.0, -33.75, 0, 25, -2.5, 200, cc[0:6]) # axes ax.set_ylim(-50, -9) - ax.set_ylabel('Elevation, in meters') - ax.set_xlim(-250., 2500.) + ax.set_ylabel("Elevation, in meters") + ax.set_xlim(-250.0, 2500.0) ax = fig.add_subplot(312) - ax.text(-0.075, 1.05, 'B', transform=ax.transAxes, va='center', - ha='center', - size='8') + ax.text( + -0.075, + 1.05, + "B", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) # confining unit - ax.fill([-600, 3400, 3400, -600], [-29, -29, -30, -30], fc=[.8, .8, .8], - ec=[.8, .8, .8]) + ax.fill( + [-600, 3400, 3400, -600], + [-29, -29, -30, -30], + fc=[0.8, 0.8, 0.8], + ec=[0.8, 0.8, 0.8], + ) # for i in range(4, 10): - zt = MergeData(ncol, - [zeta[i, 0, 0, :], zeta[i, 1, 0, :], zeta[i, 2, 0, :]], - zedge) + zt = MergeData( + ncol, [zeta[i, 0, 0, :], zeta[i, 1, 0, :], zeta[i, 2, 0, :]], zedge + ) dr = zt.copy() - ax.plot(x, dr, color=cc[i + 1], linewidth=lw, drawstyle='steps-mid') + ax.plot(x, dr, color=cc[i + 1], linewidth=lw, drawstyle="steps-mid") # Manufacture a legend bar - LegBar(ax, -200., -33.75, 1000, 25, -2.5, 200, cc[5:11]) + LegBar(ax, -200.0, -33.75, 1000, 25, -2.5, 200, cc[5:11]) # axes ax.set_ylim(-50, -9) - ax.set_ylabel('Elevation, in meters') - ax.set_xlim(-250., 2500.) + ax.set_ylabel("Elevation, in meters") + ax.set_xlim(-250.0, 2500.0) ax = fig.add_subplot(313) - ax.text(-0.075, 1.05, 'C', transform=ax.transAxes, va='center', - ha='center', - size='8') + ax.text( + -0.075, + 1.05, + "C", + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) # confining unit - ax.fill([-600, 3400, 3400, -600], [-29, -29, -30, -30], fc=[.8, .8, .8], - ec=[.8, .8, .8]) + ax.fill( + [-600, 3400, 3400, -600], + [-29, -29, -30, -30], + fc=[0.8, 0.8, 0.8], + ec=[0.8, 0.8, 0.8], + ) # - zt = MergeData(ncol, - [zeta[4, 0, 0, :], zeta[4, 1, 0, :], zeta[4, 2, 0, :]], - zedge) - ax.plot(x, zt, marker='o', markersize=3, linewidth=0.0, - markeredgecolor='blue', - markerfacecolor='None') + zt = MergeData( + ncol, [zeta[4, 0, 0, :], zeta[4, 1, 0, :], zeta[4, 2, 0, :]], zedge + ) + ax.plot( + x, + zt, + marker="o", + markersize=3, + linewidth=0.0, + markeredgecolor="blue", + markerfacecolor="None", + ) # ghyben herzberg - zeta1 = -9 - 40. * (head[0, 0, :]) + zeta1 = -9 - 40.0 * (head[0, 0, :]) gbh = np.empty(len(zeta1), np.float) gbho = np.empty(len(zeta1), np.float) for idx, z1 in enumerate(zeta1): if z1 >= -9.0 or z1 <= -50.0: gbh[idx] = np.nan - gbho[idx] = 0. + gbho[idx] = 0.0 else: gbh[idx] = z1 gbho[idx] = z1 - ax.plot(x, gbh, 'r') - np.savetxt(os.path.join(workspace, 'Ghyben-Herzberg.out'), gbho) + ax.plot(x, gbh, "r") + np.savetxt(os.path.join(workspace, "Ghyben-Herzberg.out"), gbho) # fake figures - ax.plot([-100., -100], [-100., -100], 'r', label='Ghyben-Herzberg') - ax.plot([-100., -100], [-100., -100], 'bo', markersize=3, - markeredgecolor='blue', markerfacecolor='None', label='SWI2') + ax.plot([-100.0, -100], [-100.0, -100], "r", label="Ghyben-Herzberg") + ax.plot( + [-100.0, -100], + [-100.0, -100], + "bo", + markersize=3, + markeredgecolor="blue", + markerfacecolor="None", + label="SWI2", + ) # legend - leg = ax.legend(loc='lower left', numpoints=1) + leg = ax.legend(loc="lower left", numpoints=1) leg._drawFrame = False # axes ax.set_ylim(-50, -9) - ax.set_xlabel('Horizontal distance, in meters') - ax.set_ylabel('Elevation, in meters') - ax.set_xlim(-250., 2500.) + ax.set_xlabel("Horizontal distance, in meters") + ax.set_ylabel("Elevation, in meters") + ax.set_xlim(-250.0, 2500.0) - outfig = os.path.join(workspace, 'Figure08_swi2ex3.{0}'.format(fext)) + outfig = os.path.join(workspace, "Figure08_swi2ex3.{0}".format(fext)) fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) return 0 -if __name__ == '__main__': +if __name__ == "__main__": success = run() diff --git a/examples/scripts/flopy_swi2_ex4.py b/examples/scripts/flopy_swi2_ex4.py index bacfa915bb..f1a8842dc1 100644 --- a/examples/scripts/flopy_swi2_ex4.py +++ b/examples/scripts/flopy_swi2_ex4.py @@ -11,56 +11,58 @@ import matplotlib.pyplot as plt # --modify default matplotlib settings -updates = {'font.family': ['Univers 57 Condensed', 'Arial'], - 'mathtext.default': 'regular', - 'pdf.compression': 0, - 'pdf.fonttype': 42, - 'legend.fontsize': 7, - 'axes.labelsize': 8, - 'xtick.labelsize': 7, - 'ytick.labelsize': 7} +updates = { + "font.family": ["Univers 57 Condensed", "Arial"], + "mathtext.default": "regular", + "pdf.compression": 0, + "pdf.fonttype": 42, + "legend.fontsize": 7, + "axes.labelsize": 8, + "xtick.labelsize": 7, + "ytick.labelsize": 7, +} plt.rcParams.update(updates) def LegBar(ax, x0, y0, t0, dx, dy, dt, cc): for c in cc: ax.plot([x0, x0 + dx], [y0, y0], color=c, linewidth=4) - ctxt = '{0:=3d} years'.format(t0) - ax.text(x0 + 2. * dx, y0 + dy / 2., ctxt, size=5) + ctxt = "{0:=3d} years".format(t0) + ax.text(x0 + 2.0 * dx, y0 + dy / 2.0, ctxt, size=5) y0 += dy t0 += dt return def run(): - workspace = 'swiex4' + workspace = "swiex4" cleanFiles = False skipRuns = False - fext = 'png' + fext = "png" narg = len(sys.argv) iarg = 0 if narg > 1: while iarg < narg - 1: iarg += 1 basearg = sys.argv[iarg].lower() - if basearg == '--clean': + if basearg == "--clean": cleanFiles = True - elif basearg == '--skipruns': + elif basearg == "--skipruns": skipRuns = True - elif basearg == '--pdf': - fext = 'pdf' + elif basearg == "--pdf": + fext = "pdf" if cleanFiles: - print('cleaning all files') - print('excluding *.py files') + print("cleaning all files") + print("excluding *.py files") files = os.listdir(workspace) for f in files: fpth = os.path.join(workspace, f) if os.path.isdir(fpth): continue - if '.py' != os.path.splitext(f)[1].lower(): - print(' removing...{}'.format(os.path.basename(f))) + if ".py" != os.path.splitext(f)[1].lower(): + print(" removing...{}".format(os.path.basename(f))) try: os.remove(fpth) except: @@ -68,45 +70,45 @@ def run(): return 0 # Set path and name of MODFLOW exe - exe_name = 'mf2005' - if platform.system() == 'Windows': - exe_name = 'mf2005' + exe_name = "mf2005" + if platform.system() == "Windows": + exe_name = "mf2005" ncol = 61 nrow = 61 nlay = 2 nper = 3 - perlen = [365.25 * 200., 365.25 * 12., 365.25 * 18.] + perlen = [365.25 * 200.0, 365.25 * 12.0, 365.25 * 18.0] nstp = [1000, 120, 180] save_head = [200, 60, 60] steady = True # dis data delr, delc = 50.0, 50.0 - botm = np.array([-10., -30., -50.]) + botm = np.array([-10.0, -30.0, -50.0]) # oc data ocspd = {} - ocspd[(0, 199)] = ['save head'] + ocspd[(0, 199)] = ["save head"] ocspd[(0, 200)] = [] - ocspd[(0, 399)] = ['save head'] + ocspd[(0, 399)] = ["save head"] ocspd[(0, 400)] = [] - ocspd[(0, 599)] = ['save head'] + ocspd[(0, 599)] = ["save head"] ocspd[(0, 600)] = [] - ocspd[(0, 799)] = ['save head'] + ocspd[(0, 799)] = ["save head"] ocspd[(0, 800)] = [] - ocspd[(0, 999)] = ['save head'] + ocspd[(0, 999)] = ["save head"] ocspd[(1, 0)] = [] - ocspd[(1, 59)] = ['save head'] + ocspd[(1, 59)] = ["save head"] ocspd[(1, 60)] = [] - ocspd[(1, 119)] = ['save head'] + ocspd[(1, 119)] = ["save head"] ocspd[(2, 0)] = [] - ocspd[(2, 59)] = ['save head'] + ocspd[(2, 59)] = ["save head"] ocspd[(2, 60)] = [] - ocspd[(2, 119)] = ['save head'] + ocspd[(2, 119)] = ["save head"] - modelname = 'swiex4_2d_2layer' + modelname = "swiex4_2d_2layer" # bas data # ibound - active except for the corners @@ -120,7 +122,7 @@ def run(): # lpf data laytyp = 0 - hk = 10. + hk = 10.0 vka = 0.2 # boundary condition data @@ -136,7 +138,7 @@ def run(): lrchc[:, 0] = 0 lrchc[:, 1] = rowcell[index == 1] lrchc[:, 2] = colcell[index == 1] - lrchc[:, 3] = 0. + lrchc[:, 3] = 0.0 lrchc[:, 4] = 50.0 * 50.0 / 40.0 # create ghb dictionary ghb_data = {0: lrchc} @@ -155,8 +157,8 @@ def run(): lrcqw = lrcq.copy() lrcqw[0, 3] = -250 lrcqsw = lrcq.copy() - lrcqsw[0, 3] = -250. - lrcqsw[1, 3] = -25. + lrcqsw[0, 3] = -250.0 + lrcqsw[1, 3] = -25.0 # create well dictionary base_well_data = {0: lrcq, 1: lrcqw} swwells_well_data = {0: lrcq, 1: lrcqw, 2: lrcqsw} @@ -178,78 +180,117 @@ def run(): iso[1, 30, 35] = 2 ssz = 0.2 # swi2 observations - obsnam = ['layer1_', 'layer2_'] + obsnam = ["layer1_", "layer2_"] obslrc = [[0, 30, 35], [1, 30, 35]] nobs = len(obsnam) iswiobs = 1051 - modelname = 'swiex4_s1' + modelname = "swiex4_s1" if not skipRuns: - ml = flopy.modflow.Modflow(modelname, version='mf2005', - exe_name=exe_name, - model_ws=workspace) - - discret = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol, - laycbd=0, - delr=delr, delc=delc, top=botm[0], - botm=botm[1:], - nper=nper, perlen=perlen, nstp=nstp, - steady=steady) + ml = flopy.modflow.Modflow( + modelname, version="mf2005", exe_name=exe_name, model_ws=workspace + ) + + discret = flopy.modflow.ModflowDis( + ml, + nlay=nlay, + nrow=nrow, + ncol=ncol, + laycbd=0, + delr=delr, + delc=delc, + top=botm[0], + botm=botm[1:], + nper=nper, + perlen=perlen, + nstp=nstp, + steady=steady, + ) bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=ihead) lpf = flopy.modflow.ModflowLpf(ml, laytyp=laytyp, hk=hk, vka=vka) wel = flopy.modflow.ModflowWel(ml, stress_period_data=base_well_data) ghb = flopy.modflow.ModflowGhb(ml, stress_period_data=ghb_data) rch = flopy.modflow.ModflowRch(ml, rech=rch_data) - swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, nsrf=1, istrat=1, - toeslope=toeslope, - tipslope=tipslope, nu=nu, - zeta=z, ssz=ssz, isource=iso, - nsolver=1, - nadptmx=nadptmx, nadptmn=nadptmn, - nobs=nobs, iswiobs=iswiobs, - obsnam=obsnam, - obslrc=obslrc) + swi = flopy.modflow.ModflowSwi2( + ml, + iswizt=55, + nsrf=1, + istrat=1, + toeslope=toeslope, + tipslope=tipslope, + nu=nu, + zeta=z, + ssz=ssz, + isource=iso, + nsolver=1, + nadptmx=nadptmx, + nadptmn=nadptmn, + nobs=nobs, + iswiobs=iswiobs, + obsnam=obsnam, + obslrc=obslrc, + ) oc = flopy.modflow.ModflowOc(ml, stress_period_data=ocspd) - pcg = flopy.modflow.ModflowPcg(ml, hclose=1.0e-6, rclose=3.0e-3, - mxiter=100, iter1=50) + pcg = flopy.modflow.ModflowPcg( + ml, hclose=1.0e-6, rclose=3.0e-3, mxiter=100, iter1=50 + ) # create model files ml.write_input() # run the model m = ml.run_model(silent=False) # model with saltwater wells - modelname2 = 'swiex4_s2' + modelname2 = "swiex4_s2" if not skipRuns: - ml2 = flopy.modflow.Modflow(modelname2, version='mf2005', - exe_name=exe_name, - model_ws=workspace) - - discret = flopy.modflow.ModflowDis(ml2, nlay=nlay, nrow=nrow, - ncol=ncol, - laycbd=0, - delr=delr, delc=delc, top=botm[0], - botm=botm[1:], - nper=nper, perlen=perlen, nstp=nstp, - steady=steady) + ml2 = flopy.modflow.Modflow( + modelname2, version="mf2005", exe_name=exe_name, model_ws=workspace + ) + + discret = flopy.modflow.ModflowDis( + ml2, + nlay=nlay, + nrow=nrow, + ncol=ncol, + laycbd=0, + delr=delr, + delc=delc, + top=botm[0], + botm=botm[1:], + nper=nper, + perlen=perlen, + nstp=nstp, + steady=steady, + ) bas = flopy.modflow.ModflowBas(ml2, ibound=ibound, strt=ihead) lpf = flopy.modflow.ModflowLpf(ml2, laytyp=laytyp, hk=hk, vka=vka) - wel = flopy.modflow.ModflowWel(ml2, - stress_period_data=swwells_well_data) + wel = flopy.modflow.ModflowWel( + ml2, stress_period_data=swwells_well_data + ) ghb = flopy.modflow.ModflowGhb(ml2, stress_period_data=ghb_data) rch = flopy.modflow.ModflowRch(ml2, rech=rch_data) - swi = flopy.modflow.ModflowSwi2(ml2, iswizt=55, nsrf=1, istrat=1, - toeslope=toeslope, - tipslope=tipslope, nu=nu, - zeta=z, ssz=ssz, isource=iso, - nsolver=1, - nadptmx=nadptmx, nadptmn=nadptmn, - nobs=nobs, iswiobs=iswiobs, - obsnam=obsnam, - obslrc=obslrc) + swi = flopy.modflow.ModflowSwi2( + ml2, + iswizt=55, + nsrf=1, + istrat=1, + toeslope=toeslope, + tipslope=tipslope, + nu=nu, + zeta=z, + ssz=ssz, + isource=iso, + nsolver=1, + nadptmx=nadptmx, + nadptmn=nadptmn, + nobs=nobs, + iswiobs=iswiobs, + obsnam=obsnam, + obslrc=obslrc, + ) oc = flopy.modflow.ModflowOc(ml2, stress_period_data=ocspd) - pcg = flopy.modflow.ModflowPcg(ml2, hclose=1.0e-6, rclose=3.0e-3, - mxiter=100, - iter1=50) + pcg = flopy.modflow.ModflowPcg( + ml2, hclose=1.0e-6, rclose=3.0e-3, mxiter=100, iter1=50 + ) # create model files ml2.write_input() # run the model @@ -258,33 +299,37 @@ def run(): # Load the simulation 1 `ZETA` data and `ZETA` observations. # read base model zeta zfile = flopy.utils.CellBudgetFile( - os.path.join(workspace, modelname + '.zta')) + os.path.join(workspace, modelname + ".zta") + ) kstpkper = zfile.get_kstpkper() zeta = [] for kk in kstpkper: - zeta.append(zfile.get_data(kstpkper=kk, text='ZETASRF 1')[0]) + zeta.append(zfile.get_data(kstpkper=kk, text="ZETASRF 1")[0]) zeta = np.array(zeta) # read swi obs - zobs = np.genfromtxt(os.path.join(workspace, modelname + '.zobs.out'), - names=True) + zobs = np.genfromtxt( + os.path.join(workspace, modelname + ".zobs.out"), names=True + ) # Load the simulation 2 `ZETA` data and `ZETA` observations. # read saltwater well model zeta zfile2 = flopy.utils.CellBudgetFile( - os.path.join(workspace, modelname2 + '.zta')) + os.path.join(workspace, modelname2 + ".zta") + ) kstpkper = zfile2.get_kstpkper() zeta2 = [] for kk in kstpkper: - zeta2.append(zfile2.get_data(kstpkper=kk, text='ZETASRF 1')[0]) + zeta2.append(zfile2.get_data(kstpkper=kk, text="ZETASRF 1")[0]) zeta2 = np.array(zeta2) # read swi obs - zobs2 = np.genfromtxt(os.path.join(workspace, modelname2 + '.zobs.out'), - names=True) + zobs2 = np.genfromtxt( + os.path.join(workspace, modelname2 + ".zobs.out"), names=True + ) # Create arrays for the x-coordinates and the output years x = np.linspace(-1500, 1500, 61) - xcell = np.linspace(-1500, 1500, 61) + delr / 2. + xcell = np.linspace(-1500, 1500, 61) + delr / 2.0 xedge = np.linspace(-1525, 1525, 62) years = [40, 80, 120, 160, 200, 6, 12, 18, 24, 30] @@ -304,10 +349,11 @@ def run(): # Recreate **Figure 9** from the SWI2 documentation (http://pubs.usgs.gov/tm/6a46/). - plt.rcParams.update({'legend.fontsize': 6, 'legend.frameon': False}) - fig = plt.figure(figsize=(fwid, fhgt), facecolor='w') - fig.subplots_adjust(wspace=0.25, hspace=0.25, left=flft, right=frgt, - bottom=fbot, top=ftop) + plt.rcParams.update({"legend.fontsize": 6, "legend.frameon": False}) + fig = plt.figure(figsize=(fwid, fhgt), facecolor="w") + fig.subplots_adjust( + wspace=0.25, hspace=0.25, left=flft, right=frgt, bottom=fbot, top=ftop + ) # first plot ax = fig.add_subplot(2, 2, 1) # axes limits @@ -315,27 +361,56 @@ def run(): ax.set_ylim(-50, -10) for idx in range(5): # layer 1 - ax.plot(xcell, zeta[idx, 0, 30, :], drawstyle='steps-mid', - linewidth=0.5, color=cc[idx], - label='{:2d} years'.format(years[idx])) + ax.plot( + xcell, + zeta[idx, 0, 30, :], + drawstyle="steps-mid", + linewidth=0.5, + color=cc[idx], + label="{:2d} years".format(years[idx]), + ) # layer 2 - ax.plot(xcell, zeta[idx, 1, 30, :], drawstyle='steps-mid', - linewidth=0.5, color=cc[idx], label='_None') - ax.plot([-1500, 1500], [-30, -30], color='k', linewidth=1.0) + ax.plot( + xcell, + zeta[idx, 1, 30, :], + drawstyle="steps-mid", + linewidth=0.5, + color=cc[idx], + label="_None", + ) + ax.plot([-1500, 1500], [-30, -30], color="k", linewidth=1.0) # legend - plt.legend(loc='lower left') + plt.legend(loc="lower left") # axes labels and text - ax.set_xlabel('Horizontal distance, in meters') - ax.set_ylabel('Elevation, in meters') - ax.text(0.025, .55, 'Layer 1', transform=ax.transAxes, va='center', - ha='left', - size='7') - ax.text(0.025, .45, 'Layer 2', transform=ax.transAxes, va='center', - ha='left', - size='7') - ax.text(0.975, .1, 'Recharge conditions', transform=ax.transAxes, - va='center', - ha='right', size='8') + ax.set_xlabel("Horizontal distance, in meters") + ax.set_ylabel("Elevation, in meters") + ax.text( + 0.025, + 0.55, + "Layer 1", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + ax.text( + 0.025, + 0.45, + "Layer 2", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + ax.text( + 0.975, + 0.1, + "Recharge conditions", + transform=ax.transAxes, + va="center", + ha="right", + size="8", + ) # second plot ax = fig.add_subplot(2, 2, 2) @@ -344,26 +419,56 @@ def run(): ax.set_ylim(-50, -10) for idx in range(5, len(years) - 1): # layer 1 - ax.plot(xcell, zeta[idx, 0, 30, :], drawstyle='steps-mid', - linewidth=0.5, color=cc[idx - 5], - label='{:2d} years'.format(years[idx])) + ax.plot( + xcell, + zeta[idx, 0, 30, :], + drawstyle="steps-mid", + linewidth=0.5, + color=cc[idx - 5], + label="{:2d} years".format(years[idx]), + ) # layer 2 - ax.plot(xcell, zeta[idx, 1, 30, :], drawstyle='steps-mid', - linewidth=0.5, color=cc[idx - 5], label='_None') - ax.plot([-1500, 1500], [-30, -30], color='k', linewidth=1.0) + ax.plot( + xcell, + zeta[idx, 1, 30, :], + drawstyle="steps-mid", + linewidth=0.5, + color=cc[idx - 5], + label="_None", + ) + ax.plot([-1500, 1500], [-30, -30], color="k", linewidth=1.0) # legend - plt.legend(loc='lower left') + plt.legend(loc="lower left") # axes labels and text - ax.set_xlabel('Horizontal distance, in meters') - ax.set_ylabel('Elevation, in meters') - ax.text(0.025, .55, 'Layer 1', transform=ax.transAxes, va='center', - ha='left', - size='7') - ax.text(0.025, .45, 'Layer 2', transform=ax.transAxes, va='center', - ha='left', - size='7') - ax.text(0.975, .1, 'Freshwater well withdrawal', transform=ax.transAxes, - va='center', ha='right', size='8') + ax.set_xlabel("Horizontal distance, in meters") + ax.set_ylabel("Elevation, in meters") + ax.text( + 0.025, + 0.55, + "Layer 1", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + ax.text( + 0.025, + 0.45, + "Layer 2", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + ax.text( + 0.975, + 0.1, + "Freshwater well withdrawal", + transform=ax.transAxes, + va="center", + ha="right", + size="8", + ) # third plot ax = fig.add_subplot(2, 2, 3) @@ -372,64 +477,117 @@ def run(): ax.set_ylim(-50, -10) for idx in range(5, len(years) - 1): # layer 1 - ax.plot(xcell, zeta2[idx, 0, 30, :], drawstyle='steps-mid', - linewidth=0.5, color=cc[idx - 5], - label='{:2d} years'.format(years[idx])) + ax.plot( + xcell, + zeta2[idx, 0, 30, :], + drawstyle="steps-mid", + linewidth=0.5, + color=cc[idx - 5], + label="{:2d} years".format(years[idx]), + ) # layer 2 - ax.plot(xcell, zeta2[idx, 1, 30, :], drawstyle='steps-mid', - linewidth=0.5, color=cc[idx - 5], label='_None') - ax.plot([-1500, 1500], [-30, -30], color='k', linewidth=1.0) + ax.plot( + xcell, + zeta2[idx, 1, 30, :], + drawstyle="steps-mid", + linewidth=0.5, + color=cc[idx - 5], + label="_None", + ) + ax.plot([-1500, 1500], [-30, -30], color="k", linewidth=1.0) # legend - plt.legend(loc='lower left') + plt.legend(loc="lower left") # axes labels and text - ax.set_xlabel('Horizontal distance, in meters') - ax.set_ylabel('Elevation, in meters') - ax.text(0.025, .55, 'Layer 1', transform=ax.transAxes, va='center', - ha='left', - size='7') - ax.text(0.025, .45, 'Layer 2', transform=ax.transAxes, va='center', - ha='left', - size='7') - ax.text(0.975, .1, 'Freshwater and saltwater\nwell withdrawals', - transform=ax.transAxes, - va='center', ha='right', size='8') + ax.set_xlabel("Horizontal distance, in meters") + ax.set_ylabel("Elevation, in meters") + ax.text( + 0.025, + 0.55, + "Layer 1", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + ax.text( + 0.025, + 0.45, + "Layer 2", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + ax.text( + 0.975, + 0.1, + "Freshwater and saltwater\nwell withdrawals", + transform=ax.transAxes, + va="center", + ha="right", + size="8", + ) # fourth plot ax = fig.add_subplot(2, 2, 4) # axes limits ax.set_xlim(0, 30) ax.set_ylim(-50, -10) - t = zobs['TOTIM'][999:] / 365 - 200. - tz2 = zobs['layer1_001'][999:] - tz3 = zobs2['layer1_001'][999:] + t = zobs["TOTIM"][999:] / 365 - 200.0 + tz2 = zobs["layer1_001"][999:] + tz3 = zobs2["layer1_001"][999:] for i in range(len(t)): - if zobs['layer2_001'][i + 999] < -30. - 0.1: - tz2[i] = zobs['layer2_001'][i + 999] - if zobs2['layer2_001'][i + 999] < 20. - 0.1: - tz3[i] = zobs2['layer2_001'][i + 999] - ax.plot(t, tz2, linestyle='solid', color='r', linewidth=0.75, - label='Freshwater well') - ax.plot(t, tz3, linestyle='dotted', color='r', linewidth=0.75, - label='Freshwater and saltwater well') - ax.plot([0, 30], [-30, -30], 'k', linewidth=1.0, label='_None') + if zobs["layer2_001"][i + 999] < -30.0 - 0.1: + tz2[i] = zobs["layer2_001"][i + 999] + if zobs2["layer2_001"][i + 999] < 20.0 - 0.1: + tz3[i] = zobs2["layer2_001"][i + 999] + ax.plot( + t, + tz2, + linestyle="solid", + color="r", + linewidth=0.75, + label="Freshwater well", + ) + ax.plot( + t, + tz3, + linestyle="dotted", + color="r", + linewidth=0.75, + label="Freshwater and saltwater well", + ) + ax.plot([0, 30], [-30, -30], "k", linewidth=1.0, label="_None") # legend - leg = plt.legend(loc='lower right', numpoints=1) + leg = plt.legend(loc="lower right", numpoints=1) # axes labels and text - ax.set_xlabel('Time, in years') - ax.set_ylabel('Elevation, in meters') - ax.text(0.025, .55, 'Layer 1', transform=ax.transAxes, va='center', - ha='left', - size='7') - ax.text(0.025, .45, 'Layer 2', transform=ax.transAxes, va='center', - ha='left', - size='7') - - outfig = os.path.join(workspace, 'Figure09_swi2ex4.{0}'.format(fext)) + ax.set_xlabel("Time, in years") + ax.set_ylabel("Elevation, in meters") + ax.text( + 0.025, + 0.55, + "Layer 1", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + ax.text( + 0.025, + 0.45, + "Layer 2", + transform=ax.transAxes, + va="center", + ha="left", + size="7", + ) + + outfig = os.path.join(workspace, "Figure09_swi2ex4.{0}".format(fext)) fig.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) return 0 -if __name__ == '__main__': +if __name__ == "__main__": success = run() diff --git a/examples/scripts/flopy_swi2_ex5.py b/examples/scripts/flopy_swi2_ex5.py index d4c96673f3..96e38edb71 100755 --- a/examples/scripts/flopy_swi2_ex5.py +++ b/examples/scripts/flopy_swi2_ex5.py @@ -12,57 +12,59 @@ import matplotlib.pyplot as plt # --modify default matplotlib settings -updates = {'font.family': ['Univers 57 Condensed', 'Arial'], - 'mathtext.default': 'regular', - 'pdf.compression': 0, - 'pdf.fonttype': 42, - 'legend.fontsize': 7, - 'axes.labelsize': 8, - 'xtick.labelsize': 7, - 'ytick.labelsize': 7} +updates = { + "font.family": ["Univers 57 Condensed", "Arial"], + "mathtext.default": "regular", + "pdf.compression": 0, + "pdf.fonttype": 42, + "legend.fontsize": 7, + "axes.labelsize": 8, + "xtick.labelsize": 7, + "ytick.labelsize": 7, +} plt.rcParams.update(updates) def run(): - workspace = 'swiex5' + workspace = "swiex5" if not os.path.isdir(workspace): os.mkdir(workspace) cleanFiles = False skipRuns = False - fext = 'png' + fext = "png" narg = len(sys.argv) iarg = 0 if narg > 1: while iarg < narg - 1: iarg += 1 basearg = sys.argv[iarg].lower() - if basearg == '--clean': + if basearg == "--clean": cleanFiles = True - elif basearg == '--skipruns': + elif basearg == "--skipruns": skipRuns = True - elif basearg == '--pdf': - fext = 'pdf' + elif basearg == "--pdf": + fext = "pdf" - dirs = [os.path.join(workspace, 'SWI2'), os.path.join(workspace, 'SEAWAT')] + dirs = [os.path.join(workspace, "SWI2"), os.path.join(workspace, "SEAWAT")] if cleanFiles: - print('cleaning all files') - print('excluding *.py files') + print("cleaning all files") + print("excluding *.py files") file_dict = collections.OrderedDict() file_dict[0] = os.listdir(dirs[0]) file_dict[1] = os.listdir(dirs[1]) file_dict[-1] = os.listdir(workspace) for key, files in list(file_dict.items()): - pth = '.' + pth = "." if key >= 0: pth = dirs[key] for f in files: fpth = os.path.join(pth, f) if os.path.isdir(fpth): continue - if '.py' != os.path.splitext(f)[1].lower(): - print(' removing...{}'.format(os.path.basename(f))) + if ".py" != os.path.splitext(f)[1].lower(): + print(" removing...{}".format(os.path.basename(f))) try: os.remove(fpth) except: @@ -82,7 +84,7 @@ def run(): nrow = 1 ncol = 113 delr = np.zeros((ncol), np.float) - delc = 1. + delc = 1.0 r = np.zeros((ncol), np.float) x = np.zeros((ncol), np.float) edge = np.zeros((ncol), np.float) @@ -111,11 +113,11 @@ def run(): for k in range(0, nlay): ibound[k, 0, ncol - 1] = -1 bot = np.zeros((nlay, nrow, ncol), np.float) - dz = 100. / float(nlay - 1) + dz = 100.0 / float(nlay - 1) zall = -np.arange(0, 100 + dz, dz) - zall = np.append(zall, -120.) + zall = np.append(zall, -120.0) tb = -np.arange(dz, 100 + dz, dz) - tb = np.append(tb, -120.) + tb = np.append(tb, -120.0) for k in range(0, nlay): for i in range(0, ncol): bot[k, 0, i] = tb[k] @@ -123,8 +125,8 @@ def run(): isource[:, 0, ncol - 1] = 1 isource[nlay - 1, 0, ncol - 1] = 2 - khb = (0.0000000000256 * 1000. * 9.81 / 0.001) * 60 * 60 * 24 - kvb = (0.0000000000100 * 1000. * 9.81 / 0.001) * 60 * 60 * 24 + khb = (0.0000000000256 * 1000.0 * 9.81 / 0.001) * 60 * 60 * 24 + kvb = (0.0000000000100 * 1000.0 * 9.81 / 0.001) * 60 * 60 * 24 ssb = 1e-5 sszb = 0.2 kh = np.zeros((nlay, nrow, ncol), np.float) @@ -139,14 +141,14 @@ def run(): ss[k, 0, i] = ssb * f ssz[k, 0, i] = sszb * f z = np.ones((nlay), np.float) - z = -100. * z + z = -100.0 * z nwell = 1 for k in range(0, nlay): - if zall[k] > -20. and zall[k + 1] <= -20: + if zall[k] > -20.0 and zall[k + 1] <= -20: nwell = k + 1 - print('nlay={} dz={} nwell={}'.format(nlay, dz, nwell)) - wellQ = -2400. + print("nlay={} dz={} nwell={}".format(nlay, dz, nwell)) + wellQ = -2400.0 wellbtm = -20.0 wellQpm = wellQ / abs(wellbtm) well_data = {} @@ -174,40 +176,68 @@ def run(): for j in range(0, nstp[i]): icnt += 1 if icnt == 365: - ocspd[(i, j)] = ['save head'] + ocspd[(i, j)] = ["save head"] icnt = 0 else: ocspd[(i, j)] = [] - solver2params = {'mxiter': 100, 'iter1': 20, 'npcond': 1, 'zclose': 1.0e-6, - 'rclose': 3e-3, 'relax': 1.0, - 'nbpol': 2, 'damp': 1.0, 'dampt': 1.0} + solver2params = { + "mxiter": 100, + "iter1": 20, + "npcond": 1, + "zclose": 1.0e-6, + "rclose": 3e-3, + "relax": 1.0, + "nbpol": 2, + "damp": 1.0, + "dampt": 1.0, + } # --create model file and run model - modelname = 'swi2ex5' - mf_name = 'mf2005' + modelname = "swi2ex5" + mf_name = "mf2005" if not skipRuns: - ml = flopy.modflow.Modflow(modelname, version='mf2005', - exe_name=mf_name, - model_ws=dirs[0]) - discret = flopy.modflow.ModflowDis(ml, nrow=nrow, ncol=ncol, nlay=nlay, - delr=delr, delc=delc, top=0, - botm=bot, - laycbd=0, nper=nper, perlen=perlen, - nstp=nstp, steady=steady) + ml = flopy.modflow.Modflow( + modelname, version="mf2005", exe_name=mf_name, model_ws=dirs[0] + ) + discret = flopy.modflow.ModflowDis( + ml, + nrow=nrow, + ncol=ncol, + nlay=nlay, + delr=delr, + delc=delc, + top=0, + botm=bot, + laycbd=0, + nper=nper, + perlen=perlen, + nstp=nstp, + steady=steady, + ) bas = flopy.modflow.ModflowBas(ml, ibound=ibound, strt=ihead) - lpf = flopy.modflow.ModflowLpf(ml, hk=kh, vka=kv, ss=ss, sy=ssz, - vkcb=0, - laytyp=0, layavg=1) + lpf = flopy.modflow.ModflowLpf( + ml, hk=kh, vka=kv, ss=ss, sy=ssz, vkcb=0, laytyp=0, layavg=1 + ) wel = flopy.modflow.ModflowWel(ml, stress_period_data=well_data) - swi = flopy.modflow.ModflowSwi2(ml, iswizt=55, npln=1, istrat=1, - toeslope=0.025, tipslope=0.025, - nu=[0, 0.025], zeta=z, ssz=ssz, - isource=isource, nsolver=2, - solver2params=solver2params) + swi = flopy.modflow.ModflowSwi2( + ml, + iswizt=55, + npln=1, + istrat=1, + toeslope=0.025, + tipslope=0.025, + nu=[0, 0.025], + zeta=z, + ssz=ssz, + isource=isource, + nsolver=2, + solver2params=solver2params, + ) oc = flopy.modflow.ModflowOc(ml, stress_period_data=ocspd) - pcg = flopy.modflow.ModflowPcg(ml, hclose=1.0e-6, rclose=3.0e-3, - mxiter=100, iter1=50) + pcg = flopy.modflow.ModflowPcg( + ml, hclose=1.0e-6, rclose=3.0e-3, mxiter=100, iter1=50 + ) # --write the modflow files ml.write_input() m = ml.run_model(silent=False) @@ -216,53 +246,53 @@ def run(): get_stp = [364, 729, 1094, 1459, 364, 729, 1094, 1459] get_per = [0, 0, 0, 0, 1, 1, 1, 1] nswi_times = len(get_per) - zetafile = os.path.join(dirs[0], '{}.zta'.format(modelname)) + zetafile = os.path.join(dirs[0], "{}.zta".format(modelname)) zobj = flopy.utils.CellBudgetFile(zetafile) zeta = [] for kk in zip(get_stp, get_per): - zeta.append(zobj.get_data(kstpkper=kk, text='ZETASRF 1')[0]) + zeta.append(zobj.get_data(kstpkper=kk, text="ZETASRF 1")[0]) zeta = np.array(zeta) # --seawat input - redefine input data that differ from SWI2 nlay_swt = 120 # --mt3d print times - timprs = (np.arange(8) + 1) * 365. + timprs = (np.arange(8) + 1) * 365.0 nprs = len(timprs) # -- ndecay = 4 - ibound = np.ones((nlay_swt, nrow, ncol), 'int') + ibound = np.ones((nlay_swt, nrow, ncol), "int") for k in range(0, nlay_swt): ibound[k, 0, ncol - 1] = -1 bot = np.zeros((nlay_swt, nrow, ncol), np.float) - zall = [0, -20., -40., -60., -80., -100., -120.] - dz = 120. / nlay_swt + zall = [0, -20.0, -40.0, -60.0, -80.0, -100.0, -120.0] + dz = 120.0 / nlay_swt tb = np.arange(nlay_swt) * -dz - dz sconc = np.zeros((nlay_swt, nrow, ncol), np.float) icbund = np.ones((nlay_swt, nrow, ncol), np.int) strt = np.zeros((nlay_swt, nrow, ncol), np.float) - pressure = 0. + pressure = 0.0 g = 9.81 - z = - dz / 2. # cell center + z = -dz / 2.0 # cell center for k in range(0, nlay_swt): for i in range(0, ncol): bot[k, 0, i] = tb[k] - if bot[k, 0, 0] >= -100.: - sconc[k, 0, :] = 0. / 3. * .025 * 1000. / .7143 + if bot[k, 0, 0] >= -100.0: + sconc[k, 0, :] = 0.0 / 3.0 * 0.025 * 1000.0 / 0.7143 else: - sconc[k, 0, :] = 3. / 3. * .025 * 1000. / .7143 + sconc[k, 0, :] = 3.0 / 3.0 * 0.025 * 1000.0 / 0.7143 icbund[k, 0, -1] = -1 - dense = 1000. + 0.7143 * sconc[k, 0, 0] + dense = 1000.0 + 0.7143 * sconc[k, 0, 0] pressure += 0.5 * dz * dense * g if k > 0: z = z - dz - denseup = 1000. + 0.7143 * sconc[k - 1, 0, 0] + denseup = 1000.0 + 0.7143 * sconc[k - 1, 0, 0] pressure += 0.5 * dz * denseup * g strt[k, 0, :] = z + pressure / dense / g # print z, pressure, strt[k, 0, 0], sconc[k, 0, 0] - khb = (0.0000000000256 * 1000. * 9.81 / 0.001) * 60 * 60 * 24 - kvb = (0.0000000000100 * 1000. * 9.81 / 0.001) * 60 * 60 * 24 + khb = (0.0000000000256 * 1000.0 * 9.81 / 0.001) * 60 * 60 * 24 + kvb = (0.0000000000100 * 1000.0 * 9.81 / 0.001) * 60 * 60 * 24 ssb = 1e-5 sszb = 0.2 kh = np.zeros((nlay_swt, nrow, ncol), np.float) @@ -280,12 +310,12 @@ def run(): itype = flopy.mt3d.Mt3dSsm.itype_dict() nwell = 1 for k in range(0, nlay_swt): - if bot[k, 0, 0] >= -20.: + if bot[k, 0, 0] >= -20.0: nwell = k + 1 - print('nlay_swt={} dz={} nwell={}'.format(nlay_swt, dz, nwell)) + print("nlay_swt={} dz={} nwell={}".format(nlay_swt, dz, nwell)) well_data = {} ssm_data = {} - wellQ = -2400. + wellQ = -2400.0 wellbtm = -20.0 wellQpm = wellQ / abs(wellbtm) for ip in range(0, nper): @@ -300,59 +330,90 @@ def run(): welllist[iw, 1] = 0 welllist[iw, 2] = 0 welllist[iw, 3] = q - ssmlist.append([iw, 0, 0, 0., itype['WEL']]) + ssmlist.append([iw, 0, 0, 0.0, itype["WEL"]]) well_data[ip] = welllist.copy() ssm_data[ip] = ssmlist # Define model name for SEAWAT model - modelname = 'swi2ex5_swt' - swtexe_name = 'swtv4' + modelname = "swi2ex5_swt" + swtexe_name = "swtv4" # Create the MODFLOW model data if not skipRuns: - m = flopy.seawat.Seawat(modelname, exe_name=swtexe_name, - model_ws=dirs[1]) - discret = flopy.modflow.ModflowDis(m, nrow=nrow, ncol=ncol, - nlay=nlay_swt, - delr=delr, delc=delc, top=0, - botm=bot, - laycbd=0, nper=nper, perlen=perlen, - nstp=nstp, steady=True) + m = flopy.seawat.Seawat( + modelname, exe_name=swtexe_name, model_ws=dirs[1] + ) + discret = flopy.modflow.ModflowDis( + m, + nrow=nrow, + ncol=ncol, + nlay=nlay_swt, + delr=delr, + delc=delc, + top=0, + botm=bot, + laycbd=0, + nper=nper, + perlen=perlen, + nstp=nstp, + steady=True, + ) bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt) - lpf = flopy.modflow.ModflowLpf(m, hk=kh, vka=kv, ss=ss, sy=ssz, vkcb=0, - laytyp=0, layavg=1) + lpf = flopy.modflow.ModflowLpf( + m, hk=kh, vka=kv, ss=ss, sy=ssz, vkcb=0, laytyp=0, layavg=1 + ) wel = flopy.modflow.ModflowWel(m, stress_period_data=well_data) - oc = flopy.modflow.ModflowOc(m, save_every=365, - save_types=['save head']) - pcg = flopy.modflow.ModflowPcg(m, hclose=1.0e-5, rclose=3.0e-3, - mxiter=100, - iter1=50) + oc = flopy.modflow.ModflowOc( + m, save_every=365, save_types=["save head"] + ) + pcg = flopy.modflow.ModflowPcg( + m, hclose=1.0e-5, rclose=3.0e-3, mxiter=100, iter1=50 + ) # Create the basic MT3DMS model data - adv = flopy.mt3d.Mt3dAdv(m, mixelm=-1, - percel=0.5, - nadvfd=0, - # 0 or 1 is upstream; 2 is central in space - # particle based methods - nplane=4, - mxpart=1e7, - itrack=2, - dceps=1e-4, - npl=16, - nph=16, - npmin=8, - npmax=256) - btn = flopy.mt3d.Mt3dBtn(m, icbund=icbund, prsity=ssz, ncomp=1, - sconc=sconc, - ifmtcn=-1, - chkmas=False, nprobs=10, nprmas=10, dt0=1.0, - ttsmult=1.0, - nprs=nprs, timprs=timprs, mxstrn=1e8) - dsp = flopy.mt3d.Mt3dDsp(m, al=0., trpt=1., trpv=1., dmcoef=0.) + adv = flopy.mt3d.Mt3dAdv( + m, + mixelm=-1, + percel=0.5, + nadvfd=0, + # 0 or 1 is upstream; 2 is central in space + # particle based methods + nplane=4, + mxpart=1e7, + itrack=2, + dceps=1e-4, + npl=16, + nph=16, + npmin=8, + npmax=256, + ) + btn = flopy.mt3d.Mt3dBtn( + m, + icbund=icbund, + prsity=ssz, + ncomp=1, + sconc=sconc, + ifmtcn=-1, + chkmas=False, + nprobs=10, + nprmas=10, + dt0=1.0, + ttsmult=1.0, + nprs=nprs, + timprs=timprs, + mxstrn=1e8, + ) + dsp = flopy.mt3d.Mt3dDsp(m, al=0.0, trpt=1.0, trpv=1.0, dmcoef=0.0) gcg = flopy.mt3d.Mt3dGcg(m, mxiter=1, iter1=50, isolve=1, cclose=1e-7) ssm = flopy.mt3d.Mt3dSsm(m, stress_period_data=ssm_data) # Create the SEAWAT model data - vdf = flopy.seawat.SeawatVdf(m, iwtable=0, densemin=0, densemax=0, - denseref=1000., denseslp=0.7143, - firstdt=1e-3) + vdf = flopy.seawat.SeawatVdf( + m, + iwtable=0, + densemin=0, + densemax=0, + denseref=1000.0, + denseslp=0.7143, + firstdt=1e-3, + ) # write seawat files m.write_input() @@ -361,7 +422,7 @@ def run(): # plot the results # read seawat model data - ucnfile = os.path.join(dirs[1], 'MT3D001.UCN') + ucnfile = os.path.join(dirs[1], "MT3D001.UCN") uobj = flopy.utils.UcnFile(ucnfile) times = uobj.get_times() print(times) @@ -375,16 +436,16 @@ def run(): # spatial data # swi2 bot = np.zeros((1, ncol, nlay), np.float) - dz = 100. / float(nlay - 1) + dz = 100.0 / float(nlay - 1) zall = -np.arange(0, 100 + dz, dz) - zall = np.append(zall, -120.) + zall = np.append(zall, -120.0) tb = -np.arange(dz, 100 + dz, dz) - tb = np.append(tb, -120.) + tb = np.append(tb, -120.0) for k in range(0, nlay): for i in range(0, ncol): bot[0, i, k] = tb[k] # seawat - swt_dz = 120. / nlay_swt + swt_dz = 120.0 / nlay_swt swt_tb = np.zeros((nlay_swt), np.float) zc = -swt_dz / 2.0 for klay in range(0, nlay_swt): @@ -398,32 +459,73 @@ def run(): eps = 1.0e-3 - lc = ['r', 'c', 'g', 'b', 'k'] - cfig = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H'] + lc = ["r", "c", "g", "b", "k"] + cfig = ["A", "B", "C", "D", "E", "F", "G", "H"] inc = 1.0e-3 - xsf = plt.figure(figsize=(fwid, fhgt), facecolor='w') - xsf.subplots_adjust(wspace=0.25, hspace=0.25, left=flft, right=frgt, - bottom=fbot, top=ftop) + xsf = plt.figure(figsize=(fwid, fhgt), facecolor="w") + xsf.subplots_adjust( + wspace=0.25, hspace=0.25, left=flft, right=frgt, bottom=fbot, top=ftop + ) # withdrawal and recovery titles ax = xsf.add_subplot(4, 2, 1) - ax.text(0.0, 1.03, 'Withdrawal', transform=ax.transAxes, va='bottom', - ha='left', size='8') + ax.text( + 0.0, + 1.03, + "Withdrawal", + transform=ax.transAxes, + va="bottom", + ha="left", + size="8", + ) ax = xsf.add_subplot(4, 2, 2) - ax.text(0.0, 1.03, 'Recovery', transform=ax.transAxes, va='bottom', - ha='left', - size='8') + ax.text( + 0.0, + 1.03, + "Recovery", + transform=ax.transAxes, + va="bottom", + ha="left", + size="8", + ) # dummy items for legend ax = xsf.add_subplot(4, 2, 1) - ax.plot([-1, -1], [-1, -1], 'bo', markersize=3, markeredgecolor='blue', - markerfacecolor='None', label='SWI2 interface') - ax.plot([-1, -1], [-1, -1], color='k', linewidth=0.75, linestyle='solid', - label='SEAWAT 50% seawater') - ax.plot([-1, -1], [-1, -1], marker='s', color='k', linewidth=0, - linestyle='none', markeredgecolor='w', - markerfacecolor='0.75', label='SEAWAT 5-95% seawater') - leg = ax.legend(loc='upper left', numpoints=1, ncol=1, labelspacing=0.5, - borderaxespad=1, handlelength=3) + ax.plot( + [-1, -1], + [-1, -1], + "bo", + markersize=3, + markeredgecolor="blue", + markerfacecolor="None", + label="SWI2 interface", + ) + ax.plot( + [-1, -1], + [-1, -1], + color="k", + linewidth=0.75, + linestyle="solid", + label="SEAWAT 50% seawater", + ) + ax.plot( + [-1, -1], + [-1, -1], + marker="s", + color="k", + linewidth=0, + linestyle="none", + markeredgecolor="w", + markerfacecolor="0.75", + label="SEAWAT 5-95% seawater", + ) + leg = ax.legend( + loc="upper left", + numpoints=1, + ncol=1, + labelspacing=0.5, + borderaxespad=1, + handlelength=3, + ) leg._drawFrame = False # data items for itime in range(0, nswi_times): @@ -432,8 +534,8 @@ def run(): for icol in range(0, ncol): for klay in range(0, nlay): # top and bottom of layer - ztop = float('{0:10.3e}'.format(zall[klay])) - zbot = float('{0:10.3e}'.format(zall[klay + 1])) + ztop = float("{0:10.3e}".format(zall[klay])) + zbot = float("{0:10.3e}".format(zall[klay + 1])) # fresh-salt zeta surface zt = zeta[itime, klay, 0, icol] if (ztop - zt) > eps: @@ -447,29 +549,56 @@ def run(): isp = (ic * 2) + 2 ax = xsf.add_subplot(4, 2, isp) # figure title - ax.text(-0.15, 1.025, cfig[itime], transform=ax.transAxes, va='center', - ha='center', size='8') + ax.text( + -0.15, + 1.025, + cfig[itime], + transform=ax.transAxes, + va="center", + ha="center", + size="8", + ) # swi2 - ax.plot(x, zs, 'bo', markersize=3, markeredgecolor='blue', - markerfacecolor='None', label='_None') + ax.plot( + x, + zs, + "bo", + markersize=3, + markeredgecolor="blue", + markerfacecolor="None", + label="_None", + ) # seawat - sc = ax.contour(X, Z, conc[itime, :, :], levels=[17.5], colors='k', - linestyles='solid', linewidths=0.75, zorder=30) - cc = ax.contourf(X, Z, conc[itime, :, :], levels=[0.0, 1.75, 33.250], - colors=['w', '0.75', 'w']) + sc = ax.contour( + X, + Z, + conc[itime, :, :], + levels=[17.5], + colors="k", + linestyles="solid", + linewidths=0.75, + zorder=30, + ) + cc = ax.contourf( + X, + Z, + conc[itime, :, :], + levels=[0.0, 1.75, 33.250], + colors=["w", "0.75", "w"], + ) # set graph limits ax.set_xlim(0, 500) ax.set_ylim(-100, -65) if itime < ndecay: - ax.set_ylabel('Elevation, in meters') + ax.set_ylabel("Elevation, in meters") # x labels ax = xsf.add_subplot(4, 2, 7) - ax.set_xlabel('Horizontal distance, in meters') + ax.set_xlabel("Horizontal distance, in meters") ax = xsf.add_subplot(4, 2, 8) - ax.set_xlabel('Horizontal distance, in meters') + ax.set_xlabel("Horizontal distance, in meters") # simulation time titles for itime in range(0, nswi_times): @@ -483,19 +612,25 @@ def run(): ax = xsf.add_subplot(4, 2, isp) iyr = itime + 1 if iyr > 1: - ctxt = '{} years'.format(iyr) + ctxt = "{} years".format(iyr) else: - ctxt = '{} year'.format(iyr) - ax.text(0.95, 0.925, ctxt, transform=ax.transAxes, va='top', - ha='right', - size='8') - - outfig = os.path.join(workspace, 'Figure11_swi2ex5.{0}'.format(fext)) + ctxt = "{} year".format(iyr) + ax.text( + 0.95, + 0.925, + ctxt, + transform=ax.transAxes, + va="top", + ha="right", + size="8", + ) + + outfig = os.path.join(workspace, "Figure11_swi2ex5.{0}".format(fext)) xsf.savefig(outfig, dpi=300) - print('created...', outfig) + print("created...", outfig) return 0 -if __name__ == '__main__': +if __name__ == "__main__": success = run() diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py index 2ea32548d0..013c647069 100644 --- a/flopy/discretization/structuredgrid.py +++ b/flopy/discretization/structuredgrid.py @@ -436,7 +436,7 @@ def xyzcellcenters(self): @property def grid_lines(self): """ - Get the grid lines as a list + Get the grid lines as a list """ # get edges initially in model coordinates diff --git a/flopy/export/metadata.py b/flopy/export/metadata.py index d72e2e4359..1ff420cee9 100644 --- a/flopy/export/metadata.py +++ b/flopy/export/metadata.py @@ -39,8 +39,10 @@ def __init__(self, sciencebase_id, model): self.model = model self.model_grid = model.modelgrid self.model_time = model.modeltime - self.sciencebase_url = "https://www.sciencebase.gov/catalog/item/{}".format( - sciencebase_id + self.sciencebase_url = ( + "https://www.sciencebase.gov/catalog/item/{}".format( + sciencebase_id + ) ) self.sb = self.get_sciencebase_metadata(sciencebase_id) if self.sb is None: diff --git a/flopy/export/netcdf.py b/flopy/export/netcdf.py index f70194cedf..f717f5f388 100644 --- a/flopy/export/netcdf.py +++ b/flopy/export/netcdf.py @@ -85,14 +85,18 @@ def log(self, phrase): + "\n" ) if self.echo: - print(s,) + print( + s, + ) if self.filename: self.f.write(s) self.items.pop(phrase) else: s = str(t) + " starting: " + str(phrase) + "\n" if self.echo: - print(s,) + print( + s, + ) if self.filename: self.f.write(s) self.items[phrase] = copy.deepcopy(t) @@ -109,7 +113,9 @@ def warn(self, message): """ s = str(datetime.now()) + " WARNING: " + message + "\n" if self.echo: - print(s,) + print( + s, + ) if self.filename: self.f.write(s) return @@ -624,8 +630,7 @@ def difference( self.log("processing variable {0}".format(vname)) def _dt_str(self, dt): - """ for datetime to string for year < 1900 - """ + """for datetime to string for year < 1900""" dt_str = "{0:04d}-{1:02d}-{2:02d}T{3:02d}:{4:02d}:{5:02}Z".format( dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second ) @@ -655,7 +660,7 @@ def write(self): def _initialize_attributes(self): """private method to initial the attributes - of the NetCdf instance + of the NetCdf instance """ assert ( "nc" not in self.__dict__.keys() @@ -707,8 +712,8 @@ def _initialize_attributes(self): self.nc = None def initialize_geometry(self): - """ initialize the geometric information - needed for the netcdf file + """initialize the geometric information + needed for the netcdf file """ try: import pyproj @@ -1354,7 +1359,7 @@ def create_variable( return var def add_global_attributes(self, attr_dict): - """ add global attribute to an initialized file + """add global attribute to an initialized file Parameters ---------- diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py index baae96515c..c1998dcdf5 100755 --- a/flopy/export/shapefile_utils.py +++ b/flopy/export/shapefile_utils.py @@ -418,8 +418,7 @@ def enforce_10ch_limit(names): def get_pyshp_field_info(dtypename): - """Get pyshp dtype information for a given numpy dtype. - """ + """Get pyshp dtype information for a given numpy dtype.""" fields = { "int": ("N", 18, 0), ">> data = MFOutputRequester(mfdict, path, key) >>> data.querybinarydata - """ + """ def __init__(self, mfdict, path, key): self.path = path diff --git a/flopy/mf6/utils/mfobservation.py b/flopy/mf6/utils/mfobservation.py index 2e856ebcf3..ecec8b3e0c 100644 --- a/flopy/mf6/utils/mfobservation.py +++ b/flopy/mf6/utils/mfobservation.py @@ -62,7 +62,7 @@ class Observations: idx = (int), (slice(start, stop)) integer or slice of data to be returned. corresponds to kstp*kper - 1 totim = (float) model time value to return data from - + list_records(): prints a list of all valid record names contained within the Obs.out file get_times(): (list) returns list of time values contained in Obs.out diff --git a/flopy/mf6/utils/reference.py b/flopy/mf6/utils/reference.py index 91d400519c..9e31d9baf7 100644 --- a/flopy/mf6/utils/reference.py +++ b/flopy/mf6/utils/reference.py @@ -66,7 +66,7 @@ class StructuredSpatialReference(object): xul and yul can be explicitly (re)set after SpatialReference instantiation, but only before any of the other attributes and methods are accessed - + """ def __init__( @@ -259,7 +259,7 @@ def attribute_dict(self): def set_spatialreference(self, xul=None, yul=None, rotation=0.0): """ - set spatial reference - can be called from model instance + set spatial reference - can be called from model instance """ # Set origin and rotation @@ -405,7 +405,7 @@ def get_extent(self): def get_grid_lines(self): """ - get the grid lines as a list + get the grid lines as a list """ xmin = self.xedge[0] xmax = self.xedge[-1] @@ -487,8 +487,7 @@ def get_yedge_array(self): return yedge def write_gridSpec(self, filename): - """ write a PEST-style grid specification file - """ + """write a PEST-style grid specification file""" f = open(filename, "w") f.write( "{0:10d} {1:10d}\n".format(self.delc.shape[0], self.delr.shape[0]) @@ -722,9 +721,9 @@ def __setattr__(self, key, value): def set_spatialreference(self, xadj=0.0, yadj=0.0, rotation=0.0): """ - set spatial reference - can be called from model instance - xadj, yadj should be named xadj, yadj since they represent an - adjustment factor + set spatial reference - can be called from model instance + xadj, yadj should be named xadj, yadj since they represent an + adjustment factor """ # Set origin and rotation diff --git a/flopy/modflow/mfbct.py b/flopy/modflow/mfbct.py index 038cfd00a5..428d3d1767 100644 --- a/flopy/modflow/mfbct.py +++ b/flopy/modflow/mfbct.py @@ -65,7 +65,11 @@ def __init__( self.izod = izod self.ifod = ifod self.icbund = Util3d( - model, (nlay, nrow, ncol), np.float32, icbund, "icbund", + model, + (nlay, nrow, ncol), + np.float32, + icbund, + "icbund", ) self.porosity = Util3d( model, (nlay, nrow, ncol), np.float32, porosity, "porosity" @@ -77,7 +81,11 @@ def __init__( self.dth = Util3d(model, (nlay, nrow, ncol), np.float32, dth, "dth") self.dtv = Util3d(model, (nlay, nrow, ncol), np.float32, dth, "dtv") self.sconc = Util3d( - model, (nlay, nrow, ncol), np.float32, sconc, "sconc", + model, + (nlay, nrow, ncol), + np.float32, + sconc, + "sconc", ) self.parent.add_package(self) return diff --git a/flopy/modflow/mfdis.py b/flopy/modflow/mfdis.py index 3d1b85c2a8..210be95a23 100644 --- a/flopy/modflow/mfdis.py +++ b/flopy/modflow/mfdis.py @@ -599,16 +599,16 @@ def get_node(self, lrc_list): def get_layer(self, i, j, elev): """Return the layer for an elevation at an i, j location. - Parameters - ---------- - i : row index (zero-based) - j : column index - elev : elevation (in same units as model) - - Returns - ------- - k : zero-based layer index - """ + Parameters + ---------- + i : row index (zero-based) + j : column index + elev : elevation (in same units as model) + + Returns + ------- + k : zero-based layer index + """ return get_layer(self, i, j, elev) def gettop(self): diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py index fb9fb05a05..7bd24930d7 100644 --- a/flopy/modflow/mfsfr2.py +++ b/flopy/modflow/mfsfr2.py @@ -3156,11 +3156,13 @@ def elevations(self, min_strtop=-10, max_strtop=15000): def slope(self, minimum_slope=1e-4, maximum_slope=1.0): """Checks that streambed slopes are greater than or equal to a specified minimum value. - Low slope values can cause "backup" or unrealistic stream stages with icalc options - where stage is computed. - """ - headertxt = "Checking for streambed slopes of less than {}...\n".format( - minimum_slope + Low slope values can cause "backup" or unrealistic stream stages with icalc options + where stage is computed. + """ + headertxt = ( + "Checking for streambed slopes of less than {}...\n".format( + minimum_slope + ) ) txt = "" if self.verbose: @@ -3191,8 +3193,10 @@ def slope(self, minimum_slope=1e-4, maximum_slope=1.0): passed = True self._txt_footer(headertxt, txt, "minimum slope", passed) - headertxt = "Checking for streambed slopes of greater than {}...\n".format( - maximum_slope + headertxt = ( + "Checking for streambed slopes of greater than {}...\n".format( + maximum_slope + ) ) txt = "" if self.verbose: diff --git a/flopy/modpath/mp7.py b/flopy/modpath/mp7.py index 13d16090bc..e4df040694 100644 --- a/flopy/modpath/mp7.py +++ b/flopy/modpath/mp7.py @@ -171,7 +171,9 @@ def __init__( shape = (nlay, ncpl) elif dis.package_name.lower() == "disu": nodes = dis.nodes.array - shape = tuple(nodes,) + shape = tuple( + nodes, + ) else: msg = ( "DIS, DISV, or DISU packages must be " diff --git a/flopy/modpath/mpsim.py b/flopy/modpath/mpsim.py index 7fc958b70c..aa4a7a4140 100644 --- a/flopy/modpath/mpsim.py +++ b/flopy/modpath/mpsim.py @@ -444,7 +444,7 @@ def __init__(self, model, inputstyle=1, extension="loc", verbose=False): @staticmethod def get_dtypes(): """ - Build numpy dtype for the MODPATH 6 starting locations file. + Build numpy dtype for the MODPATH 6 starting locations file. """ dtype = np.dtype( [ diff --git a/flopy/mt3d/mt.py b/flopy/mt3d/mt.py index 663f4acfba..c5e2469ae3 100644 --- a/flopy/mt3d/mt.py +++ b/flopy/mt3d/mt.py @@ -936,8 +936,10 @@ def load_obs(fname): msg = "First line in file must be \n{}\nFound {}".format( firstline, line.strip() ) - msg += "\n{} does not appear to be a valid MT3D OBS file".format( - fname + msg += ( + "\n{} does not appear to be a valid MT3D OBS file".format( + fname + ) ) raise Exception(msg) diff --git a/flopy/mt3d/mtcts.py b/flopy/mt3d/mtcts.py index 440f8fa462..aa49ab8988 100644 --- a/flopy/mt3d/mtcts.py +++ b/flopy/mt3d/mtcts.py @@ -140,7 +140,9 @@ class Mt3dCts(Package): """ - def __init__(self,): + def __init__( + self, + ): raise NotImplementedError() # # unit number # if unitnumber is None: diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py index ba5f5700e5..82e6eda1da 100755 --- a/flopy/utils/datafile.py +++ b/flopy/utils/datafile.py @@ -206,7 +206,9 @@ def __init__(self, filename, precision, verbose, kwargs): if self.mg is None: self.mg = StructuredGrid( delc=np.ones((self.nrow,)), - delr=np.ones(self.ncol,), + delr=np.ones( + self.ncol, + ), xoff=0.0, yoff=0.0, angrot=0.0, @@ -222,41 +224,41 @@ def to_shapefile( attrib_name="lf_data", ): """ - Export model output data to a shapefile at a specific location - in LayerFile instance. - - Parameters - ---------- - filename : str - Shapefile name to write - kstpkper : tuple of ints - A tuple containing the time step and stress period (kstp, kper). - These are zero-based kstp and kper values. - totim : float - The simulation time. - mflay : integer - MODFLOW zero-based layer number to return. If None, then layer 1 - will be written - attrib_name : str - Base name of attribute columns. (default is 'lf_data') - - Returns - ---------- - None - - See Also - -------- - - Notes - ----- - - Examples - -------- - >>> import flopy - >>> hdobj = flopy.utils.HeadFile('test.hds') - >>> times = hdobj.get_times() - >>> hdobj.to_shapefile('test_heads_sp6.shp', totim=times[-1]) - """ + Export model output data to a shapefile at a specific location + in LayerFile instance. + + Parameters + ---------- + filename : str + Shapefile name to write + kstpkper : tuple of ints + A tuple containing the time step and stress period (kstp, kper). + These are zero-based kstp and kper values. + totim : float + The simulation time. + mflay : integer + MODFLOW zero-based layer number to return. If None, then layer 1 + will be written + attrib_name : str + Base name of attribute columns. (default is 'lf_data') + + Returns + ---------- + None + + See Also + -------- + + Notes + ----- + + Examples + -------- + >>> import flopy + >>> hdobj = flopy.utils.HeadFile('test.hds') + >>> times = hdobj.get_times() + >>> hdobj.to_shapefile('test_heads_sp6.shp', totim=times[-1]) + """ plotarray = np.atleast_3d( self.get_data( diff --git a/flopy/utils/gridintersect.py b/flopy/utils/gridintersect.py index 8dbab97891..a5ceda6a0e 100644 --- a/flopy/utils/gridintersect.py +++ b/flopy/utils/gridintersect.py @@ -1223,12 +1223,16 @@ def _intersect_polygon_structured(self, shp): v_realworld = [] if intersect.geom_type.startswith("Multi"): for ipoly in intersect: - v_realworld += self._transform_geo_interface_polygon( - ipoly + v_realworld += ( + self._transform_geo_interface_polygon( + ipoly + ) ) else: - v_realworld += self._transform_geo_interface_polygon( - intersect + v_realworld += ( + self._transform_geo_interface_polygon( + intersect + ) ) intersect_realworld = rotate( intersect, self.mfgrid.angrot, origin=(0.0, 0.0) diff --git a/flopy/utils/mflistfile.py b/flopy/utils/mflistfile.py index 89edc1cca6..235c60f8c5 100644 --- a/flopy/utils/mflistfile.py +++ b/flopy/utils/mflistfile.py @@ -241,7 +241,7 @@ def get_cumulative(self, names=None): >>> mf_list = MfListBudget("my_model.list") >>> cumulative = mf_list.get_cumulative() - """ + """ if not self._isvalid: return None if names is None: @@ -1005,9 +1005,7 @@ def _parse_time_line(self, line): class SwtListBudget(ListBudget): - """ - - """ + """""" def set_budget_key(self): self.budgetkey = "MASS BUDGET FOR ENTIRE MODEL" @@ -1015,9 +1013,7 @@ def set_budget_key(self): class MfListBudget(ListBudget): - """ - - """ + """""" def set_budget_key(self): self.budgetkey = "VOLUMETRIC BUDGET FOR ENTIRE MODEL" @@ -1025,9 +1021,7 @@ def set_budget_key(self): class Mf6ListBudget(ListBudget): - """ - - """ + """""" def set_budget_key(self): self.budgetkey = "VOLUME BUDGET FOR ENTIRE MODEL" @@ -1035,9 +1029,7 @@ def set_budget_key(self): class MfusgListBudget(ListBudget): - """ - - """ + """""" def set_budget_key(self): self.budgetkey = "VOLUMETRIC BUDGET FOR ENTIRE MODEL" @@ -1045,9 +1037,7 @@ def set_budget_key(self): class SwrListBudget(ListBudget): - """ - - """ + """""" def set_budget_key(self): self.budgetkey = "VOLUMETRIC SURFACE WATER BUDGET FOR ENTIRE MODEL" diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py index 60b504ae48..3e16ed81e0 100644 --- a/flopy/utils/modpathfile.py +++ b/flopy/utils/modpathfile.py @@ -88,7 +88,7 @@ def __init__(self, filename, verbose=False): def _build_index(self): """ - Set position of the start of the pathline data. + Set position of the start of the pathline data. """ self.skiprows = 0 self.file = open(self.fname, "r") @@ -123,7 +123,7 @@ def _build_index(self): def _get_dtypes(self): """ - Build numpy dtype for the MODPATH 6 pathline file. + Build numpy dtype for the MODPATH 6 pathline file. """ if self.version == 3 or self.version == 5: dtype = np.dtype( @@ -706,7 +706,7 @@ def __init__(self, filename, verbose=False): def _build_index(self): """ - Set position of the start of the pathline data. + Set position of the start of the pathline data. """ self.skiprows = 0 self.file = open(self.fname, "r") @@ -750,7 +750,7 @@ def _build_index(self): def _get_dtypes(self): """ - Build numpy dtype for the MODPATH 6 endpoint file. + Build numpy dtype for the MODPATH 6 endpoint file. """ if self.version == 3 or self.version == 5: dtype = self._get_mp35_dtype() @@ -1252,7 +1252,7 @@ def __init__(self, filename, verbose=False): def _build_index(self): """ - Set position of the start of the timeseries data. + Set position of the start of the timeseries data. """ compact = False self.skiprows = 0 @@ -1295,7 +1295,7 @@ def _build_index(self): def _get_dtypes(self): """ - Build numpy dtype for the MODPATH 6 timeseries file. + Build numpy dtype for the MODPATH 6 timeseries file. """ if self.version == 3 or self.version == 5: if self.compact: diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py index 903d986204..c530828719 100755 --- a/flopy/utils/reference.py +++ b/flopy/utils/reference.py @@ -972,8 +972,7 @@ def get_yedge_array(self): return yedge def write_gridSpec(self, filename): - """ write a PEST-style grid specification file - """ + """write a PEST-style grid specification file""" f = open(filename, "w") f.write( "{0:10d} {1:10d}\n".format(self.delc.shape[0], self.delr.shape[0]) diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py index 790490e896..bb42fd0267 100644 --- a/flopy/utils/util_list.py +++ b/flopy/utils/util_list.py @@ -140,7 +140,7 @@ def export(self, f, **kwargs): return export.utils.mflist_export(f, self, **kwargs) def append(self, other): - """ append the recarrays from one MfList to another + """append the recarrays from one MfList to another Parameters ---------- other: variable: an item that can be cast in to an MfList @@ -866,7 +866,7 @@ def __find_last_kper(self, kper): def get_indices(self): """ - a helper function for plotting - get all unique indices + a helper function for plotting - get all unique indices """ names = self.dtype.names lnames = [] @@ -1242,7 +1242,7 @@ def from_4d(cls, model, pak_name, m4ds): @staticmethod def masked4D_arrays_to_stress_period_data(dtype, m4ds): - """ convert a dictionary of 4-dim masked arrays to + """convert a dictionary of 4-dim masked arrays to a stress_period_data style dict of recarray Parameters ---------- diff --git a/flopy/utils/zonbud.py b/flopy/utils/zonbud.py index 5dad800153..133f3aa329 100644 --- a/flopy/utils/zonbud.py +++ b/flopy/utils/zonbud.py @@ -2517,7 +2517,13 @@ def dataframe_to_netcdf_fmt(self, df, flux=True): data[col] = np.zeros((totim.size, zones.size), dtype=float) for i, time in enumerate(totim): - tdf = df.loc[df.totim.isin([time,])] + tdf = df.loc[ + df.totim.isin( + [ + time, + ] + ) + ] tdf = tdf.sort_values(by=["zone"]) for col in df.columns: diff --git a/release/make-release.py b/release/make-release.py index 2212262406..705010f6ff 100644 --- a/release/make-release.py +++ b/release/make-release.py @@ -10,21 +10,27 @@ # update files and paths so that there are the same number of # path and file entries in the paths and files list. Enter '.' # as the path if the file is in the root repository directory -paths = ['../flopy', '../', - '../docs', '../docs', - '../', '../', '../docs'] -files = ['version.py', 'README.md', - 'USGS_release.md', 'PyPi_release.md', - 'code.json', 'DISCLAIMER.md', 'notebook_examples.md'] +paths = ["../flopy", "../", "../docs", "../docs", "../", "../", "../docs"] +files = [ + "version.py", + "README.md", + "USGS_release.md", + "PyPi_release.md", + "code.json", + "DISCLAIMER.md", + "notebook_examples.md", +] # check that there are the same number of entries in files and paths if len(paths) != len(files): - msg = 'The number of entries in paths ' + \ - '({}) must equal '.format(len(paths)) + \ - 'the number of entries in files ({})'.format(len(files)) + msg = ( + "The number of entries in paths " + + "({}) must equal ".format(len(paths)) + + "the number of entries in files ({})".format(len(files)) + ) assert False, msg -pak = 'flopy' +pak = "flopy" # local import of package variables in flopy/version.py # imports author_dict @@ -34,12 +40,12 @@ authors = [] for key in author_dict.keys(): t = key.split() - author = '{}'.format(t[-1]) + author = "{}".format(t[-1]) for str in t[0:-1]: - author += ' {}'.format(str) + author += " {}".format(str) authors.append(author) -approved = '''Disclaimer +approved = """Disclaimer ---------- This software has been approved for release by the U.S. Geological Survey @@ -51,9 +57,9 @@ software is released on condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from its authorized or unauthorized use. -''' +""" -preliminary = '''Disclaimer +preliminary = """Disclaimer ---------- This software is preliminary or provisional and is subject to revision. It is @@ -64,14 +70,14 @@ constitute any such warranty. The software is provided on the condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from the authorized or unauthorized use of the software. -''' +""" def get_disclaimer(): # get current branch branch = get_branch() - if branch.lower().startswith('release') or 'master' in branch.lower(): + if branch.lower().startswith("release") or "master" in branch.lower(): disclaimer = approved is_approved = True else: @@ -86,79 +92,79 @@ def get_branch(): # determine if branch defined on command line for argv in sys.argv: - if 'master' in argv: - branch = 'master' - elif 'develop' in argv.lower(): - branch = 'develop' + if "master" in argv: + branch = "master" + elif "develop" in argv.lower(): + branch = "develop" if branch is None: try: # determine current branch - b = subprocess.Popen(("git", "status"), - stdout=subprocess.PIPE, - stderr=subprocess.STDOUT).communicate()[0] + b = subprocess.Popen( + ("git", "status"), + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + ).communicate()[0] if isinstance(b, bytes): - b = b.decode('utf-8') + b = b.decode("utf-8") for line in b.splitlines(): - if 'On branch' in line: - branch = line.replace('On branch ', '').rstrip() + if "On branch" in line: + branch = line.replace("On branch ", "").rstrip() except: - msg = 'Could not determine current branch. Is git installed?' + msg = "Could not determine current branch. Is git installed?" raise ValueError(msg) return branch def get_version_str(v0, v1, v2): - version_type = ('{}'.format(v0), - '{}'.format(v1), - '{}'.format(v2)) - version = '.'.join(version_type) + version_type = ("{}".format(v0), "{}".format(v1), "{}".format(v2)) + version = ".".join(version_type) return version def get_tag(v0, v1, v2): - tag_type = ('{}'.format(v0), - '{}'.format(v1), - '{}'.format(v2)) - tag = '.'.join(tag_type) + tag_type = ("{}".format(v0), "{}".format(v1), "{}".format(v2)) + tag = ".".join(tag_type) return tag def get_software_citation(version, is_approved): now = datetime.datetime.now() - sb = '' + sb = "" if not is_approved: - sb = ' — release candidate' + sb = " — release candidate" # format author names - line = '[' + line = "[" for ipos, author in enumerate(authors): if ipos > 0: - line += ', ' + line += ", " if ipos == len(authors) - 1: - line += 'and ' + line += "and " sv = author.split() - tauthor = '{}'.format(sv[0]) + tauthor = "{}".format(sv[0]) if len(sv) < 3: gname = sv[1] if len(gname) > 1: - tauthor += ', {}'.format(gname) + tauthor += ", {}".format(gname) else: - tauthor += ', {}.'.format(gname[0]) + tauthor += ", {}.".format(gname[0]) else: - tauthor += ', {}. {}.'.format(sv[1][0], sv[2][0]) + tauthor += ", {}. {}.".format(sv[1][0], sv[2][0]) # add formatted author name to line line += tauthor # add the rest of the citation - line += ', {}, '.format(now.year) + \ - 'FloPy v{}{}: '.format(version, sb) + \ - 'U.S. Geological Survey Software Release, ' + \ - '{}, '.format(now.strftime('%d %B %Y')) + \ - 'http://dx.doi.org/10.5066/F7BK19FH]' + \ - '(http://dx.doi.org/10.5066/F7BK19FH)' + line += ( + ", {}, ".format(now.year) + + "FloPy v{}{}: ".format(version, sb) + + "U.S. Geological Survey Software Release, " + + "{}, ".format(now.strftime("%d %B %Y")) + + "http://dx.doi.org/10.5066/F7BK19FH]" + + "(http://dx.doi.org/10.5066/F7BK19FH)" + ) return line @@ -171,44 +177,49 @@ def update_version(): vmajor = 0 vminor = 0 vmicro = 0 - lines = [line.rstrip('\n') for line in open(fpth, 'r')] + lines = [line.rstrip("\n") for line in open(fpth, "r")] for idx, line in enumerate(lines): t = line.split() - if 'major =' in line: + if "major =" in line: vmajor = int(t[2]) - elif 'minor =' in line: + elif "minor =" in line: vminor = int(t[2]) - elif 'micro =' in line: + elif "micro =" in line: vmicro = int(t[2]) - elif '__version__' in line: + elif "__version__" in line: name_pos = idx + 1 except: - msg = 'There was a problem updating the version file' + msg = "There was a problem updating the version file" raise IOError(msg) try: # write new version file - f = open(fpth, 'w') - f.write('# {} version file automatically '.format(pak) + - 'created using...{0}\n'.format(os.path.basename(__file__))) - f.write('# created on...' + - '{0}\n'.format( - datetime.datetime.now().strftime("%B %d, %Y %H:%M:%S"))) - f.write('\n') - f.write('major = {}\n'.format(vmajor)) - f.write('minor = {}\n'.format(vminor)) - f.write('micro = {}\n'.format(vmicro)) + f = open(fpth, "w") + f.write( + "# {} version file automatically ".format(pak) + + "created using...{0}\n".format(os.path.basename(__file__)) + ) + f.write( + "# created on..." + + "{0}\n".format( + datetime.datetime.now().strftime("%B %d, %Y %H:%M:%S") + ) + ) + f.write("\n") + f.write("major = {}\n".format(vmajor)) + f.write("minor = {}\n".format(vminor)) + f.write("micro = {}\n".format(vmicro)) f.write("__version__ = '{:d}.{:d}.{:d}'.format(major, minor, micro)\n") # write the remainder of the version file if name_pos is not None: for line in lines[name_pos:]: - f.write('{}\n'.format(line)) + f.write("{}\n".format(line)) f.close() - print('Successfully updated version.py') + print("Successfully updated version.py") except: - msg = 'There was a problem updating the version file' + msg = "There was a problem updating the version file" raise IOError(msg) # update README.md with new version information @@ -235,24 +246,24 @@ def update_codejson(vmajor, vminor, vmicro): version = get_tag(vmajor, vminor, vmicro) # load and modify json file - with open(json_fname, 'r') as f: + with open(json_fname, "r") as f: data = json.load(f, object_pairs_hook=OrderedDict) # modify the json file data now = datetime.datetime.now() - sdate = now.strftime('%Y-%m-%d') - data[0]['date']['metadataLastUpdated'] = sdate - if branch.lower().startswith('release') or 'master' in branch.lower(): - data[0]['version'] = version - data[0]['status'] = 'Production' + sdate = now.strftime("%Y-%m-%d") + data[0]["date"]["metadataLastUpdated"] = sdate + if branch.lower().startswith("release") or "master" in branch.lower(): + data[0]["version"] = version + data[0]["status"] = "Production" else: - data[0]['version'] = version - data[0]['status'] = 'Release Candidate' + data[0]["version"] = version + data[0]["status"] = "Release Candidate" # rewrite the json file - with open(json_fname, 'w') as f: + with open(json_fname, "w") as f: json.dump(data, f, indent=4) - f.write('\n') + f.write("\n") return @@ -263,54 +274,61 @@ def update_readme_markdown(vmajor, vminor, vmicro): # define branch if is_approved: - branch = 'master' + branch = "master" else: - branch = 'develop' + branch = "develop" # create version version = get_tag(vmajor, vminor, vmicro) # read README.md into memory fpth = os.path.join(paths[1], files[1]) - with open(fpth, 'r') as file: + with open(fpth, "r") as file: lines = [line.rstrip() for line in file] # rewrite README.md terminate = False - f = open(fpth, 'w') + f = open(fpth, "w") for line in lines: - if '### Version ' in line: - line = '### Version {}'.format(version) + if "### Version " in line: + line = "### Version {}".format(version) if not is_approved: - line += ' — release candidate' - elif '[Build Status]' in line: - line = '[![Build Status](https://travis-ci.org/modflowpy/' + \ - 'flopy.svg?branch={})]'.format(branch) + \ - '(https://travis-ci.org/modflowpy/flopy)' - elif '[Coverage Status]' in line: - line = '[![Coverage Status](https://coveralls.io/repos/github/' + \ - 'modflowpy/flopy/badge.svg?branch={})]'.format(branch) + \ - '(https://coveralls.io/github/modflowpy/' + \ - 'flopy?branch={})'.format(branch) - elif '[Binder]' in line: + line += " — release candidate" + elif "[Build Status]" in line: + line = ( + "[![Build Status](https://travis-ci.org/modflowpy/" + + "flopy.svg?branch={})]".format(branch) + + "(https://travis-ci.org/modflowpy/flopy)" + ) + elif "[Coverage Status]" in line: + line = ( + "[![Coverage Status](https://coveralls.io/repos/github/" + + "modflowpy/flopy/badge.svg?branch={})]".format(branch) + + "(https://coveralls.io/github/modflowpy/" + + "flopy?branch={})".format(branch) + ) + elif "[Binder]" in line: # [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/modflowpy/flopy.git/develop) - line = '[![Binder](https://mybinder.org/badge_logo.svg)]' + \ - '(https://mybinder.org/v2/gh/modflowpy/flopy.git/' + \ - '{}'.format(branch) + ')' - elif 'http://dx.doi.org/10.5066/F7BK19FH' in line: + line = ( + "[![Binder](https://mybinder.org/badge_logo.svg)]" + + "(https://mybinder.org/v2/gh/modflowpy/flopy.git/" + + "{}".format(branch) + + ")" + ) + elif "http://dx.doi.org/10.5066/F7BK19FH" in line: line = get_software_citation(version, is_approved) - elif 'Disclaimer' in line: + elif "Disclaimer" in line: line = disclaimer terminate = True - f.write('{}\n'.format(line)) + f.write("{}\n".format(line)) if terminate: break f.close() # write disclaimer markdown file - fpth = os.path.join(paths[0], 'DISCLAIMER.md') - f = open(fpth, 'w') + fpth = os.path.join(paths[0], "DISCLAIMER.md") + f = open(fpth, "w") f.write(disclaimer) f.close() @@ -323,25 +341,28 @@ def update_notebook_examples_markdown(): # define branch if is_approved: - branch = 'master' + branch = "master" else: - branch = 'develop' + branch = "develop" # read notebook_examples.md into memory fpth = os.path.join(paths[6], files[6]) - with open(fpth, 'r') as file: + with open(fpth, "r") as file: lines = [line.rstrip() for line in file] # rewrite notebook_examples.md terminate = False - f = open(fpth, 'w') + f = open(fpth, "w") for line in lines: - if '[Binder]' in line: + if "[Binder]" in line: # [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/modflowpy/flopy.git/develop) - line = '[![Binder](https://mybinder.org/badge_logo.svg)]' + \ - '(https://mybinder.org/v2/gh/modflowpy/flopy.git/' + \ - '{}'.format(branch) + ')' - f.write('{}\n'.format(line)) + line = ( + "[![Binder](https://mybinder.org/badge_logo.svg)]" + + "(https://mybinder.org/v2/gh/modflowpy/flopy.git/" + + "{}".format(branch) + + ")" + ) + f.write("{}\n".format(line)) f.close() @@ -354,108 +375,110 @@ def update_USGSmarkdown(vmajor, vminor, vmicro): # set repo branch string based on is_approved if is_approved: - repo_str = 'master' + repo_str = "master" else: - repo_str = 'develop' + repo_str = "develop" # create version version = get_tag(vmajor, vminor, vmicro) # read README.md into memory fpth = os.path.join(paths[1], files[1]) - with open(fpth, 'r') as file: + with open(fpth, "r") as file: lines = [line.rstrip() for line in file] # write USGS_release.md fpth = os.path.join(paths[2], files[2]) - f = open(fpth, 'w') + f = open(fpth, "w") # write PyPi_release.md fpth = os.path.join(paths[3], files[3]) - f2 = open(fpth, 'w') + f2 = open(fpth, "w") # date and branch information now = datetime.datetime.now() sdate = now.strftime("%m/%d/%Y") # write header information - f.write('---\n') - f.write('title: FloPy Release Notes\n') - f.write('author:\n') + f.write("---\n") + f.write("title: FloPy Release Notes\n") + f.write("author:\n") for author in authors: sv = author.split() - tauthor = '{}'.format(sv[1]) + tauthor = "{}".format(sv[1]) if len(sv) > 2: - tauthor += ' {}.'.format(sv[2][0]) - tauthor += ' {}'.format(sv[0]) - f.write(' - {}\n'.format(tauthor)) - f.write('header-includes:\n') - f.write(' - \\usepackage{fancyhdr}\n') - f.write(' - \\usepackage{lastpage}\n') - f.write(' - \\pagestyle{fancy}\n') - f.write(' - \\fancyhf{{}}\n') - f.write(' - \\fancyhead[LE, LO, RE, RO]{}\n') - f.write(' - \\fancyhead[CE, CO]{FloPy Release Notes}\n') - f.write(' - \\fancyfoot[LE, RO]{{FloPy version {}}}\n'.format(version)) - f.write(' - \\fancyfoot[CO, CE]{\\thepage\\ of \\pageref{LastPage}}\n') - f.write(' - \\fancyfoot[RE, LO]{{{}}}\n'.format(sdate)) - f.write('geometry: margin=0.75in\n') - f.write('---\n\n') + tauthor += " {}.".format(sv[2][0]) + tauthor += " {}".format(sv[0]) + f.write(" - {}\n".format(tauthor)) + f.write("header-includes:\n") + f.write(" - \\usepackage{fancyhdr}\n") + f.write(" - \\usepackage{lastpage}\n") + f.write(" - \\pagestyle{fancy}\n") + f.write(" - \\fancyhf{{}}\n") + f.write(" - \\fancyhead[LE, LO, RE, RO]{}\n") + f.write(" - \\fancyhead[CE, CO]{FloPy Release Notes}\n") + f.write(" - \\fancyfoot[LE, RO]{{FloPy version {}}}\n".format(version)) + f.write(" - \\fancyfoot[CO, CE]{\\thepage\\ of \\pageref{LastPage}}\n") + f.write(" - \\fancyfoot[RE, LO]{{{}}}\n".format(sdate)) + f.write("geometry: margin=0.75in\n") + f.write("---\n\n") # write select information from README.md writeline = False for line in lines: - if line == 'Introduction': + if line == "Introduction": writeline = True - elif line == 'Installation': + elif line == "Installation": writeline = False - elif line == 'Documentation': + elif line == "Documentation": writeline = True - elif line == 'Getting Started': + elif line == "Getting Started": writeline = False - elif line == 'Contributing': + elif line == "Contributing": writeline = True - elif line == 'How to Cite': + elif line == "How to Cite": writeline = True - elif 'http://dx.doi.org/10.5066/F7BK19FH' in line: + elif "http://dx.doi.org/10.5066/F7BK19FH" in line: writeline = True line = get_software_citation(version, is_approved) - elif line == 'MODFLOW Resources': + elif line == "MODFLOW Resources": writeline = False - elif 'Installing the latest FloPy release candidate' in line: + elif "Installing the latest FloPy release candidate" in line: writeline = False - elif line == 'Disclaimer': + elif line == "Disclaimer": writeline = True - elif '[MODFLOW 6](docs/mf6.md)' in line: - line = line.replace('[MODFLOW 6](docs/mf6.md)', 'MODFLOW 6') + elif "[MODFLOW 6](docs/mf6.md)" in line: + line = line.replace("[MODFLOW 6](docs/mf6.md)", "MODFLOW 6") if writeline: # USGS_release.md - f.write('{}\n'.format(line)) + f.write("{}\n".format(line)) # PyPi_release.md - line = line.replace('***', '*') - line = line.replace('##### ', '') - if 'CONTRIBUTING.md' in line: - rstr = 'https://github.com/modflowpy/flopy/blob/' + \ - '{}/CONTRIBUTING.md'.format(repo_str) - line = line.replace('CONTRIBUTING.md', rstr) - f2.write('{}\n'.format(line)) + line = line.replace("***", "*") + line = line.replace("##### ", "") + if "CONTRIBUTING.md" in line: + rstr = ( + "https://github.com/modflowpy/flopy/blob/" + + "{}/CONTRIBUTING.md".format(repo_str) + ) + line = line.replace("CONTRIBUTING.md", rstr) + f2.write("{}\n".format(line)) # write installation information - cweb = 'https://water.usgs.gov/ogw/flopy/flopy-{}.zip'.format(version) - line = '' - line += 'Installation\n' - line += '-----------------------------------------------\n' - line += 'To install FloPy version {} '.format(version) - line += 'from the USGS FloPy website:\n' - line += '```\n' - line += 'pip install {}\n'.format(cweb) - line += '```\n\n' - line += 'To update to FloPy version {} '.format(version) - line += 'from the USGS FloPy website:\n' - line += '```\n' - line += 'pip install {} --upgrade\n'.format(cweb) - line += '```\n' + cweb = "https://water.usgs.gov/ogw/flopy/flopy-{}.zip".format(version) + line = "" + line += "Installation\n" + line += "-----------------------------------------------\n" + line += "To install FloPy version {} ".format(version) + line += "from the USGS FloPy website:\n" + line += "```\n" + line += "pip install {}\n".format(cweb) + line += "```\n\n" + line += "To update to FloPy version {} ".format(version) + line += "from the USGS FloPy website:\n" + line += "```\n" + line += "pip install {} --upgrade\n".format(cweb) + line += "```\n" # write installation instructions of USGS_release.md f.write(line) @@ -471,4 +494,4 @@ def update_USGSmarkdown(vmajor, vminor, vmicro): if __name__ == "__main__": update_version() - get_software_citation('3.1.1', True) \ No newline at end of file + get_software_citation("3.1.1", True) diff --git a/release/update-version_changes.py b/release/update-version_changes.py index 78eff258b1..0b3b96290b 100644 --- a/release/update-version_changes.py +++ b/release/update-version_changes.py @@ -4,6 +4,7 @@ PY3 = sys.version_info[0] >= 3 + def process_Popen(cmdlist, verbose=False): """Generic function to initialize a Popen process. @@ -34,152 +35,161 @@ def process_Popen(cmdlist, verbose=False): # catch non-zero return code if process.returncode != 0: - errmsg = '{} failed\n'.format(' '.join(process.args)) + \ - '\tstatus code {}\n'.format(process.returncode) + errmsg = "{} failed\n".format( + " ".join(process.args) + ) + "\tstatus code {}\n".format(process.returncode) raise ChildProcessError(errmsg) return stdout + def get_version(): major = 0 minor = 0 micro = 0 # read existing version file - fpth = os.path.join('..', 'flopy', 'version.py') - f = open(fpth, 'r') - lines = [line.rstrip('\n') for line in open(fpth, 'r')] + fpth = os.path.join("..", "flopy", "version.py") + f = open(fpth, "r") + lines = [line.rstrip("\n") for line in open(fpth, "r")] for idx, line in enumerate(lines): t = line.split() - if 'major =' in line: + if "major =" in line: major = int(t[2]) - elif 'minor =' in line: + elif "minor =" in line: minor = int(t[2]) - elif 'micro =' in line: + elif "micro =" in line: micro = int(t[2]) f.close() - return '{:d}.{:d}.{:d}'.format(major, minor, micro) + return "{:d}.{:d}.{:d}".format(major, minor, micro) + def get_branch(): branch = None # determine if branch defined on command line for argv in sys.argv: - if '--master' in argv: - branch = 'master' - elif '--develop' in argv.lower(): - branch = 'develop' + if "--master" in argv: + branch = "master" + elif "--develop" in argv.lower(): + branch = "develop" # use git to determine branch if branch is None: cmdlist = ("git", "status") stdout = process_Popen(cmdlist) for line in stdout.splitlines(): - if 'On branch' in line: - branch = line.replace('On branch ', '').rstrip() + if "On branch" in line: + branch = line.replace("On branch ", "").rstrip() break return branch + def get_last_tag_date(): - cmdlist = ("git", "log", "--tags", "--no-walk", "--date=iso-local", - "--pretty='%cd %D'") + cmdlist = ( + "git", + "log", + "--tags", + "--no-walk", + "--date=iso-local", + "--pretty='%cd %D'", + ) stdout = process_Popen(cmdlist) line = stdout.splitlines()[0] ipos = line.find("tag") if ipos > -1: tag_date = line[1:20] - tag = line[ipos+4:].split(',')[0].strip() + tag = line[ipos + 4 :].split(",")[0].strip() return tag, tag_date + def get_hash_dict(branch): tag, tag_date = get_last_tag_date() - # get hash and fmt = '--pretty="%H"' since = '--since="{}"'.format(tag_date) - hash_dict = {'fix': {}, 'feat': {}} + hash_dict = {"fix": {}, "feat": {}} cmdlist = ("git", "log", branch, fmt, since) stdout = process_Popen(cmdlist) fix_dict = {} feat_dict = {} - fix_tags = ['fix', 'bug'] - feat_tags = ['feat'] + fix_tags = ["fix", "bug"] + feat_tags = ["feat"] # parse stdout for line in stdout.splitlines(): - hash = line.split()[0].replace('"', '') - url = 'https://github.com/modflowpy/flopy/commit/{}'.format(hash) + hash = line.split()[0].replace('"', "") + url = "https://github.com/modflowpy/flopy/commit/{}".format(hash) fmt = '--pretty="%s."' cmdlist = ("git", "log", fmt, "-n1", hash) - subject = process_Popen(cmdlist).strip().replace('"', '') + subject = process_Popen(cmdlist).strip().replace('"', "") fmt = '--pretty="Committed by %aN on %cd."' cmdlist = ("git", "log", "--date=short", fmt, "-n1", hash) - cdate = process_Popen(cmdlist).strip().replace('"', '') - ipos = subject.find(':') + cdate = process_Popen(cmdlist).strip().replace('"', "") + ipos = subject.find(":") key = None if ipos > -1: type = subject[0:ipos] - subject = subject.replace(type + ':', '').strip().capitalize() + subject = subject.replace(type + ":", "").strip().capitalize() for tag in fix_tags: if type.lower().startswith(tag): - key = 'fix' + key = "fix" type = key + type[3:] break if key is None: for tag in feat_tags: if type.lower().startswith(tag): - key = 'feat' + key = "feat" break else: type = None slower = subject.lower() for tag in fix_tags: if slower.startswith(tag): - key = 'fix' - type = 'fix()' + key = "fix" + type = "fix()" break if key is None: for tag in feat_tags: if slower.startswith(tag): - key = 'feat' - type = 'feat()' + key = "feat" + type = "feat()" break if key is not None: - message = '[{}]({}): '.format(type, url) - message += subject + ' ' + cdate - if key == 'fix': + message = "[{}]({}): ".format(type, url) + message += subject + " " + cdate + if key == "fix": fix_dict[hash] = message - elif key == 'feat': + elif key == "feat": feat_dict[hash] = message - # add dictionaries to the hash dictionary - hash_dict['fix'] = fix_dict - hash_dict['feat'] = feat_dict + hash_dict["fix"] = fix_dict + hash_dict["feat"] = feat_dict return hash_dict + def create_md(hash_dict): # get current version information version = get_version() - tag = '### Version' - version_text = '{} {}'.format(tag, version) + tag = "### Version" + version_text = "{} {}".format(tag, version) # # read the lines in the existing version_changes.md - fpth = os.path.join('..', 'docs', 'version_changes.md') + fpth = os.path.join("..", "docs", "version_changes.md") with open(fpth) as f: lines = f.read().splitlines() - # rewrite the existing version_changes.md - f = open(fpth, 'w') + f = open(fpth, "w") write_line = True write_update = True @@ -192,35 +202,35 @@ def create_md(hash_dict): # write the changes for the latest comment if write_update: write_update = False - f.write('{}\n\n'.format(version_text)) + f.write("{}\n\n".format(version_text)) write_version_changes(f, hash_dict) if write_line: - f.write('{}\n'.format(line)) - + f.write("{}\n".format(line)) f.close() + def write_version_changes(f, hash_dict): # features istart = True - for key, val in hash_dict['feat'].items(): + for key, val in hash_dict["feat"].items(): if istart: istart = False - f.write('* New features:\n\n') - f.write(4*' ' + '* ' + val + '\n') + f.write("* New features:\n\n") + f.write(4 * " " + "* " + val + "\n") if not istart: - f.write('\n\n') + f.write("\n\n") # bug fixes istart = True - for key, val in hash_dict['fix'].items(): + for key, val in hash_dict["fix"].items(): if istart: istart = False - f.write('* Bug fixes:\n\n') - f.write(4 * ' ' + '* ' + val + '\n') + f.write("* Bug fixes:\n\n") + f.write(4 * " " + "* " + val + "\n") if not istart: - f.write('\n') + f.write("\n") return @@ -234,4 +244,4 @@ def main(): if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/setup.py b/setup.py index efe3b8a9f3..6283194693 100644 --- a/setup.py +++ b/setup.py @@ -4,7 +4,7 @@ # ensure minimum version of Python is running if sys.version_info[0:2] < (3, 5): - raise RuntimeError('Flopy requires Python >= 3.5') + raise RuntimeError("Flopy requires Python >= 3.5") # local import of package variables in flopy/version.py # imports __version__, __pakname__, __author__, __author_email__ @@ -13,27 +13,42 @@ try: import pypandoc - fpth = os.path.join('docs', 'PyPi_release.md') - long_description = pypandoc.convert(fpth, 'rst') + fpth = os.path.join("docs", "PyPi_release.md") + long_description = pypandoc.convert(fpth, "rst") except ImportError: - long_description = '' + long_description = "" -setup(name=__pakname__, - description='FloPy is a Python package to create, run, and ' + - 'post-process MODFLOW-based models.', - long_description=long_description, - author=__author__, - author_email=__author_email__, - url='https://github.com/modflowpy/flopy/', - license='CC0', - platforms='Windows, Mac OS-X, Linux', - install_requires=['numpy'], - packages=['flopy', 'flopy.modflow', 'flopy.modflowlgr', 'flopy.modpath', - 'flopy.mt3d', 'flopy.seawat', 'flopy.utils', 'flopy.plot', - 'flopy.pest', 'flopy.export', 'flopy.discretization', - 'flopy.mf6', 'flopy.mf6.coordinates', 'flopy.mf6.data', - 'flopy.mf6.modflow', 'flopy.mf6.utils'], - include_package_data=True, # includes files listed in MANIFEST.in - # use this version ID if .svn data cannot be found - version=__version__, - classifiers=['Topic :: Scientific/Engineering :: Hydrology']) +setup( + name=__pakname__, + description="FloPy is a Python package to create, run, and " + + "post-process MODFLOW-based models.", + long_description=long_description, + author=__author__, + author_email=__author_email__, + url="https://github.com/modflowpy/flopy/", + license="CC0", + platforms="Windows, Mac OS-X, Linux", + install_requires=["numpy"], + packages=[ + "flopy", + "flopy.modflow", + "flopy.modflowlgr", + "flopy.modpath", + "flopy.mt3d", + "flopy.seawat", + "flopy.utils", + "flopy.plot", + "flopy.pest", + "flopy.export", + "flopy.discretization", + "flopy.mf6", + "flopy.mf6.coordinates", + "flopy.mf6.data", + "flopy.mf6.modflow", + "flopy.mf6.utils", + ], + include_package_data=True, # includes files listed in MANIFEST.in + # use this version ID if .svn data cannot be found + version=__version__, + classifiers=["Topic :: Scientific/Engineering :: Hydrology"], +)