Merge pull request #1 from GeoStat-Examples/theis_analysis

Theis analysis
GeoStat-Examples · Apr 27, 2021 · 47c2694 · 47c2694
2 parents 7fd2bda + aae1f13
commit 47c2694
Show file tree

Hide file tree

Showing 93 changed files with 236 additions and 8 deletions.
diff --git a/README.md b/README.md
@@ -60,16 +60,20 @@ The workflow is organized by the following structure:
 
 - `data/`
   - contains the campaign files for both sites in the `welltestpy` format
+  - contains time series for diagnostic plots
 - `src/` - contains the scripts to produce the results
   - `00_wtp_plot.py` - plotting well-constellation and campaign overviews
   - `01_est_run.sh` - bash file running `02_para_estimation.py` in parallel
+  - `01b_est_run.sh` - bash file running `02b_para_estimation.py` in parallel
   - `02_para_estimation.py` - estimate parameters of heterogeneity from the pumping tests
+  - `02b_para_estimation.py` - estimate equivalent parameters of homogeneity the pumping tests
   - `03_postpro_results.py` - plotting the estimation results for both sites
   - `04_postpro_sensitivity.py` - plotting the sensitivity results for both sites
   - `05_est_radial_sens.sh` - bash file running `06_rad_sens.py` in parallel
   - `06_rad_sens.py` - estimate parameter sensitivites depending on the radial distance
     to the pumping well. when run in serial, results will be plotted.
   - `07_comparison_len_scale.py` - generate comparison plot for different length scales
+  - `08_check_unconfined_effect.py` - generate diagnostic plots
 - `results/` - all produced results
 
 

diff --git a/results/01_estimate/HH_results.pdf b/results/01_estimate/HH_results.pdf
diff --git a/results/01_estimate/LW_results.pdf b/results/01_estimate/LW_results.pdf
diff --git a/results/01b_estimate/B2/estimate.txt b/results/01b_estimate/B2/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.041889357980989139e+00
+-2.135803320654749093e+00
diff --git a/results/01b_estimate/B2/fit.pdf b/results/01b_estimate/B2/fit.pdf
diff --git a/results/01b_estimate/B2/parainteract.pdf b/results/01b_estimate/B2/parainteract.pdf
diff --git a/results/01b_estimate/B2/paratrace.pdf b/results/01b_estimate/B2/paratrace.pdf
diff --git a/results/01b_estimate/B2/sens.pdf b/results/01b_estimate/B2/sens.pdf
diff --git a/results/01b_estimate/B2/sens_estimate.txt b/results/01b_estimate/B2/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.245823995020188146e-01
+3.719771898965359425e-01
diff --git a/results/01b_estimate/B2/sens_estimate_S1.txt b/results/01b_estimate/B2/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+5.038909757147155277e-01
+4.920174665307080336e-02
diff --git a/results/01b_estimate/B2/sens_trace.pdf b/results/01b_estimate/B2/sens_trace.pdf
diff --git a/results/01b_estimate/B3/estimate.txt b/results/01b_estimate/B3/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.599175133715375452e+00
+-2.087130985265374772e+00
diff --git a/results/01b_estimate/B3/fit.pdf b/results/01b_estimate/B3/fit.pdf
diff --git a/results/01b_estimate/B3/parainteract.pdf b/results/01b_estimate/B3/parainteract.pdf
diff --git a/results/01b_estimate/B3/paratrace.pdf b/results/01b_estimate/B3/paratrace.pdf
diff --git a/results/01b_estimate/B3/sens.pdf b/results/01b_estimate/B3/sens.pdf
diff --git a/results/01b_estimate/B3/sens_estimate.txt b/results/01b_estimate/B3/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.225084243636654247e-01
+3.793486775373275943e-01
diff --git a/results/01b_estimate/B3/sens_estimate_S1.txt b/results/01b_estimate/B3/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+4.869093878870895398e-01
+4.923811270271956364e-02
diff --git a/results/01b_estimate/B3/sens_trace.pdf b/results/01b_estimate/B3/sens_trace.pdf
diff --git a/results/01b_estimate/B4/estimate.txt b/results/01b_estimate/B4/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.469906947665515773e+00
+-2.558349156725309381e+00
diff --git a/results/01b_estimate/B4/fit.pdf b/results/01b_estimate/B4/fit.pdf
diff --git a/results/01b_estimate/B4/parainteract.pdf b/results/01b_estimate/B4/parainteract.pdf
diff --git a/results/01b_estimate/B4/paratrace.pdf b/results/01b_estimate/B4/paratrace.pdf
diff --git a/results/01b_estimate/B4/sens.pdf b/results/01b_estimate/B4/sens.pdf
diff --git a/results/01b_estimate/B4/sens_estimate.txt b/results/01b_estimate/B4/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.248095705848720094e-01
+3.488981377084626745e-01
diff --git a/results/01b_estimate/B4/sens_estimate_S1.txt b/results/01b_estimate/B4/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+5.003059063840357146e-01
+4.880588355868517580e-02
diff --git a/results/01b_estimate/B4/sens_trace.pdf b/results/01b_estimate/B4/sens_trace.pdf
diff --git a/results/01b_estimate/B5/estimate.txt b/results/01b_estimate/B5/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-3.687795829848876217e+00
+-3.341085358753064227e+00
diff --git a/results/01b_estimate/B5/fit.pdf b/results/01b_estimate/B5/fit.pdf
diff --git a/results/01b_estimate/B5/parainteract.pdf b/results/01b_estimate/B5/parainteract.pdf
diff --git a/results/01b_estimate/B5/paratrace.pdf b/results/01b_estimate/B5/paratrace.pdf
diff --git a/results/01b_estimate/B5/sens.pdf b/results/01b_estimate/B5/sens.pdf
diff --git a/results/01b_estimate/B5/sens_estimate.txt b/results/01b_estimate/B5/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.432098411371168112e-01
+2.745695629845287256e-01
diff --git a/results/01b_estimate/B5/sens_estimate_S1.txt b/results/01b_estimate/B5/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+6.384341518194548426e-01
+4.890268754194360828e-02
diff --git a/results/01b_estimate/B5/sens_trace.pdf b/results/01b_estimate/B5/sens_trace.pdf
diff --git a/results/01b_estimate/HH_compare_S.pdf b/results/01b_estimate/HH_compare_S.pdf
diff --git a/results/01b_estimate/HH_compare_T.pdf b/results/01b_estimate/HH_compare_T.pdf
diff --git a/results/01b_estimate/HH_results.pdf b/results/01b_estimate/HH_results.pdf
diff --git a/results/01b_estimate/LW_compare_S.pdf b/results/01b_estimate/LW_compare_S.pdf
diff --git a/results/01b_estimate/LW_compare_T.pdf b/results/01b_estimate/LW_compare_T.pdf
diff --git a/results/01b_estimate/LW_results.pdf b/results/01b_estimate/LW_results.pdf
diff --git a/results/01b_estimate/all_hh/estimate.txt b/results/01b_estimate/all_hh/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.616495101489531905e+00
+-1.814607315343006677e+00
diff --git a/results/01b_estimate/all_hh/fit.pdf b/results/01b_estimate/all_hh/fit.pdf
diff --git a/results/01b_estimate/all_hh/parainteract.pdf b/results/01b_estimate/all_hh/parainteract.pdf
diff --git a/results/01b_estimate/all_hh/paratrace.pdf b/results/01b_estimate/all_hh/paratrace.pdf
diff --git a/results/01b_estimate/all_hh/sens.pdf b/results/01b_estimate/all_hh/sens.pdf
diff --git a/results/01b_estimate/all_hh/sens_estimate.txt b/results/01b_estimate/all_hh/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.269326284369937552e-01
+3.549613290790790687e-01
diff --git a/results/01b_estimate/all_hh/sens_estimate_S1.txt b/results/01b_estimate/all_hh/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+5.410627948848325541e-01
+5.298159575586086634e-02
diff --git a/results/01b_estimate/all_hh/sens_trace.pdf b/results/01b_estimate/all_hh/sens_trace.pdf
diff --git a/results/01b_estimate/all_lw/estimate.txt b/results/01b_estimate/all_lw/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.254505512021658831e+00
+-2.398455627950027313e+00
diff --git a/results/01b_estimate/all_lw/fit.pdf b/results/01b_estimate/all_lw/fit.pdf
diff --git a/results/01b_estimate/all_lw/parainteract.pdf b/results/01b_estimate/all_lw/parainteract.pdf
diff --git a/results/01b_estimate/all_lw/paratrace.pdf b/results/01b_estimate/all_lw/paratrace.pdf
diff --git a/results/01b_estimate/all_lw/sens.pdf b/results/01b_estimate/all_lw/sens.pdf
diff --git a/results/01b_estimate/all_lw/sens_estimate.txt b/results/01b_estimate/all_lw/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.237182492940093637e-01
+3.122240325370111425e-01
diff --git a/results/01b_estimate/all_lw/sens_estimate_S1.txt b/results/01b_estimate/all_lw/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+4.959178374954668067e-01
+4.916851923696172189e-02
diff --git a/results/01b_estimate/all_lw/sens_trace.pdf b/results/01b_estimate/all_lw/sens_trace.pdf
diff --git a/results/01b_estimate/p05/estimate.txt b/results/01b_estimate/p05/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.108162874320472113e+00
+-2.363775493525403171e+00
diff --git a/results/01b_estimate/p05/fit.pdf b/results/01b_estimate/p05/fit.pdf
diff --git a/results/01b_estimate/p05/parainteract.pdf b/results/01b_estimate/p05/parainteract.pdf
diff --git a/results/01b_estimate/p05/paratrace.pdf b/results/01b_estimate/p05/paratrace.pdf
diff --git a/results/01b_estimate/p05/sens.pdf b/results/01b_estimate/p05/sens.pdf
diff --git a/results/01b_estimate/p05/sens_estimate.txt b/results/01b_estimate/p05/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.160116039704321311e-01
+3.383748292193309304e-01
diff --git a/results/01b_estimate/p05/sens_estimate_S1.txt b/results/01b_estimate/p05/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+4.587503507086950205e-01
+5.934418044628082811e-02
diff --git a/results/01b_estimate/p05/sens_trace.pdf b/results/01b_estimate/p05/sens_trace.pdf
diff --git a/results/01b_estimate/p40/estimate.txt b/results/01b_estimate/p40/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-5.022233358482504073e+00
+-1.668808864276526194e+00
diff --git a/results/01b_estimate/p40/fit.pdf b/results/01b_estimate/p40/fit.pdf
diff --git a/results/01b_estimate/p40/parainteract.pdf b/results/01b_estimate/p40/parainteract.pdf
diff --git a/results/01b_estimate/p40/paratrace.pdf b/results/01b_estimate/p40/paratrace.pdf
diff --git a/results/01b_estimate/p40/sens.pdf b/results/01b_estimate/p40/sens.pdf
diff --git a/results/01b_estimate/p40/sens_estimate.txt b/results/01b_estimate/p40/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.222609064596241080e-01
+3.693173255790749288e-01
diff --git a/results/01b_estimate/p40/sens_estimate_S1.txt b/results/01b_estimate/p40/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+5.012227180660423498e-01
+5.500582275881459038e-02
diff --git a/results/01b_estimate/p40/sens_trace.pdf b/results/01b_estimate/p40/sens_trace.pdf
diff --git a/results/01b_estimate/p42/estimate.txt b/results/01b_estimate/p42/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.184148473481728026e+00
+-2.259962264311639046e+00
diff --git a/results/01b_estimate/p42/fit.pdf b/results/01b_estimate/p42/fit.pdf
diff --git a/results/01b_estimate/p42/parainteract.pdf b/results/01b_estimate/p42/parainteract.pdf
diff --git a/results/01b_estimate/p42/paratrace.pdf b/results/01b_estimate/p42/paratrace.pdf
diff --git a/results/01b_estimate/p42/sens.pdf b/results/01b_estimate/p42/sens.pdf
diff --git a/results/01b_estimate/p42/sens_estimate.txt b/results/01b_estimate/p42/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.468237177786873460e-01
+3.401689381372080545e-01
diff --git a/results/01b_estimate/p42/sens_estimate_S1.txt b/results/01b_estimate/p42/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+6.919797689549598285e-01
+5.663505775749034382e-02
diff --git a/results/01b_estimate/p42/sens_trace.pdf b/results/01b_estimate/p42/sens_trace.pdf
diff --git a/results/01b_estimate/p44/estimate.txt b/results/01b_estimate/p44/estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+-4.846352643774878288e+00
+-1.702396323325426231e+00
diff --git a/results/01b_estimate/p44/fit.pdf b/results/01b_estimate/p44/fit.pdf
diff --git a/results/01b_estimate/p44/parainteract.pdf b/results/01b_estimate/p44/parainteract.pdf
diff --git a/results/01b_estimate/p44/paratrace.pdf b/results/01b_estimate/p44/paratrace.pdf
diff --git a/results/01b_estimate/p44/sens.pdf b/results/01b_estimate/p44/sens.pdf
diff --git a/results/01b_estimate/p44/sens_estimate.txt b/results/01b_estimate/p44/sens_estimate.txt
@@ -0,0 +1,3 @@
+# mu lnS
+9.160779984425324596e-01
+4.159615368057123508e-01
diff --git a/results/01b_estimate/p44/sens_estimate_S1.txt b/results/01b_estimate/p44/sens_estimate_S1.txt
@@ -0,0 +1,3 @@
+# mu lnS
+4.582454345407614427e-01
+5.645156912666832866e-02
diff --git a/results/01b_estimate/p44/sens_trace.pdf b/results/01b_estimate/p44/sens_trace.pdf
diff --git a/src/01b_est_run.sh b/src/01b_est_run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+mpiexec -n 4 python3 02b_para_estimation_theis.py
diff --git a/src/02b_para_estimation_theis.py b/src/02b_para_estimation_theis.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""Estimate HH+LW: single test and combined test estimation."""
+
+import os
+from mpi4py import MPI
+import welltestpy as wtp
+import matplotlib.pyplot as plt
+
+# increase fontsize of plots, prevent type 3 fonts in pdf output
+plt.rcParams.update({"font.size": 16, "pdf.fonttype": 42, "ps.fonttype": 42})
+
+# rank is the actual core-number, size is total number of cores
+rank = MPI.COMM_WORLD.Get_rank()
+size = MPI.COMM_WORLD.Get_size()
+
+# saving estimation results in results
+here = os.path.abspath(os.path.dirname(__file__))
+results = os.path.normpath(os.path.join(here, "..", "results", "01b_estimate"))
+data = os.path.normpath(os.path.join(here, "..", "data"))
+os.makedirs(results, exist_ok=True)
+
+# setups for single test and combined estimation
+setups = [
+    [os.path.join(data, "lauswiesen.cmp"), "all_lw"],
+    [os.path.join(data, "lauswiesen.cmp"), "B2"],
+    [os.path.join(data, "lauswiesen.cmp"), "B3"],
+    [os.path.join(data, "lauswiesen.cmp"), "B4"],
+    [os.path.join(data, "lauswiesen.cmp"), "B5"],
+    [os.path.join(data, "horkheim.cmp"), "all_hh"],
+    [os.path.join(data, "horkheim.cmp"), "p05"],
+    [os.path.join(data, "horkheim.cmp"), "p40"],
+    [os.path.join(data, "horkheim.cmp"), "p42"],
+    [os.path.join(data, "horkheim.cmp"), "p44"],
+]
+
+for i, setup in enumerate(setups):
+    if i % size != rank:
+        continue
+    campaign_file, well = setup
+    print(well + " on core {}/{} started".format(rank + 1, size))
+    # load the lauswiesen pumpingtest
+    cmp = wtp.load_campaign(campaign_file)
+    # get the included tests
+    is_all = well.startswith("all")
+    testinclude = {well: cmp.tests[well].wells} if not is_all else None
+    # set up the estimation
+    estimation = wtp.estimate.Theis(
+        "est", cmp, testinclude=testinclude, generate=True
+    )
+    # run the estimation
+    estimation.run(
+        rep=2500,
+        folder=os.path.join(results, well),
+        dbname="db",
+        traceplotname="paratrace.pdf",
+        fittingplotname="fit.pdf",
+        interactplotname="parainteract.pdf",
+        estname="estimate.txt",
+        # run=False,
+    )
+    # estimate the sensitivites
+    estimation.sensitivity(
+        folder=os.path.join(results, well),
+        dbname="sens_db",
+        plotname="sens.pdf",
+        traceplotname="sens_trace.pdf",
+        sensname="sens_estimate.txt",
+    )
diff --git a/src/03_postpro_results.py b/src/03_postpro_results.py
@@ -29,24 +29,30 @@ def plot(site, root=None):
     ptests = wells[site]
     estim = os.path.join(root, "all_" + site, "estimate.txt")
     testdict = {}
-    print(estim)
+    # print(estim)
     # load data from file
     testdict["all"] = np.loadtxt(estim)
+    val_cnt = 4
+    if len(testdict["all"]) == 2:
+        val_cnt = 2
+        testdict["all"] = np.insert(testdict["all"], 1, [0, 0])
 
     for p in ptests:
         estim = os.path.join(root, p, "estimate.txt")
-        print(estim)
+        # print(estim)
         # load data from file
         testdict[p] = np.loadtxt(estim)
+        if len(testdict[p]) == 2:
+            testdict[p] = np.insert(testdict[p], 1, [0, 0])
         # if var is 0 --> len_scale is 0
         if testdict[p][1] < 0.02:
-            print(testdict[p][1])
+            # print(testdict[p][1])
             testdict[p][2] = 0
 
     lin = np.ones(2)
     keys = ["all", "mean"] + ptests  # +["mean"]#, "geo-mean"]
     varnames = [
-        r"$T^G_{all}$",
+        r"$T^G_{all}$" if val_cnt == 4 else r"$T_{all}$",
         r"$\sigma^2_{all}$",
         r"$\ell_{all}$",
         r"$S_{all}$",
@@ -57,11 +63,21 @@ def plot(site, root=None):
         r" $m$",
         r"",
     ]
-
+    ret1 = []  # T[_H]
+    ret2 = []  # S
     for res in testdict:
         testdict[res][0] = np.exp(testdict[res][0])
         testdict[res][3] = np.exp(testdict[res][3])
 
+    for res in keys:
+        if res == "mean":
+            continue
+        if val_cnt == 2:
+            ret1.append(testdict[res][0])
+        else:
+            ret1.append(testdict[res][0] * np.exp(-testdict[res][1] / 2))
+        ret2.append(testdict[res][-1])
+
     # rescale to common y-range
     max_y = 2.0
     for res in ptests:
@@ -91,7 +107,9 @@ def plot(site, root=None):
 
     fig = plt.figure(dpi=75, figsize=[7.5, 4])
     for j, var in enumerate(varnames):
-        ax = fig.add_subplot(1, len(varnames), j + 1)
+        if val_cnt < 4 and j in [1, 2]:
+            continue
+        ax = fig.add_subplot(1, val_cnt, (min(j, 1) if val_cnt < 4 else j) + 1)
         for i, res in enumerate(keys):
             if i < 2:
                 ax.plot(
@@ -141,7 +159,51 @@ def plot(site, root=None):
         os.path.join(root, site.upper() + "_results" + file_ext),
         bbox_inches="tight",
     )
+    return ret1, ret2, ["all"] + ptests
+
+
+def compare_vars(
+    v1, v2, names, site, root, lab1="T_H", lab2="T", unit="$m^2/s$"
+):
+    """Compare variables."""
+    fig, ax = plt.subplots(dpi=75, figsize=[7.5, 4])
+    x = np.arange(len(names))
+    w = 0.4
+    ax.bar(x - w / 2, v1, w, alpha=0.8, label=f"${lab1}$ from ext. Theis")
+    ax.bar(x + w / 2, v2, w, alpha=0.8, label=f"${lab2}$ from Theis")
+    ax.set_xticks(x)
+    ax.set_xticklabels(names)
+    ax.set_ylabel(f"{lab2} in [{unit}]", fontsize=16)
+    ax.grid(axis="x", which="both")
+    legend = ax.get_legend_handles_labels()
+    fig.legend(
+        *legend,
+        loc="lower center",
+        ncol=2,
+        bbox_to_anchor=(0.5, 0),
+        handlelength=1,
+        columnspacing=1,
+        handletextpad=0.5,
+    )
+    fig.tight_layout()
+    fig.subplots_adjust(wspace=0.1, bottom=0.3)
+    fig.show()
+    fig.savefig(
+        os.path.join(root, site.upper() + "_compare_" + lab2 + file_ext),
+        bbox_inches="tight",
+    )
+    return fig
+
+
+root = os.path.normpath(os.path.join(here, "..", "results", "01_estimate"))
+TH_HH, S1_HH, names_HH = plot("HH", root)
+TH_LW, S1_LW, names_LW = plot("LW", root)
+root = os.path.normpath(os.path.join(here, "..", "results", "01b_estimate"))
+TG_HH, S2_HH, __ = plot("HH", root)
+TG_LW, S2_LW, __ = plot("LW", root)
 
+compare_vars(TH_HH, TG_HH, names_HH, "HH", root)
+compare_vars(S1_HH, S2_HH, names_HH, "HH", root, lab1="S", lab2="S", unit="-")
 
-plot("HH")
-plot("LW")
+compare_vars(TH_LW, TG_LW, names_LW, "LW", root)
+compare_vars(S1_LW, S2_LW, names_LW, "LW", root, lab1="S", lab2="S", unit="-")
diff --git a/src/08_check_unconfined_effect.py b/src/08_check_unconfined_effect.py
@@ -3,6 +3,7 @@
 import os
 import numpy as np
 from matplotlib import pyplot as plt
+
 # plotting style
 plt.style.use("ggplot")
 # increase fontsize of plots, prevent type 3 fonts in pdf output
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,3 @@
		#!/bin/bash

		mpiexec -n 4 python3 02b_para_estimation_theis.py