Skip to content
Snippets Groups Projects
Select Git revision
  • c01c699f16926b6fb0c9a810c20ae700a0988bb2
  • main default protected
  • YL
  • NS
4 results

postprocessing.py

Blame
    • plot_mass_point_model.py 2.92 KiB
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107 
    

    

    

    """

Mass Point Model
================
"""

# %% md
# 
# This example shows how to simulate the movement of a masspoint on a topography
# using :class:`.MassPointModel`.
#
# 

# %% md
# 
# First, import the class :class:`.MassPointModel` and create an instance by

from psimpy.simulator import MassPointModel

mpm = MassPointModel()

# %% md
# 
# Required inputs for a simulation using :class:`.MassPointModel` include:
#     1. topographic data: digital elevation model (in `ESRI ascii` format)
#     2. friction coefficients: coulomb friction and turbulent friction coefficients
#     3. initial state: initial location and initial velocity of the masspoint
#     4. computational parameters: such as time step, end time, etc.
#
# The synthetic topography ``synthetic_topo.asc`` is used here for illustration.
# It is located at the `tests/data/` folder.

import os
import linecache
import numpy as np

dir_data = os.path.abspath('../../../tests/data/')
elevation = os.path.join(dir_data, 'synthetic_topo.asc')


# %% md
# .. note:: You may need to modify ``dir_data`` according to where you save
#    ``synthetic_topo.asc`` on your local machine.

# %% md
# We can load the elevation data and visulize it. The figure below shows how
# the topography looks like, as well as the initial location of the masspoint
# (noted by the red dot).

header = [linecache.getline(elevation, i) for i in range(1,6)]
header_values = [float(h.split()[-1].strip()) for h in header]
ncols, nrows, xll, yll, cellsize = header_values
ncols = int(ncols)
nrows = int(nrows)

x_values = np.arange(xll, xll+(cellsize*ncols), cellsize)
y_values = np.arange(yll, yll+(cellsize*nrows), cellsize)
        
z_values = np.loadtxt(elevation, skiprows=6)
z_values = np.rot90(np.transpose(z_values))

# initial location
x0 = 200
y0 = 2000
z0 = z_values[0, int(x0/cellsize)]

import matplotlib.pyplot as plt

fig, ax = plt.subplots(2, 1, figsize=(10,6), height_ratios=[3,2])

fig0 = ax[0].contourf(x_values, y_values, z_values, levels=20)
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')
ax[0].set_title('synthetic topography')
ax[0].scatter(x0, y0, s=10, c='r', marker='o')
cbar = plt.colorbar(fig0, ax=ax[0], format='%d', orientation='horizontal',
    fraction=0.1, pad=0.2)
cbar.ax.set_ylabel('z')


ax[1].plot(x_values, z_values[0, :])
ax[1].scatter(x0, z0, s=10, c='r', marker='o')
ax[1].set_xlabel('x')
ax[1].set_ylabel('z')
ax[1].set_xlim(0, 5000)
ax[1].set_title('cross section at any y')

plt.tight_layout()

# %% md
# We set the friction coefficients as

mu = 0.15
xi = 1000

# %% md
# Given above topography, initial location, and friction coefficients, we can
# call the :py:meth:`.MassPointModel.run` method to perform a simulation. Other
# parameters are set to their default values. (we suppress raised warnings )
# (other parameters are set to their default values).

import warnings

warnings.filterwarnings("ignore")
output = mpm.run(elevation=elevation, coulomb_friction=mu, turbulent_friction=xi, x0=x0, y0=y0)