|
1 | | -#!/usr/bin/env python3 |
2 | 1 | # -*- coding: utf-8 -*- |
3 | 2 | """ |
4 | 3 | Subsurface flow and tranport simulation in a random binary inclusion structure |
5 | 4 | with the FEM solver of the groundwater flow equation and advection-dispersion- |
6 | 5 | equation OpenGeoSys. |
7 | 6 |
|
8 | | -This script allows to perform a 2D simulation according to the flow and transport |
9 | | -situation of the MADE-1 tracer tests (Columbus, Mississippi) in a 2D setting |
10 | | -with a binary inclusion structure of hydraulic conductivity. |
| 7 | +This script allows to perform a 2D simulation according to the flow and transport |
| 8 | +situation of the MADE-1 tracer tests (Columbus, Mississippi) in a 2D setting |
| 9 | +with a binary inclusion structure of hydraulic conductivity. |
11 | 10 |
|
12 | | -Simulation are performed with the FEM solver OpenGeoSys (OGS5) for subsurface |
13 | | -flow and transport controlled via the API ogs5py: Initialization of the project, |
14 | | -writing files, simulation run and reading of simulation results. Setting-addapted |
| 11 | +Simulation are performed with the FEM solver OpenGeoSys (OGS5) for subsurface |
| 12 | +flow and transport controlled via the API ogs5py: Initialization of the project, |
| 13 | +writing files, simulation run and reading of simulation results. Setting-addapted |
15 | 14 | conductivity fields are generated with the script binary_inclusions. |
16 | | - |
| 15 | +
|
17 | 16 | @author: A. Zech |
18 | 17 | Licence MIT, A.Zech, 2020 |
19 | 18 | """ |
|
115 | 114 | sim.gli.add_polyline(points=[2, 3], name="top") |
116 | 115 | sim.gli.add_polyline(points=[1, 2], name="right") |
117 | 116 | sim.gli.add_polyline( |
118 | | - points=[[0, domain["z_0"] + 5.2, 0], [0, domain["z_0"] + 5.8, 0]], |
| 117 | + points=[[0, domain["z_0"] + 5.2, 0], [0, domain["z_0"] + 5.8, 0]], |
119 | 118 | name="source" |
120 | 119 | ) |
121 | 120 | sim.gli.swap_axis("y", "z") |
122 | 121 |
|
123 | 122 | ### ----------------------- properties------------------------------------- # |
124 | 123 | sim.mfp.add_block( # FLUID_PROPERTIES |
125 | | - FLUID_TYPE="LIQUID", |
126 | | - PCS_TYPE="HEAD", |
127 | | - DENSITY=[1, 1000.0], |
| 124 | + FLUID_TYPE="LIQUID", |
| 125 | + PCS_TYPE="HEAD", |
| 126 | + DENSITY=[1, 1000.0], |
128 | 127 | VISCOSITY=[1, 0.001] |
129 | 128 | ) |
130 | 129 | sim.mcp.add_block( |
131 | | - NAME="CONCENTRATION1", |
132 | | - MOBILE=1, |
| 130 | + NAME="CONCENTRATION1", |
| 131 | + MOBILE=1, |
133 | 132 | DIFFUSION=[1, 1.0e-08] |
134 | 133 | ) |
135 | 134 | sim.msp.add_block(DENSITY=[1, 2000]) |
|
142 | 141 | STORAGE="1 {}".format(storage), |
143 | 142 | MASS_DISPERSION=[1, dispersivity_long, dispersivity_trans], |
144 | 143 | ) |
145 | | -# sim.mmp.update_block(PERMEABILITY_TENSOR=['ISOTROPIC', '1e-6']) # for homogeneous K-distribution |
146 | 144 |
|
147 | 145 | ### --------------------heterogeneous K: binary structure ------------------ # |
148 | 146 |
|
149 | 147 | sim.mpd.add(name="conductivity") |
150 | 148 | sim.mpd.add_block( |
151 | | - MSH_TYPE="GROUNDWATER_FLOW", |
152 | | - MMP_TYPE="PERMEABILITY", |
| 149 | + MSH_TYPE="GROUNDWATER_FLOW", |
| 150 | + MMP_TYPE="PERMEABILITY", |
153 | 151 | DIS_TYPE="ELEMENT" |
154 | 152 | ) |
155 | 153 |
|
|
269 | 267 | print("Write Tranport Simulation files to directory: {}".format(sim.task_root)) |
270 | 268 | sim.write_input() |
271 | 269 |
|
272 | | -sim.msh.export_mesh(# export binary conductivity structure for vtk-visualization |
| 270 | +sim.msh.export_mesh( # export binary conductivity structure for vtk-visualization |
273 | 271 | filepath="{}/conductivity.vtu".format(sim.task_root), |
274 | 272 | file_format="vtk", |
275 | 273 | cell_data_by_id={"transmissivity": BI.kk_mesh}, |
|
284 | 282 | if success: |
285 | 283 | print("Simulation finished successfully.") |
286 | 284 | it_out = -1 |
287 | | - |
| 285 | + |
288 | 286 | out = sim.readvtk() |
289 | 287 | times = out[""]["TIME"] |
290 | 288 | points = out[""]["DATA"][1]["points"] |
|
301 | 299 | ax = fig.add_subplot(1, 1, 1) |
302 | 300 | colormap = plt.get_cmap("bone_r") |
303 | 301 |
|
304 | | -# ixmax = int(100 / domain["dx"]) # limit plot in x-direction to 80m |
305 | | - ixmax = int(200 / domain["dx"]) # limit plot in x-direction to 80m |
| 302 | + ixmax = int(200 / domain["dx"]) |
306 | 303 | ax.contourf( |
307 | 304 | xmesh[:, :ixmax], |
308 | 305 | zmesh[:, :ixmax], |
|
0 commit comments