1+ import flopy4
2+ ws = './mymodel'
3+ name = 'mymodel'
4+ sim = flopy4 .mf6 .MFSimulation (sim_name = name , sim_ws = ws , exe_name = 'mf6' )
5+ tdis = flopy4 .mf6 .ModflowTdis (sim )
6+ ims = flopy4 .mf6 .ModflowIms (sim )
7+ gwf = flopy4 .mf6 .ModflowGwf (sim , modelname = name , save_flows = True )
8+ dis = flopy4 .mf6 .ModflowGwfdis (gwf , nrow = 10 , ncol = 10 )
9+ ic = flopy4 .mf6 .ModflowGwfic (gwf )
10+ npf = flopy4 .mf6 .ModflowGwfnpf (gwf , save_specific_discharge = True )
11+ chd = flopy4 .mf6 .ModflowGwfchd (gwf ) # self.stress_period_data = sparse.DOK((ncol, nrow, nlay))
12+ chd .stress_period_data [0 ,0 ,0 ] = 1.
13+ chd .stress_period_data [0 ,9 ,9 ] = 0.
14+
15+ # chd = flopy4.mf6.ModflowGwfchd(gwf, stress_period_data=sparse.COO([[0,0], [0,9], [0,9]], [1., 0.]))
16+ # chd = flopy4.mf6.ModflowGwfchd(gwf, stress_period_data=[Chd.PeriodData((0, 0, 0), 1.),
17+ # Chd.PeriodData[(0, 9, 9), 0.]])
18+ budget_file = name + '.bud'
19+ head_file = name + '.hds'
20+ oc = flopy4 .mf6 .ModflowGwfoc (gwf ,
21+ budget_filerecord = budget_file ,
22+ head_filerecord = head_file ,
23+ save = {tis .get_period ("10-20-2020" ): {"head" : "ALL" , "budget" : "ALL" }}
24+ print = {"budget" : np .ones ()}
25+ # save_head="ALL",
26+ # save_budget=np.ones(nstp),
27+ # print_head="",
28+ # print_budget="",
29+ oc .saverecord ["HEAD" ]= "ALL"
30+ oc .saverecord ["BUDGET" ]= np .ones ((nstp ))
31+
32+ # Normally you would write and run the simulation
33+ # sim.write_simulation()
34+ # sim.run_simulation()
35+
36+
37+ head = gwf .output .head ().get_data ()
38+ bud = gwf .output .budget ()
39+
40+ spdis = bud .get_data (text = 'DATA-SPDIS' )[0 ]
41+ qx , qy , qz = flopy .utils .postprocessing .get_specific_discharge (spdis , gwf )
42+ pmv = flopy .plot .PlotMapView (gwf )
43+ pmv .plot_array (head )
44+ pmv .plot_grid (colors = 'white' )
45+ pmv .contour_array (head , levels = [.2 , .4 , .6 , .8 ], linewidths = 3. )
46+ pmv .plot_vector (qx , qy , normalize = True , color = "white" )
47+
48+
49+
50+
51+ # from demo:
52+
53+ # # Create a simulation.
54+ #
55+ # sim = Simulation()
56+ # tdis = Tdis(sim=sim, nper=1, perioddata=[Tdis.PeriodData()])
57+ # gwf = Gwf(sim=sim)
58+ # dis = Dis(model=gwf)
59+ # ic = Ic(model=gwf, strt=1.0)
60+ # oc = Oc(model=gwf, perioddata=[Oc.Steps()])
61+ # npf = Npf(model=gwf, icelltype=0, k=1.0)
62+ #
63+ # # View the data tree.
64+ # sim.data
0 commit comments