datagenerator.fluvsim

  1########################################################################
  2# Code to create numpy 3D array with output from channel simulations
  3# - runs fortran program fluvsim.for for a single channel
  4# -- creates fluvsim parameter file (fluvsim.par)
  5# -- reads parameter file
  6# -- runs fluvsim simulation
  7# -- returns 3D array
  8#    --- Codes for facies output:
  9#    --- 0  is floodplain shale
 10#    --- 1  is channel fill (sand)
 11#    --- 2  is shale channel drape
 12#    --- 3  is levee (mid quality sand)
 13#    --- 4  is crevasse (low quality sand)
 14########################################################################
 15
 16import os
 17import numpy as np
 18from datetime import datetime
 19from datagenerator.util import import_matplotlib
 20
 21plt = import_matplotlib()
 22
 23
 24def create_fluvsim_params(
 25    nx=100,
 26    ny=100,
 27    maxthickness=50,
 28    work_folder="/scratch"
 29) -> None:
 30    """
 31    Create Fluvsim parameter file
 32    -----------------------------
 33
 34    Create fluvsim.par parameter file with 
 35    parameters for fluvsim simulation
 36
 37    Parameters
 38    ----------
 39    nx : int, optional
 40        _description_, by default 100
 41    ny : int, optional
 42        _description_, by default 100
 43    maxthickness : int, optional
 44        Maximum thickness, by default 50
 45    work_folder : str, optional
 46        Default work folder, by default "/scratch"
 47    
 48    Returns
 49    -------
 50    None
 51    """
 52
 53    # create fluvsim.par parameter file with randomly chosen values for some parameters
 54    seed = datetime.now().microsecond
 55
 56    # number of channels
 57    # - compute number of channels (needs 10 to work properly)
 58    num_chans = int(np.random.triangular(left=3, mode=5, right=25) + 0.5)
 59    print("...num_chans = ", num_chans)
 60
 61    # channel depth
 62    # - compute number that range from 1+ to 7+, centered at ~3
 63    # - re-scale for maxthickness/50.
 64    vert_scale_factor = maxthickness / 50.0
 65    # random_chan_thickness = np.random.gamma(shape=15.,scale=3.)/13.5       *  vert_scale_factor
 66    # random_chan_thickness = np.random.gamma(shape=5.,scale=7.)/13.5       *  vert_scale_factor
 67    random_chan_thickness = (
 68        np.random.gamma(shape=5.0, scale=5.65) / 13.5 * vert_scale_factor
 69    )
 70    random_chan_thickness = format(random_chan_thickness, "2.1f")
 71    print("...random_chan_thickness = ", random_chan_thickness)
 72
 73    # channel lateral departure (meander width)
 74    # - re-scale
 75    lateral_scale_factor = nx / 100.0
 76    avg_departure = 175.0 * lateral_scale_factor
 77    avg_departure = format(avg_departure, "5.0f")
 78
 79    # channel sinuosity (meander) wavelength
 80    # - re-scale for low/mid/high (L/M/H)
 81    wavelength_factor = nx / 100.0
 82    avg_length = 800.0 * wavelength_factor
 83    avg_length_l = format(avg_length * 0.7, "5.0f")
 84    avg_length_m = format(avg_length, "5.0f")
 85    avg_length_h = format(avg_length * 1.3, "5.0f")
 86
 87    # levee width
 88    levee_onoff = np.random.binomial(1, 0.6)
 89    # levee width
 90    # - re-scale for low/mid/high (L/M/H)
 91    avg_levee = 240.0 * lateral_scale_factor
 92    avg_levee_l = format(avg_levee * 0.67, "5.0f")
 93    avg_levee_m = format(avg_levee, "5.0f")
 94    avg_levee_h = format(avg_levee * 2.0, "5.0f")
 95    print("...levee widths = ", avg_levee_l, avg_levee_m, avg_levee_h)
 96
 97    param_file_text = (
 98        "                   Parameters for FLUVSIM"
 99        + "\n"
100        + "                   **********************"
101        + "\n"
102        + ""
103        + "\n"
104        + "START OF PARAMETERS:"
105        + "\n"
106        + "data/well01.dat                -file with well conditioning data"
107        + "\n"
108        + "1  2  3  4  5                  -  columns for X, Y, Z, well #, facies"
109        + "\n"
110        + "-1.0       1.0e21              -  trimming limits"
111        + "\n"
112        + "1                              -debugging level: 0,1,2,3"
113        + "\n"
114        + "output/fluvsim.dbg             -file for debugging output"
115        + "\n"
116        + "output/fluvsim.geo             -file for geometric specification"
117        + "\n"
118        + "output/fluvsim.out             -file for simulation output"
119        + "\n"
120        + "output/fluvsim.vp              -file for vertical prop curve output"
121        + "\n"
122        + "output/fluvsim.ap              -file for areal prop map output"
123        + "\n"
124        + "output/fluvsim.wd              -file for well data output"
125        + "\n"
126        + "1                              -number of realizations to generate"
127        + "\n"
128        + format(nx, "<3d")
129        + "   0.0    40.0             -nx,xmn,xsiz - geological coordinates"
130        + "\n"
131        + format(nx, "<3d")
132        + "   0.0    40.0             -ny,ymn,ysiz - geological coordinates"
133        + "\n"
134        + format(int(maxthickness), "<3d")
135        + "          50.0              -nz, average thickness in physical units"
136        + "\n"
137        + format(seed, "<6d")
138        + "                         -random number seed"
139        + "\n"
140        + "1   0   0   1                  -1=on,0=off: global, vert, areal, wells"
141        + "\n"
142        + "1.  1.  1.  1.                 -weighting : global, vert, areal, wells"
143        + "\n"
144        + "1     10   0.05                -maximum iter, max no change, min. obj."
145        + "\n"
146        + "0.0   0.10   3  1  8           -annealing schedule: t0,redfac,ka,k,num"
147        + "\n"
148        + "0.1 0.1 0.1 1.0                -Pert prob: 1on+1off, 1on, 1off, fix well"
149        + "\n"
150        + "   1    1    "
151        + format(levee_onoff, "1d")
152        + "                 -Facies(on): channel, levee, crevasse"
153        + "\n"
154        + "0.30 0.15 0.03                 -Proportion: channel, levee, crevasse"
155        + "\n"
156        + "pcurve.dat                     -  vertical proportion curves"
157        + "\n"
158        + "0                              -     0=net-to-gross, 1=all facies"
159        + "\n"
160        + "1  7  8                        -     column numbers"
161        + "\n"
162        + "arealprop.dat                  -  areal proportion map"
163        + "\n"
164        + "1                              -     0=net-to-gross, 1=all facies"
165        + "\n"
166        + "2  3  4                        -     column numbers"
167        + "\n"
168        + format(num_chans, "<3d")
169        + "                            -maximum number of channels"
170        + "\n"
171        + "-45.0    0.0    45.0           -channel:  orientation (degrees)"
172        + "\n"
173        + avg_departure
174        + "  "
175        + avg_departure
176        + "   "
177        + avg_departure
178        + "           -channel:  sinuosity: average departure"
179        + "\n"
180        + avg_length_l
181        + "  "
182        + avg_length_m
183        + "   "
184        + avg_length_h
185        + "           -channel:  sinuosity: length scale"
186        + "\n"
187        + "  "
188        + random_chan_thickness
189        + "    "
190        + random_chan_thickness
191        + "     "
192        + random_chan_thickness
193        + "           -channel:  thickness"
194        + "\n"
195        + "  1.0    1.0     1.0           -channel:  thickness undulation"
196        + "\n"
197        + "250.0  400.0   650.0           -channel:  thickness undul. length scale"
198        + "\n"
199        + " 50.0  150.0   300.0           -channel:  width/thickness ratio"
200        + "\n"
201        + "  1.0    1.0     1.0           -channel:  width: undulation"
202        + "\n"
203        + "250.0  400.0   550.0           -channel:  width: undulation length scale"
204        + "\n"
205        + "500.0 1500.0  3500.0           -levee:    average width"
206        + "\n"
207        + "  0.1    0.15    0.25          -levee:    average height"
208        + "\n"
209        + "  0.05   0.1     0.3           -levee:    depth below top"
210        + "\n"
211        + " 80.0  200.0   500.0           -crevasse: attachment length"
212        + "\n"
213        + "  0.25   0.5     0.75          -crevasse: relative thickness by channel"
214        + "\n"
215        + "200.0  500.0  4500.0           -crevasse: areal size (diameter)"
216        + "\n"
217    )
218
219    # avg_levee_l+'  '+avg_levee_m+'   '+avg_levee_h+'           -levee:    average width'+'\n'+\
220
221    # write param file to disk
222    params_file = os.path.abspath(os.path.join(work_folder, "fluvsim.par"))
223    with open(params_file, "w") as paramfile:
224        paramfile.write(param_file_text)
225    return
226
227
228def read_fluvsim_output(nx, ny, nz, work_folder="/scratch"):
229    # from scipy.ndimage.morphology import grey_closing
230
231    output_file = os.path.abspath(os.path.join(work_folder, "output", "fluvsim.out"))
232    f = open(output_file)
233    a = f.read().split("\n")
234
235    kk = 0
236    facies = np.zeros((nx, ny, nz), "int")
237    for k in range(nz):
238        for j in range(ny):
239            for i in range(nx):
240                facies[i, j, k] = int(a[kk])
241                kk = kk + 1
242
243    return facies
244
245
246def run_fluvsim(nx=100, ny=100, maxthickness=50, work_folder="/scratch", quiet=True):
247    # set up to work on scratch folder
248    current_dir = os.getcwd()
249    code_dir = os.path.join(os.path.abspath(os.path.dirname(__file__)))
250    code_file = os.path.join(code_dir, "fluvsim.f90")
251    workfolder_code_file = os.path.join(work_folder, "fluvsim.f90")
252    os.system("cp " + code_file + " " + workfolder_code_file)
253    try:
254        workfolder_output = os.path.abspath(os.path.join(work_folder, "output"))
255        os.system("mkdir -p " + workfolder_output + "&>/scratch/outputfile")
256    except OSError:
257        pass
258    os.chdir(work_folder)
259    fortran_code_file = os.path.abspath(os.path.join(work_folder, "fluvsim.f90"))
260    compiled_fortran_code_file = "./fluvsim.o"
261    return_code_2 = os.system(
262        "gfortran " + fortran_code_file + " -o " + compiled_fortran_code_file
263    )
264    if return_code_2 != 0:
265        print("...retrieving fluvsim.o from code repository....")
266        os.system(
267            "cp -p "
268            + os.path.join(code_dir, "fluvsim2.o")
269            + " "
270            + compiled_fortran_code_file
271        )
272
273    # create params file
274    create_fluvsim_params(
275        nx=nx, ny=ny, maxthickness=maxthickness, work_folder=work_folder
276    )
277
278    # run fluvsim
279    os.chdir(work_folder)
280    os.system(
281        compiled_fortran_code_file
282        + "&>"
283        + os.path.join(work_folder, "fluvsim_output.txt")
284    )
285
286    # read fluvsim.out ascii file into numpy array
287    facies = read_fluvsim_output(nx, ny, maxthickness, work_folder=work_folder)
288
289    if quiet:  # TODO save plots to output dir?
290        # plot
291        xsection = 10
292        ysection = 10
293        zsection = 30
294
295        # xsection
296        plt.figure(1, figsize=(15, 9))
297        plt.subplot(2, 1, 1)
298        temp = facies[xsection, :, :]
299        plt.imshow(np.rot90(temp))
300        plt.title("xsection = " + str(xsection))
301
302        # ysection
303        plt.subplot(2, 1, 2)
304        temp = facies[:, ysection, :]
305        plt.imshow(np.rot90(temp))
306        plt.title("ysection = " + str(ysection))
307
308        # zsection
309        plt.figure(2, figsize=(10, 10))
310        temp = facies[:, :, zsection]
311        plt.imshow(temp)
312        plt.title("zsection = " + str(zsection))
313
314        plt.show()
315
316    os.chdir(current_dir)
317
318    return facies
319
320
321if __name__ == "__main__":
322
323    # nx=300
324    # ny=300
325    # maxthickness=800
326    # quiet=False
327    model = run_fluvsim(nx=300, ny=300, maxthickness=800, quiet=False)
328    print("model shape = ", model.shape)