#!/usr/bin/env python3

import grass.script as gs

# set region
gs.run_command("g.region", rows=100, cols=100, res3=1, t=50, b=0, n=1000,s=0, w=0,\
 e=1000, overwrite=True)

# generate simple raster maps
gs.run_command("r.mapcalc", expression="phead = if(col() == 1, 39, 39-0.001*x())", overwrite=True)
gs.run_command("r.mapcalc", expression="status = if(col() == 1 || col() == 100,\
 2, 1)", overwrite=True)
gs.run_command("r.mapcalc", expression="null = 0.0",overwrite=True)
gs.run_command("r.mapcalc", expression="top = 25",overwrite=True)

# generate artificial heterogeneity
#gs.run_command("r.random.surface", output="heterogenMap", distance=30.0, 
#flat=10.0, seed=500, high=100, overwrite=True)

# apply heterogeneity on maps
# use the statistics tool to calculate the mean of heterogenMap (46.61), 
# subtract the mean and divide by 100 to obtain a distribution of cell values
# between -0.5 and 0.5. Then apply a factor (6) to achieve the desired 
# standard deviation of 0.000075 m/s. Use the statistics tool to verify
gs.run_command("r.mapcalc", expression="hydcond = 0.0015",overwrite=True)

gs.run_command("r.mapcalc", expression="poros = 0.0001",overwrite=True)

gs.run_command("r.mapcalc", expression="q = if (row()==50 && col()==50, -0.01, 0)",overwrite=True)

# Compute a steady state groundwater flow without source
gs.run_command("r.gwflow", top="top", bottom="null",phead="phead",
    status="status", hc_x="hydcond", hc_y="hydcond", s="poros",
    output="gwresult_0", dt=8640000000, type="confined",
    vx="velx", vy="vely", overwrite=True)

gs.run_command("r.gwflow", top="top", bottom="null",phead="gwresult_0",
    status="status", hc_x="hydcond", hc_y="hydcond", s="poros", q="q",
    output="gwresult_q", dt=864000000000, type="confined",
    vx="velx", vy="vely", 
    error="0.00000000000001",
    maxit="100000000000", overwrite=True)

gs.run_command("r.gwflow", top="top", bottom="null",phead="gwresult_q",
    status="status", hc_x="hydcond", hc_y="hydcond", s="poros", q="q",
    output="gwresult_q", dt=8640000000, type="confined",
    vx="velx", vy="vely", 
    error="0.00000000001",
    maxit="100000000000", overwrite=True)