Revised plotters.

dapper/mods/Stommel/__init__.py

@@ -17,10 +17,7 @@
 #Directory to store figures. 
 fig_dir = "<dir for storing figures>"
 DIR = "<topdir containing files downloaded from server or containing dirs with those files>"
-
-fig_dir = '/home/ivo/Figures/stommel/hadley0830'
-DIR = '/home/ivo/dpr_data/stommel'
-FILE_NAME = 'boxed_hadley_inverted0830.pkl'
+FILE_NAME = 'boxed_hadley_inverted0906.pkl'
 
 
 hadley_file = os.path.join(DIR, FILE_NAME)
@@ -32,7 +29,7 @@
     hadley = pkl.load(stream)
     #dy = .5*hadley['geo_pole']['dy'] + .5*hadley['geo_eq']['dy']
     #hadley['geo_pole']['dy'], hadley['geo_eq']['dy'] = dy, dy
-    #hadley['geo_pole']['dz'], hadley['geo_eq']['dz'] = 4.2e3, 4.2e3
+    #hadley['geo_pole']['dz'], hadley['geo_eq']['dz'] = 1.5e3, 1.5e3
 
     ref = hadley['yy'][0] #np.mean(hadley['yy'], axis=0)
 
@@ -449,12 +446,10 @@ def ens_member(self, member):
 
     def __post_init__(self):
         """Part of object initialization to be carried out after __init__ provided by dataclass."""
-        #self.init_state.gamma = self.default_gamma(self.init_state)
-        #self.state.gamma = self.default_gamma(self.state)
 
-        self.init_state.temp = np.array([[2.57,9.09]])
-        self.init_state.salt = np.array([[34.738,35.455]])
-        self.init_state = self.default_parameters1(self.init_state, Q_overturning)
+        self.init_state.temp = np.array([hadley['yy'][0,0:2]])
+        self.init_state.salt = np.array([hadley['yy'][0,2:4]])
+        self.init_state = self.default_parameters3(self.init_state, Q_overturning)
         self.state = copy(self.init_state)
         self.fluxes = self.default_fluxes(self.fluxes) 
 
@@ -511,9 +506,7 @@ def equi(params):
         trans_eq = self.trans_eq(state)
         temp_eq = self.temp_eq(state, trans_eq)
         salt_eq = self.salt_eq(state, trans_eq)
-        print('EQ ',trans_eq*1e-6, temp_eq, salt_eq)
-        #state.salt[0] = np.mean(state.salt[0]) + np.array([-.5,.5]) * salt_eq 
-        #state.temp[0] = np.mean(state.temp[0]) + np.array([-.5,.5]) * temp_eq 
+        print('Initial Q, diff T, diff S',trans_eq*1e-6, temp_eq, salt_eq)
 
         return state
 
@@ -553,13 +546,12 @@ def equi(params):
         trans_eq = self.trans_eq(state)
         temp_eq = self.temp_eq(state, trans_eq)
         salt_eq = self.salt_eq(state, trans_eq)
-        print('EQ ',trans_eq*1e-6, temp_eq, salt_eq)
+        print('Initial Q, diff T, diff S',trans_eq*1e-6, temp_eq, salt_eq)
 
         return state
 
     def default_parameters1(self, state, Q=Q_overturning):
-        """ Find parameter values such that current state is an equilibrium
-        with transport Q. """
+        """ Set parameters.  """
 
         state.gamma = self.default_gamma(state, Q)
 
@@ -571,7 +563,7 @@ def default_parameters1(self, state, Q=Q_overturning):
         trans_eq = self.trans_eq(state)
         temp_eq = self.temp_eq(state, trans_eq)
         salt_eq = self.salt_eq(state, trans_eq)
-        print('EQ ',trans_eq*1e-6, temp_eq, salt_eq)
+        print('Initial Q, diff T, diff S',trans_eq*1e-6, temp_eq, salt_eq)
 
         return state
 

dapper/mods/Stommel/clima_forcing_da.py

@@ -22,7 +22,7 @@ def exp_clima_forcing_da(N=100, seed=1000):
     Tda = 20 * stommel.year #time period over which DA takes place. 
     kko = [] #np.arange(1, len(hadley['yy'])+1)
     tseq = modelling.Chronology(stommel.year/12, kko=kko, 
-                                T=10*stommel.year, BurnIn=0)  # 1 observation/month
+                                T=100*stommel.year, BurnIn=0)  # 1 observation/month
     T_warming=100*stommel.year
     # Timestepping. Timesteps of 1 day, running for 200 year.
     #Create model

dapper/mods/Stommel/experiment.py

@@ -0,0 +1,201 @@
+# -*- coding: utf-8 -*-
+
+""" 
+Run model over a range of climate parameters. 
+"""
+
+import numpy as np
+import dapper.mods as modelling
+import dapper.mods.Stommel as stommel
+from dapper.da_methods.ensemble import EnKF
+from dapper.mods.Stommel import hadley
+from scipy.interpolate import interp1d
+import dill, os
+import multiprocessing
+import concurrent.futures
+
+def climate_temp_forcing(N, T_da, T_warming=0.0, A_polar=0.0):
+    """ Create surface forcing for temperatur flux due to global warming. """
+    #Heat air fluxes 
+    functions = stommel.hadley_air_temp(N)
+    #Add linear warming with 6C/T_warming over the pole and 3C/T_warming over the equatior. 
+    trend = interp1d(np.array([0.,T_da,T_warming+T_da]), 
+                     np.array([[0.,0.,A_polar],[0.,0.,0.5*A_polar]]), 
+                     fill_value='extrapolate', axis=1)
+    trended = [stommel.add_functions(func, trend) for func in functions]
+    #For time<Tda all ensemble member n uses noised[0] after that 
+    #functions = [stommel.merge_functions(T_da, func, trend) 
+    #             for func,trend in zip(functions,trended)]
+    return trended
+
+def climate_ep_forcing(N, model, T_da, T_melt):
+    """ Create evaporation-precipitation forcing due to Greenland ice sheet melt."""
+    #Melt flux 
+    melt_rate = -stommel.V_ice * np.array([1.0/(model.dx[0,0]*model.dy[0,0]), 0.0]) / T_melt #ms-1
+    #Default evaporation-percipitation flux (=0)
+    functions = stommel.default_air_ep(N)
+    #Add effect Greenland melt with annual rate melt_rate
+    def add_melting(func):
+        def func_with_melt(t):
+            if t<T_da:
+                return func(t)
+            elif t<T_da+T_melt:
+                return func(t)+melt_rate
+            else:
+                return func(t)
+        return func_with_melt
+
+    functions = [add_melting(func) for func in functions]
+    return functions
+
+
+def experiment(N=100, seed=1100, do_da=True, 
+               A_polar=0.0, T_warming = np.inf, T_melt=np.inf):
+    #Set seed
+    np.random.seed(seed)
+    # Timestepping. Timesteps of 1 day, running for 200 year.
+    dt = stommel.year/12
+    kko = np.arange(1, len(hadley['yy'][1:])+1) if do_da else np.array([])    
+    tseq = modelling.Chronology(dt, kko=kko, T=100*stommel.year)  # 1 observation/month
+    T_da = np.size(hadley['yy'],0) * dt
+
+    #Activate surface heat flux. functions[n] contains atm. temperature for ensemble member n. 
+    temp_forcings = climate_temp_forcing(N, T_da, T_warming, A_polar)
+
+    #Salinity air fluxes 
+    salt_forcings = stommel.hadley_air_salt(N)
+
+    #Activate EP flux. 
+    ref_model = stommel.StommelModel(fluxes=[stommel.TempAirFlux(temp_forcings),
+                                         stommel.SaltAirFlux(salt_forcings)])
+    ep_forcings = climate_ep_forcing(N, ref_model, T_da, T_melt)
+
+    #Create model
+    model = stommel.StommelModel(fluxes=[stommel.TempAirFlux(temp_forcings),
+                                         stommel.SaltAirFlux(salt_forcings),
+                                         stommel.EPFlux(ep_forcings)])
+
+
+    # Initial conditions
+    x0 = model.x0
+
+    # Variance Ocean temp[2], ocean salinity[2], temp diffusion parameter,
+    # salt diffusion parameter, transport parameter
+    B = stommel.State().zero()
+    B.temp += hadley['R'][:2]  # C2
+    B.salt += hadley['R'][2:]  # ppt2
+    B.temp_diff += np.log(1.3)**2  # (0.0*model.init_state.temp_diff)**2
+    B.salt_diff += np.log(1.3)**2  # (0.0*model.init_state.salt_diff)**2
+    B.gamma += np.log(1.3)**2  # (0.0*model.init_state.gamma)**
+
+    # Transform modus value in x0 to mean value.
+    x0[4] += B.temp_diff
+    x0[5] += B.salt_diff
+    x0[6] += B.gamma
+
+    #Create sampler for initial conditions. 
+    X0 = modelling.GaussRV(C=B.to_vector(), mu=x0)
+    # Dynamisch model. All model error is assumed to be in forcing.
+    Dyn = {'M': model.M,
+           'model': model.step,
+           'noise': 0
+           }
+    # Observation. If factor=1 no correction to obs. error std. dev. is applied. 
+    Obs = model.obs_hadley(factor=np.array([1.,1.,1.,1.]))
+    # Create model.
+    HMM = modelling.HiddenMarkovModel(Dyn, Obs, tseq, X0)
+    # Create DA with inflation
+    xp = EnKF('Sqrt', N, infl=[1.0]*4+[1.01]*3)
+    return xp, HMM, model
+
+def flip_probability(HMM,E):
+    """ Calculate probability of flip at some time. """
+    times = HMM.tseq.times
+    N_flips = 0
+    end_da = np.max(HMM.tseq.kko)
+    for e in E.transpose((1,0,2)):
+        states = stommel.array2states(e[end_da:], times[end_da:])
+        SA = np.array([s.regime=='SA' for s in states], dtype=bool)
+        N_flips += np.any(SA)
+    return float(N_flips) / float(np.size(E,1))
+
+#%% Run the experiment
+
+def run_list(filename, observations=[None], T_melts = [np.inf], 
+             A_polars=[0.0], labels=['nCC-DA']):
+    """
+    Run all experiments in sequence. 
+    """   
+    filepath = os.path.join(stommel.fig_dir, filename+'.pkl')
+    if os.path.exists(filepath):
+        with open(filepath,'rb') as stream:
+            data = dill.load(stream)
+    else:
+        data = []
+
+    for T_melt, A_polar, obs, label in zip(T_melts, A_polars, observations, labels):
+        labels = [data1['label'] for data1 in data]
+        if label in labels:
+            continue
+
+        print('Running ',label)
+        #Create experiment
+        do_da = len(obs)>0
+        xp, HMM, model = experiment(T_melt=T_melt, A_polar=A_polar, 
+                                    T_warming=100.*stommel.year, do_da=do_da)
+        #Run
+        xx, yy = HMM.simulate()
+        if do_da and str(obs[0])!='twin':
+            assert len(yy)==len(obs)
+            yy = obs
+        Efor, Eana = xp.assimilate(HMM, xx, yy)
+
+        data1 = {'model':model,'xp':xp, 'HMM':HMM,
+                 'yy': yy, 'Efor':Efor, 'Eana':Eana,
+                 'T_melt':T_melt, 'A_polar':A_polar, 'label':label}
+
+        if len(obs)==0 or str(obs[0])=='twin':
+            data1 = {**data1, 'xx':xx}
+
+        data.append(data1)
+
+    with open(filepath,'wb') as stream:
+        dill.dump(data, stream)
+
+    return data 
+
+def experiment1(param_in):
+    from time import sleep
+    T_melt, A_polar = param_in
+    #Create experiment
+    xp, HMM, model = experiment(T_melt=T_melt, A_polar=A_polar,
+                                T_warming=100.*stommel.year)
+    # #Run
+    xx, yy = HMM.simulate()
+    yy = hadley['yy'][HMM.tseq.kko]
+    Efor, Eana = xp.assimilate(HMM, xx, yy)
+    #Calculate flip probility
+    return flip_probability(HMM, Efor)
+
+
+def run_hypercube(filename, T_melts = [np.inf], 
+                  A_polars=[0.0]):
+    """ 
+    Run all combinations of settings. 
+    """ 
+
+    filepath = os.path.join(stommel.fig_dir, filename+'.pkl')
+    climates = [(T_melt, A_polar) for T_melt in T_melts for A_polar in A_polars]
+
+    pool = multiprocessing.Pool(processes = 8)
+    probs = pool.map(experiment1, climates)
+    print('PROBS ',probs)
+
+    data = {'climas':climates, 'probabilities':probs}
+    with open(filepath,'wb') as stream:
+        dill.dump(data, stream)
+
+    return data
+
+
+

dapper/mods/Stommel/plot_boxes.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Sep 13 11:25:08 2024
+
+Plot boxes ET4
+
+@author: ivo
+"""
+import os, dill
+import dapper.mods.Stommel.plotters as plotters
+import dapper.mods.Stommel as stommel
+import matplotlib.pyplot as plt
+
+cluster_file = os.path.join(stommel.DIR, 'clusters_boxed_hadley_inverted0906.pkl')
+with open(cluster_file,'rb') as stream:
+    data = dill.load(stream)
+
+#%% 
+
+plt.close('all')
+boxes_plot = plotters.BoxesPlot(data['indices'], data['output'])
+boxes_plot.plot('boxes.png')    

dapper/mods/Stommel/plot_scenarios.py

@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Plot single case scenarios. 
+"""
+
+import dapper.mods.Stommel.plotters as plotters
+import dapper.mods.Stommel.experiment as exp
+import dapper.mods.Stommel as stommel
+import numpy as np
+
+#%% Run twin experiments
+
+data = exp.run_list('twin3', 
+                observations=[[],['twin'],[],['twin'],[]],
+                T_melts=np.array([np.inf,np.inf,1e4,1e4,10e2])*stommel.year,
+                A_polars = [0.0,0.0,6.0,6.0,6.0],
+                labels=['nCC-nDA','nCC-yDA','yCC-nDA','yCC-yDA','yCC-nDA 1000y']
+                )
+
+P = plotters.DiffPhase(data[:-1])
+P.plot('fig04_diff_phase_10000.png')
+P = plotters.ParameterPlot(**data[1])
+P.plot('fig05_parameters_yDA.png')
+P = plotters.RmseStdPlot(data[:-1])
+P.plot('fig06_rmse_std.png')
+P = plotters.DiffPhase(data[-1:])
+P.plot_flip('fig07_diff_phase_1000.png')
+
+#%% Run experiments with real observations.
+
+data = exp.run_list('scenarios3',
+                    observations= [stommel.hadley['yy'][1:]]*2,
+                    T_melts=[np.inf, 1e4*stommel.year],
+                    A_polars = [0.0, 6.0], 
+                    labels=['nCC-EN4','yCC-EN4'])
+
+#Plot output
+P = plotters.ParameterPlot(**data[0])
+P.plot('fig08_parameters_yCC-EN4.png')
+P = plotters.EtaPlot(**data[-1])
+P.plot('fig09_etas_yCC-EN4.png')
+P = plotters.DiffPhase(data[:-1])
+P.plot('fig10_diff_phase_nCC_EN4.png')
+P = plotters.DiffPhase(data[-1:])
+P.plot('fig11_diff_phase_yCC_EN4.png')
+P = plotters.TransportPlot(**data[-1])
+P.plot('fig12_transport_yCC-EN4.png')
+
+

dapper/mods/Stommel/plot_tipping.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Plot probability tipping.
+"""
+
+import dapper.mods.Stommel.plotters as plotters
+from dapper.mods.Stommel.experiment import run_hypercube
+import dapper.mods.Stommel as stommel
+import numpy as np
+import dill,os
+import matplotlib.pyplot as plt
+import dapper.mods.Stommel.plotters as plotter
+
+#%% Run experiment with different levels of forcing. 
+
+#data = run_hypercube('tipping3', T_melts=np.arange(400,4001,400)*stommel.year,
+#                      A_polars = np.arange(0,7.1,.5))
+
+#%% Run experiments with real observations.
+
+filepath = os.path.join(stommel.fig_dir,'tipping3.pkl')
+with open(filepath,'rb') as stream:
+    data = dill.load(stream)
+
+#%% 
+
+plt.close('all')
+P = plotter.TippingPlot(data)
+P.plot(filename='fig12_tipping_percentage')
+