Commitin

notebooks/ES_2D_Heat_Equation.py

@@ -53,7 +53,6 @@
 
 # Number of grid-cells in x and y direction
 nx = 10
-ny = 10
 
 # time steps
 k_start = 0
@@ -71,7 +70,7 @@
 
 # %%
 def sample_prior_perm(N):
-    mesh = np.meshgrid(np.linspace(0, 1, nx), np.linspace(0, 1, ny))
+    mesh = np.meshgrid(np.linspace(0, 1, nx), np.linspace(0, 1, nx))
     lperms = np.exp(geostat.gaussian_fields(mesh, rng, N, r=0.8))
     return lperms
 
@@ -84,7 +83,7 @@ def sample_prior_perm(N):
 # A list is also a bit easier to interactively visualize.
 alphas = []
 for j in range(N):
-    alphas.append(A[:, j].reshape(nx, ny))
+    alphas.append(A[:, j].reshape(nx, nx))
 
 
 # %%
@@ -110,8 +109,9 @@ def interactive_prior_fields(n):
 
 # %%
 # Set the coefficient of heat transfer for each grid cell.
-# Let's set the truth to be one realization from the prior to make it 
-# easier for the ensemble smoother to find a solution.
+# alpha_t = sample_prior_perm(1).T.reshape(nx, nx)
+# Let's use as true parameter field one relization from the prior
+# to make it easier for the Ensemble Smoother to find the solution.
 alpha_t = alphas[0]
 
 # Resultion on the `x` direction (nothing to worry about really)
@@ -125,7 +125,7 @@ def interactive_prior_fields(n):
 # Define initial condition, i.e., the initial temperature distribution.
 # How you define initial conditions will effect the spread of results,
 # i.e., how similar different realisations are.
-u_init = np.zeros((k_end, nx, ny))
+u_init = np.zeros((k_end, nx, nx))
 u_init[:, 5:7, 5:7] = 100
 
 # How much noise to add to heat equation, also called model noise.
@@ -174,11 +174,11 @@ def interactive_truth(k):
 # %%
 # We'll place the sensors close to heat-sources because the plate cools of quite quickly.
 obs_coordinates = [
-    utils.Coordinate(5, 5),
+    utils.Coordinate(5, 3),
     utils.Coordinate(3, 5),
-    utils.Coordinate(5, 2),
-    utils.Coordinate(8, 8),
     utils.Coordinate(5, 7),
+    utils.Coordinate(7, 5),
+    utils.Coordinate(2, 2),
     utils.Coordinate(7, 2),
 ]
 
@@ -214,12 +214,10 @@ def interactive_truth(k):
 fig, ax = plt.subplots()
 time_to_plot_at = 20
 p = ax.pcolormesh(u_t[time_to_plot_at].T, cmap=plt.cm.jet, vmin=0, vmax=100)
-#ax.invert_yaxis()
+# ax.invert_yaxis()
 ax.set_title(f"True temperature field at t = {time_to_plot_at} with sensor placement")
 utils.colorbar(p)
-ax.plot(
-    [i + 0.5 for i in x], [j + 0.5 for j in y], "s", color="white", markersize=5
-)
+ax.plot([i + 0.5 for i in x], [j + 0.5 for j in y], "s", color="white", markersize=5)
 fig.tight_layout()
 
 # %% [markdown]
@@ -311,6 +309,7 @@ def interactive_realisations(k, n):
 # Ideally, we want our model to be so good that its responses do not deviate much from what our sensors are telling us.
 # When this is the case, we say that we have achieved good `coverage`.
 
+
 # %%
 def plot_responses(
     time_steps,