Remove integrated stress in lithosphere

The way it is computed is too specific of one problem to be useful
in the general case.

Author: amorison

stagpy/config.py

@@ -217,8 +217,6 @@ def _index_collection(arg):
                        True, 'min stress in plots')),
     ('stressmax', Conf(800, False, None, {},
                        True, 'max stress in plots')),
-    ('lstressmax', Conf(50, False, None, {},
-                        True, 'max lithospheric stress in plots')),
 ))
 
 CONF_DEF['info'] = OrderedDict((

stagpy/plates.py

@@ -448,90 +448,6 @@ def plot_plates(step, time, vrms_surface, trench, ridge, agetrench,
         ))
 
 
-def lithospheric_stress(step, trench, ridge, time):
-    """Calculate stress in the lithosphere."""
-    timestep = step.isnap
-    base_lith = step.geom.rcmb + 1 - 0.105
-
-    stressfld = step.fields['sII'].values[0, :, :, 0]
-    r_centers = np.outer(np.ones(stressfld.shape[0]), step.geom.r_centers)
-    stressfld = np.ma.masked_where(r_centers < base_lith, stressfld)
-
-    # stress integration in the lithosphere
-    dzm = (step.geom.r_centers[1:] - step.geom.r_centers[:-1])
-    stress_lith = np.sum((stressfld[:, 1:] * dzm.T), axis=1)
-    ph_coord = step.geom.p_centers  # probably doesn't need alias
-
-    # plot stress in the lithosphere
-    fig, axis, _, _ = field.plot_scalar(step, 'sII', stressfld,
-                                        cmap='plasma_r', vmin=0, vmax=300)
-    # Annotation with time and step
-    axis.text(1., 0.9, str(round(time, 0)) + ' My', transform=axis.transAxes)
-    axis.text(1., 0.1, str(timestep), transform=axis.transAxes)
-    saveplot(fig, 'lith', timestep)
-
-    # velocity
-    vphi = step.fields['v2'].values[0, :, :, 0]
-
-    # position of continents
-    concfld = step.fields['c'].values[0, :, :, 0]
-    _, indcont = _isurf_icont(step)
-    if step.sdat.par['boundaries']['air_layer'] and\
-            not step.sdat.par['continents']['proterozoic_belts']:
-        continents = np.ma.masked_where(
-            np.logical_or(concfld[:-1, indcont] < 3,
-                          concfld[:-1, indcont] > 4),
-            concfld[:-1, indcont])
-    elif step.sdat.par['boundaries']['air_layer'] and\
-            step.sdat.par['continents']['proterozoic_belts']:
-        continents = np.ma.masked_where(
-            np.logical_or(concfld[:-1, indcont] < 3,
-                          concfld[:-1, indcont] > 5),
-            concfld[:-1, indcont])
-    elif step.sdat.par['tracersin']['tracers_weakcrust']:
-        continents = np.ma.masked_where(
-            concfld[:-1, indcont] < 3, concfld[:-1, indcont])
-    else:
-        continents = np.ma.masked_where(
-            concfld[:-1, indcont] < 2, concfld[:-1, indcont])
-
-    # masked array, only continents are true
-    continentsall = continents / continents
-
-    # plot integrated stress in the lithosphere
-    fig0, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
-    ax1.plot(step.geom.p_walls, vphi[:, -1], label='Vel')
-    ax1.axhline(y=0, xmin=0, xmax=2 * np.pi,
-                color='black', ls='solid', alpha=0.2)
-    ax1.set_ylabel("Velocity")
-    ax1.text(0.95, 1.07, str(round(time, 0)) + ' My',
-             transform=ax1.transAxes)
-    ax1.text(0.01, 1.07, str(round(step.time, 8)),
-             transform=ax1.transAxes)
-
-    intstr_scale = step.sdat.scales.stress * step.sdat.scales.length / 1.e12
-    ax2.plot(ph_coord, stress_lith * intstr_scale, color='k', label='Stress')
-    ax2.set_ylabel(r"Integrated stress [$TN\,m^{-1}$]")
-
-    plot_plate_limits(ax1, ridge, trench, conf.plates.vmin,
-                      conf.plates.vmax)
-    plot_plate_limits(ax2, ridge, trench, conf.plates.stressmin,
-                      conf.plates.lstressmax)
-    ax1.set_xlim(0, 2 * np.pi)
-    ax1.set_title(timestep)
-
-    ax1.fill_between(
-        ph_coord, continentsall * conf.plates.vmin,
-        conf.plates.vmax, facecolor='#8b6914', alpha=0.2)
-    ax1.set_ylim(conf.plates.vmin, conf.plates.vmax)
-    ax2.fill_between(
-        ph_coord, continentsall * conf.plates.stressmin,
-        conf.plates.lstressmax, facecolor='#8b6914', alpha=0.2)
-    ax2.set_ylim(conf.plates.stressmin, conf.plates.lstressmax)
-
-    saveplot(fig0, 'svelslith', timestep)
-
-
 def set_of_vars(arg_plot):
     """Build set of needed variables.