commit 3a8080534dbb60ccc47435d15ed6ce9fd909fc59
parent e1e23e617f46e30a99b7ef75ff8010c3813ab1b2
Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Wed, 1 Feb 2017 14:22:31 -0800
improve output figures
Diffstat:
| M | 1d-test.py | | | 51 | ++++++++++++++++++++++++++++++++++++--------------- |
1 file changed, 36 insertions(+), 15 deletions(-)
diff --git a/1d-test.py b/1d-test.py
@@ -17,10 +17,9 @@ import matplotlib.pyplot as plt
import sys
## Model parameters
-#Ns = 10 # Number of nodes [-]
-Ns = 25 # Number of nodes [-]
+Ns = 25 # Number of nodes [-]
Ls = 100e3 # Model length [m]
-t_end = 24.*60.*60.*7. # Total simulation time [s]
+t_end = 24.*60.*60.*1.9 # Total simulation time [s]
tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
tol_P_c = 1e-3 # Tolerance criteria for the normalized max. residual for P_c
max_iter = 1e2*Ns # Maximum number of solver iterations before failure
@@ -222,16 +221,37 @@ def pressure_solver(psi, F, Q, S):
def plot_state(step, time):
# Plot parameters along profile
- plt.plot(s_c/1000., S, '-k', label='$S$')
- plt.plot(s_c/1000., S_max, '--k', label='$S_{max}$')
- plt.plot(s_c/1000., Q, '-b', label='$Q$')
- #plt.plot(s_c/1000., N_c/1000., '--r', label='$N_c$')
- plt.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
- #plt.plot(s, b, ':k', label='$b$')
- #plt.plot(s, b, ':k', label='$b$')
- plt.legend()
- plt.xlabel('Distance from terminus [km]')
- plt.title('Day: {:.2}'.format(time/(60.*60.*24.)))
+ fig = plt.gcf()
+ fig.set_size_inches(3.3, 3.3)
+
+ ax_Pa = plt.subplot(2, 1, 1) # axis with Pascals as y-axis unit
+ ax_Pa.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
+
+ ax_m3s = ax_Pa.twinx() # axis with meters as y-axis unit
+ ax_m3s.plot(s_c/1000., Q, '-b', label='$Q$')
+
+ plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
+ ax_Pa.legend(loc=2)
+ ax_m3s.legend(loc=1)
+ ax_Pa.set_ylabel('[kPa]')
+ ax_m3s.set_ylabel('[m$^3$/s]')
+
+ ax_m = plt.subplot(2, 1, 2, sharex=ax_Pa)
+ ax_m.plot(s_c/1000., S, '-k', label='$S$')
+ ax_m.plot(s_c/1000., S_max, '--k', label='$S_{max}$')
+
+ ax_ms = ax_m.twinx()
+ ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
+ ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
+
+ ax_m.legend(loc=2)
+ ax_ms.legend(loc=1)
+ ax_m.set_xlabel('$s$ [km]')
+ ax_m.set_ylabel('[m]')
+ ax_ms.set_ylabel('[m/s]')
+
+
+ plt.setp(ax_Pa.get_xticklabels(), visible=False)
plt.tight_layout()
if step == -1:
plt.savefig('chan-0.init.pdf')
@@ -241,7 +261,8 @@ def plot_state(step, time):
def find_new_timestep(ds, Q, S):
# Determine the timestep using the Courant-Friedrichs-Lewy condition
- return numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
+ safety = 0.8
+ return safety*numpy.minimum(60.*60.*24., numpy.min(numpy.abs(ds/(Q*S))))
def print_status_to_stdout(time, dt):
sys.stdout.write('\rt = {:.2} s or {:.4} d, dt = {:.2} s'.format(\
@@ -269,7 +290,7 @@ Q = channel_water_flux(S, hydro_pot_grad)
psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
# Prepare figure object for plotting during the simulation
-fig = plt.figure('channel', figsize=(3.3, 2.))
+fig = plt.figure('channel')
plot_state(-1, 0.0)
#import ipdb; ipdb.set_trace()
