Skip to content

Commit

Permalink
Veg growth funcs
Browse files Browse the repository at this point in the history
  • Loading branch information
reevesi_laptop committed Oct 25, 2022
1 parent a768e71 commit 6fd487a
Show file tree
Hide file tree
Showing 2 changed files with 180 additions and 15 deletions.
89 changes: 85 additions & 4 deletions dubeveg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 16,7 @@
import math
import random
import matplotlib.pyplot as plt
import time

# __________________________________________________________________________________________________________________________________
# USER INPUTS
Expand Down Expand Up @@ -183,14 184,18 @@
landward_transport = np.zeros([len(timeits)])
avalanches = np.zeros([len(timeits)])

start_time = time.time() # Record time at start of simulation to track duration


# __________________________________________________________________________________________________________________________________
# MAIN ITERATION LOOP

for it in range(iterations):
year = math.ceil(it / iterations_per_cycle)

# Print time step to screen
print("\r", "Time Step: ", (it / iterations_per_cycle), end="")
print("\r", "Time Step: ", (it / iterations_per_cycle), "yrs", end="")

year = math.ceil(it / iterations_per_cycle)

if eqtopo_initial.ndim == 3:
eqtopo = np.squeeze(eqtopo_initial[:, :, it]) / slabheight_m
Expand Down Expand Up @@ -234,6 239,11 @@
balance = balance (topo - before) # Update the sedimentation balance map
stability = stability abs(topo - before)


# ADD ADDITIONAL DATA SAVING



# --------------------------------------
# BEACH UPDATE

Expand All @@ -260,21 270,92 @@
pcurr,
)

seainput = topo - before1 # Sand added to the beach by the sea
# seainput_slabs[it] = np.nansum(seainput) # IRBR 24OCt22: Bug here. Why is this even needed?

if 'seainput_total' in locals():
seainput_total[:, :, it] = seainput

if 'diss_total' in locals():
diss_total[:, :, it] = diss

if 'cumdiss_total' in locals():
cumdiss_total[:, :, it] = cumdiss

if 'pwave_total' in locals():
pwave_total[:, :, it] = pwave

if 'pbeachupdate_total' in locals():
pbeachupdate_total[:, :, it] = pbeachupdate

if 'seainput_sum' in locals():
seainput_sum = seainput_sum seainput

inundated = np.logical_or(inundated, inundatedold) # Combine updated area from individual loops
pbeachupdatecum = pbeachupdate pbeachupdatecum # Cumulative beachupdate probabilities
topo = routine.enforceslopes3(topo, veg, slabheight, repose_bare, repose_veg, repose_threshold)[0] # Enforce angles of repose again
balance = balance (topo - before1)
stability = stability abs(topo - before1)

if 'balanceb_sum' in locals():
balanceb_sum = balanceb_sum balance

if 'balanceb_tot' in locals():
balanceb_tot[:, :, it] = balance

# Temp Plot
if 'stabilityb_sum' in locals():
stabilityb_sum = stabilityb_sum stability

if 'stabilityb' in locals():
stabilityb[:, :, it] = balance

beachcount = 1 # Update counter


# --------------------------------------
# VEGETATION

if it % iterations_per_cycle == 0 and it > 0: # Update the vegetation

veg_multiplier = (1 growth_reduction_timeseries[vegcount]) # For the long term reduction.
sp1_peak = sp1_peak_at0 * veg_multiplier
sp2_peak = sp2_peak_at0 * veg_multiplier
spec1_old = spec1
spec2_old = spec2
spec1 = routine.growthfunction1_sens(spec1, balance * slabheight_m, sp1_a, sp1_b, sp1_c, sp1_d, sp1_e, sp1_peak)
spec2 = routine.growthfunction2_sens(spec2, balance * slabheight_m, sp2_a, sp2_b, sp2_d, sp1_e, sp2_peak)

# Lateral Expansion
lateral1 = routine.lateral_expansion(spec1_old, 1, lateral_probability * veg_multiplier) # Species 1








# --------------------------------------
# RESET DOMAINS
balance_copy = balance
balance[:] = 0 # Reset the balance map
inundated[:] = 0 # Reset part of beach with wave/current action
pbeachupdatecum[:] = 0
stability[:] = 0 # Reset the balance map


# __________________________________________________________________________________________________________________________________
# SIMULATION END

# Print elapsed time of simulation
print()
SimDuration = time.time() - start_time
print()
print("Elapsed Time: ", SimDuration, "sec")

# Temp Plot
plt.matshow(topo * slabheight,
cmap='terrain',
)
plt.title("Elev TMAX")

plt.show()
106 changes: 95 additions & 11 deletions routines_dubeveg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -315,7 315,9 @@ def marine_processes3_diss3e(total_tide, msl, slabheight, cellsizef, topof, eqto
L = max(0, -30.59 46.74 * tide_m)

# Runup as a function of wave conditions (H, L) and foreshore slope (b)
if b / math.sqrt(H / L) < 0.3:
if L == 0 and H > 0:
runup_m = 0.043 * math.sqrt(H * L)
elif (L > 0 and H > 0) and b / math.sqrt(H / L) < 0.3:
runup_m = 0.043 * math.sqrt(H * L)
else:
runup_m = 1.1 * (0.35 * math.tan(b) * math.sqrt(H * L) math.sqrt(H * L * (0.563 * math.tan(b**2) 0.004)) / 2)
Expand All @@ -331,8 333,10 @@ def marine_processes3_diss3e(total_tide, msl, slabheight, cellsizef, topof, eqto
toolow = topof < totalwater # [0 1] Give matrix with cells that are potentially under water
pexposed = np.ones(topof.shape) # Initialise matrix
for m20 in range(len(topof[:, 0])): # Run along all the rows
m21 = np.argwhere(topof[m20, :] >= totalwater)[0][0]
pexposed[m20, m21: -1] = 1 - shelterf
twlloc = np.argwhere(topof[m20, :] >= totalwater)
if len(twlloc) > 0: # If there are any sheltered cells
m21 = twlloc[0][0] # Finds for every row the first instance where the topography exceeds the sea level
pexposed[m20, m21: -1] = 1 - shelterf # Subtract shelter from exposedness: sheltered cells are less exposed

# --------------------------------------
# FILL TOPOGRAPHY TO EQUILIBRIUM
Expand Down Expand Up @@ -421,14 425,94 @@ def marine_processes3_diss3e(total_tide, msl, slabheight, cellsizef, topof, eqto
return topof, inundatedf, pbeachupdate, diss, cumdiss, pwave










def growthfunction1_sens(spec, sed, A1, B1, C1, D1, E1, P1):
"""Input is the vegetation map, and the sedimentation balance in units of cell dimension (!), i.e. already adjusted by the slabheight
Output is the change in vegetation effectiveness
Ammophilia-like vegetation"""

# Physiological range (needs to be specified)
minimum = 0.0
maximum = 1.0

# Vertices (these need to be specified)
x1 = A1 # was x1 = -1.4
x2 = B1 # was x2 = 0.1
x3 = C1 # was x3 = 0.45
x4 = D1 # was x4 = 0.85
x5 = E1 # was x5 = 1.4
y1 = -1.0 # y1 = -1.0; y1 = -1.0
y2 = 0.0 # y2 = 0.0; y2 = 0.0
y3 = P1 # y3 = 0.4; y3 = P1;
y4 = 0.0 # y4 = 0.0; y4 = 0.0
y5 = -1.0 # y5 = -1.0; y5 = -1.0

# Slopes between vertices (calculated from vertices)
s12 = (y2 - y1) / (x2 - x1)
s23 = (y3 - y2) / (x3 - x2)
s34 = (y4 - y3) / (x4 - x3)
s45 = (y5 - y4) / (x5 - x4)

leftextension = (sed < x1) * -1
firstleg = (sed >= x1) * (sed < x2) * ((sed - x1) * s12 y1)
secondleg = (sed >= x2) * (sed < x3) * ((sed - x2) * s23 y2)
thirdleg = (sed >= x3) * (sed < x4) * ((sed - x3) * s34 y3)
fourthleg = (sed >= x4) * (sed < x5) * ((sed - x4) * s45 y4)
rightextension = (sed >= x5) * -1

spec = spec leftextension firstleg secondleg thirdleg fourthleg rightextension

spec[spec < minimum] = minimum
spec[spec > maximum] = maximum

return spec


def growthfunction2_sens(spec, sed, A2, B2, D2, E2, P2):
"""Input is the vegetation map, and the sedimentation balance in units of cell dimension (!), i.e. already adjusted by the slabheight
Output is the change in vegetation effectiveness
Buckthorn-type vegetation"""

# Physiological range (needs to be specified)
minimum = 0.0 # was -0.5
maximum = 1.0 # was 1.5

# Vertices (these need to be specified)
x1 = A2 # was x1 = -1.3
x2 = B2 # was x2 = -0.65
x3 = 0 # was x3 = 0
x4 = D2 # was x4 = 0.2
x5 = E2 # was x5 = 2.2
y1 = -1.0 # y1 = -1.0
y2 = 0.0 # y2 = 0.0
y3 = P2 # y3 = 0.1
y4 = 0.0 # y4 = 0.0
y5 = -1.0 # y5 = -1.0

# Slopes between vertices (calculated from vertices)
s12 = (y2 - y1) / (x2 - x1)
s23 = (y3 - y2) / (x3 - x2)
s34 = (y4 - y3) / (x4 - x3)
s45 = (y5 - y4) / (x5 - x4)

leftextension = (sed < x1) * -1
firstleg = (sed >= x1) * (sed < x2) * ((sed - x1) * s12 y1)
secondleg = (sed >= x2) * (sed < x3) * ((sed - x2) * s23 y2)
thirdleg = (sed >= x3) * (sed < x4) * ((sed - x3) * s34 y3)
fourthleg = (sed >= x4) * (sed < x5) * ((sed - x4) * s45 y4)
rightextension = (sed >= x5) * -1

spec = spec leftextension firstleg secondleg thirdleg fourthleg rightextension

spec[spec < minimum] = minimum
spec[spec > maximum] = maximum

return spec


def lateral_expansion(veg, dist, prob):
"""LATERAL_EXPANSION implements lateral expansion of existing vegetation patches
1) mark cells that lie within <dist> of existing patches: probability for new vegetated cells = 1
2) cells not adjacent to existing patches get probability depending on elevation: pioneer most likely to establish between 3 and 5 m sea level"""



Expand Down

0 comments on commit 6fd487a

Please sign in to comment.