Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Some question about FlowAccumulator, StreamPowerEroder and LinearDiffuser #170

Open
kookil opened this issue Apr 1, 2022 · 5 comments
Open

Comments

@kookil
Copy link

kookil commented Apr 1, 2022

I am using real terrain and precipitation to calculate soil erosion and sediment transport by landlab. And here is part of my code.
(mg, z) = read_esri_ascii(path, name="topographic__elevation") fr = FlowAccumulator(mg) Q = mg.at_node['surface_water__discharge'] ld = LinearDiffuser(mg, linear_diffusivity=0.0001) sp = StreamPowerEroder(mg, K_sp=1., sp_type='Unit', a_sp=1., b_sp=0.5, c_sp=1., use_Q=Q)
The problem I'm having is as follows:

  1. Because I want to change the runoff, so I add
    runoff_rate = np.arange(mg.number_of_nodes, dtype=float) rnff = mg.add_field(" water__unit_flux_in ", runoff_rate, at=" node", clobber=True)
    like the document, but I find the result mg.at_node['surface_water__discharge'] does not change, and I can't find the reason;

  2. I follow the method StreamPowerEroder in landlab to add existing precipitation values in use_Q, but the code reports an error that there is no use_Q in the function, is this a version problem? my landlab version is 2.2.0;

  3. The grid data I am using has several thousand rows and columns, and I run LinerDiffuser for many days without results, is this because the number of data is too large?

I would be grateful if you can give any assistance.

@mcflugen
Copy link
Member

mcflugen commented Apr 1, 2022

@kookil Thank you for your interest in landlab!

  1. Could you please provide a link to the document you are referring to? In your sample code you use " water__unit_flux_in " and " node". Are those extra spaces just typos introduced when you made this issue or are they also in the code you are trying to run? The spaces would definitely result in unexpected behavior.
  2. It looks like the use_Q keyword was removed from the StreamPowerEroder several years ago when we moved from version 1 to 2. I wonder if you are following old documentation?
  3. Yes, it could definitely be the case that your grid is too large for the LinearDiffuser. If this poses a problem for you, we could look into speeding it up.

@kookil
Copy link
Author

kookil commented Apr 4, 2022 via email

@mcflugen
Copy link
Member

mcflugen commented Apr 5, 2022

@kookil In the future, if possible, please reply to GitHub issues through the website rather than by email as GitHub doesn't properly format the emails.

I'm having trouble reproducing your problem. If I run the following code, which I think is what you are also using, the resulting mg.at_node["topographic__elevation"] is different from its initial value. If I then rerun the code with a different random value of lst, I get a different value for the elevations.

import numpy as np
from numpy import random
from landlab import RasterModelGrid
from landlab.components import FlowAccumulator, StreamPowerEroder

mg = RasterModelGrid((9, 9))
mg.add_field('topographic__elevation', mg.node_x**0.5 + mg.node_y**2, at='node')

for edge in (mg.nodes_at_left_edge, mg.nodes_at_right_edge):
    mg.status_at_node[edge] = mg.BC_NODE_IS_FIXED_VALUE
for edge in (mg.nodes_at_top_edge, mg.nodes_at_bottom_edge):
    mg.status_at_node[edge] = mg.BC_NODE_IS_CLOSED

lst = random.random(size=(9,9))
mg.at_node['water__unit_flux_in'] = lst
fr = FlowAccumulator(mg, flow_director='D8')

z_initial = mg.at_node["topographic__elevation"].copy()

sp = StreamPowerEroder(
    mg,
    K_sp=1.,
    sp_type='Unit',
    a_sp=1.,
    b_sp=0.5,
    c_sp=1.,
    discharge_field='surface_water__discharge'
)

fr.run_one_step()
sp.run_one_step(1.)

assert not np.allcose(z_initial, mg.at_node["topographic__elevation"])

Are you doing something differently? Or perhaps I misunderstand what you are expecting?

@kookil
Copy link
Author

kookil commented Apr 7, 2022

I'm sorry for the inconvenience caused by replying by mail earlier.
And I'm surprised that you run my code and get different elevation values ​​because that's exactly what I want. Maybe I'm doing something wrong somewhere else. Thanks for being able to solve my problem!

@mcflugen
Copy link
Member

mcflugen commented Apr 7, 2022

@kookil Could you please paste a minimal code example that demonstrates the issue you are seeing?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants