2D surface processes implementation

[bknight1]

Implementation of surface processes in 2D. (1) erosion and sedimentation rate (2) diffusive surface.

Updated tutorial list with 6.1 (erosion and sedimentation rate) and 6.2 (diffusive surface).

Updated user guide with info on how to include these processes in models.

[julesghub]

Nice work @bknight1. We'll review this week and get back to you.

[rbeucher]

Great! I'm gonna look at that. Thanks @bknight1 !

[bknight1]

Thanks @julesghub and @rbeucher

I realised I pulled the tutorials from v2.10.2, so the Model.init_model will need updating.

I have also just noticed that the way to include passive tracers has changed so I'll have to update some parts of the code in the surface processes file. What's the best way to make these changes?

EDIT: I can see that I just add the updated files to my branch, I'll get those updates done today

[bknight1]

I've updated the passive tracers to be accessed through Model, however when trying to update the z coordinates I get the following error: ValueError: assignment destination is read-only

Example to recreate error:

import numpy as np

npoints = int(Model.length.magnitude)

coords = np.ndarray((npoints, 2))
coords[:, 0] = np.linspace(GEO.nd(Model.minCoord[0]), GEO.nd(Model.maxCoord[0]), npoints)
coords[:, 1] = GEO.nd(uppercrust.top)
Model.add_passive_tracers(name="PT", vertices=coords, advect=False)
Model.PT_tracers.data[:,0] = 0.

Using: UW version - 2.11.0b UWGeodynamics version - 2.11.0-dev-485ad28(master)

[rbeucher]

@NengLu can you have a look at this please?

[NengLu]

@NengLu can you have a look at this please?

Sure, by the way, thanks ben for sharing this.

[NengLu]

ValueError: assignment destination is read-only

Hi, Ben, the PT_tracer should be set like:

Model.PT_tracers = Model.add_passive_tracers(name="PT", vertices=coords, advect=False)
Model.PT_tracers.data[:,0] = 0.

or independently, not the Model's attribute, like:

PT_tracers = Model.add_passive_tracers(name="PT", vertices=coords, advect=False)
PT_tracers.data[:,0] = 0.

And here PT_tracer.data[:,0] is initially set by "add_passive _tracters" be coords[:, 0], then be set to 0. If you want to trace another field, it could be done like:

PT_tracers.add_tracked_field(Model.velocityField,
                              name="tracers_vel",
                              units=u.meter/u.year,
                              dataType="float")

[rbeucher]

No Neng. The passive tracers access method has been updated. You can't do that anymore. @julesghub why do we get a read only error?

Get Outlook for Androidhttps://aka.ms/AAb9ysg

---

From: Neng Lu @.> Sent: Thursday, September 30, 2021 7:47:03 PM To: underworldcode/UWGeodynamics @.> Cc: Romain Beucher @.>; Mention @.> Subject: Re: [underworldcode/UWGeodynamics] 2D surface processes implementation (#249)

ValueError: assignment destination is read-only

Hi, Ben, the PT_tracer should be set like:

Model.PT_tracers = Model.add_passive_tracers(name="PT", vertices=coords, advect=False) Model.PT_tracers.data[:,0] = 0.

or independently, not the Model's attribute, like:

PT_tracers = Model.add_passive_tracers(name="PT", vertices=coords, advect=False) PT_tracers.data[:,0] = 0.

And here PT_tracer.data[:,0] is initially set by "add_passive _tracters" be coords[:, 0], then be set to 0. If you want to tracer another field, it could be done like:

PT_tracers.add_tracked_field(Model.velocityField, name="tracers_vel", units=u.meter/u.year, dataType="float")

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHubhttps://aus01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Funderworldcode%2FUWGeodynamics%2Fpull%2F249%23issuecomment-931157606&data=04%7C01%7Cromain.beucher%40anu.edu.au%7C8c6f8f7d832949a0eeb208d983f74389%7Ce37d725cab5c46249ae5f0533e486437%7C0%7C0%7C637685920275415561%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=xtCJeOXMsIbnuNYFCjPJzBYmfkRSnbiPrlh%2BZ0r3fQw%3D&reserved=0, or unsubscribehttps://aus01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FAA4CVFXXO6FHBNDN5JF7BF3UEQWZPANCNFSM5ELKJCEQ&data=04%7C01%7Cromain.beucher%40anu.edu.au%7C8c6f8f7d832949a0eeb208d983f74389%7Ce37d725cab5c46249ae5f0533e486437%7C0%7C0%7C637685920275415561%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=Tm14OimHH8PCaxW7GC2YbKiUO6ZJpCEfSOsfrzxk%2BAU%3D&reserved=0. Triage notifications on the go with GitHub Mobile for iOShttps://aus01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fapps.apple.com%2Fapp%2Fapple-store%2Fid1477376905%3Fct%3Dnotification-email%26mt%3D8%26pt%3D524675&data=04%7C01%7Cromain.beucher%40anu.edu.au%7C8c6f8f7d832949a0eeb208d983f74389%7Ce37d725cab5c46249ae5f0533e486437%7C0%7C0%7C637685920275425518%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=Xg6zhkMhPtxzVopyIXouoDpVjgKp2G8NLbD8VR%2FRgFU%3D&reserved=0 or Androidhttps://aus01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fplay.google.com%2Fstore%2Fapps%2Fdetails%3Fid%3Dcom.github.android%26referrer%3Dutm_campaign%253Dnotification-email%2526utm_medium%253Demail%2526utm_source%253Dgithub&data=04%7C01%7Cromain.beucher%40anu.edu.au%7C8c6f8f7d832949a0eeb208d983f74389%7Ce37d725cab5c46249ae5f0533e486437%7C0%7C0%7C637685920275425518%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=EYyFLVyIIdvlxbfBXlcAOw%2BCxKDXRrj9pyjsKX%2F0ftY%3D&reserved=0.

[rbeucher]

Tracers are swarms. So to access the coordinates you need to use:

Model.PT_tracers.particleCoordinates.data[:, 0] = 0.

[bknight1]

I have just tried that and getting the same error @rbeucher

'ValueError: assignment destination is read-only'

[julesghub]

I believe the read only error is because the swarm is under a context manager, ie. one needs to use the following to write values to the data structure.

with Model.PT_tracers.deform_swarm():
   Model.PT_tracers.particleCoordinates.data[:,0] = 0.

[bknight1]

Thanks @julesghub that works! Will re-run the tests and update the files today

[rbeucher]

Great @NengLu will go through your examples and then we can merge.

[bknight1]

Should be good to go now @NengLu

[NengLu]

Should be good to go now @NengLu

Hi, Romain and Ben @rbeucher @bknight1, I have gone through the codes and examples, there are small typos in the Tutorial_6.1 and Tutorial_6.2, Model.surfacetracers_tracersshould be Model.surface_tracers

The topography from sedimentation&erosion seems incorrect, may need a further check with velocitySurface_2D(SurfaceProcesses)

[rbeucher]

What's the reference for the elevations?

[bknight1]

I believe this is due to differences in the diffusivity rate (D = 150 m2/yr) and erosion and sedimentation rates (2 mm/yr). That's my fault for leaving in a very high erosion and sedimentation rate, which would flatten the topography very quickly (as is seen in the figure for elevation).

A more appropriate erosion and sedimentation rate should be ~0.15 mm/yr, which should then produce a similar topography to the no erosion and diffusive surface models.

Model.surfaceProcesses = velocitySurface_2D(
    airIndex=air.index,
    sedimentIndex=sediment.index,
    sedimentationRate= 0.15*u.millimeter / u.year,
    erosionRate= 0.15*u.millimeter / u.year,
    surfaceElevation=uppercrust.top,
    surfaceArray = coords
)

@rbeucher the reference height for the erosion & sedimentation rate is 0 km (uppercrust.top), which is why the topography is flattened out at 0 km except for the parts that experience the most uplift (at around x = 120 km and x = 240 km).

[NengLu]

I believe this is due to differences in the diffusivity rate (D = 150 m2/yr) and erosion and sedimentation rates (2 mm/yr). That's my fault for leaving in a very high erosion and sedimentation rate, which would flatten the topography very quickly (as is seen in the figure for elevation).

A more appropriate erosion and sedimentation rate should be ~0.15 mm/yr, which should then produce a similar topography to the no erosion and diffusive surface models.

Model.surfaceProcesses = velocitySurface_2D(
    airIndex=air.index,
    sedimentIndex=sediment.index,
    sedimentationRate= 0.15*u.millimeter / u.year,
    erosionRate= 0.15*u.millimeter / u.year,
    surfaceElevation=uppercrust.top,
    surfaceArray = coords
)

@rbeucher the reference height for the erosion & sedimentation rate is 0 km (uppercrust.top), which is why the topography is flattened out at 0 km except for the parts that experience the most uplift (at around x = 120 km and x = 240 km).

Merge is good to go! @rbeucher

[rbeucher]

Awesome guys. Thanks @bknight1 @NengLu

[bknight1]

Nice one, thanks everyone!

[rbeucher]

@bknight1 can you send me your ORCID so I can add you to the list of contributors?

[bknight1]

It's 0000-0001-7919-2575, thanks @rbeucher