underworldcode / underworld2

underworld2: A parallel, particle-in-cell, finite element code for Geodynamics.
http://www.underworldcode.org/
Other
168 stars 58 forks source link

Surface processes #644

Open PatriceFRey opened 1 year ago

PatriceFRey commented 1 year ago

Hello all, Is there the 3D equivalent to following 2D functions: GEO.surfaceProcesses.velocitySurface_2D GEO.surfaceProcesses.diffusiveSurface_2D ? Thanks

bknight1 commented 1 year ago

Hey @patrice-rey,

I wrote these functions and also implemented a 3D velocity surface too. Let me know if you'd like the functions, I have two, one with a constant erosional and sedimentation rate and another that varies laterally. I think I can implement the diffusive surface in 3D but I haven't got round to it yet. Let me know what you need.

Cheers

PatriceFRey commented 1 year ago

G'day Ben, Great, any 3D function would do really, though I like the idea where the sedimentation rate varies laterally. In a rift environment, this would allow the partitioning of sedimentation along the margins. Although tracking the position of the margins through time may not be that easy. Thanks Ben

bknight1 commented 1 year ago

All good Patrice, let me see what I can do in the next couple of days to make both situations feasible. If you'd like something ASAP I can upload my original implementation.

Cheers, Ben

PatriceFRey commented 1 year ago

No rush Ben, I can wait a few days or weeks. Thanks

bknight1 commented 1 year ago

Okay great, let me see what I can do and I'll get back to you.

Cheers

bknight1 commented 1 year ago

Hey @patrice-rey,

I've managed to put something together, just got to run a final couple of tests on the HPC to make sure it works. I'll share it once completed.

Thanks

PatriceFRey commented 1 year ago

Thanks for the update Ben. I am on annual leave until next Wednesday 7th December, so no rush. Cheers

bknight1 commented 1 year ago

Hey @PatriceFRey sorry it's taken so long to get back to you. I've implemented a 3D surface processes module for you to use, for now you'll have to have the file attached in the same directory as the running script to use the sub module, it'll also need to be renamed so it has a .py extension. You'll then need to set up the surface processes as following in your script:

from UW3DsurfaceProcesses import *

coords = None

if uw.mpi.rank == 0:
    ### set up the x and y coords of the surface
    x = np.linspace(Model.minCoord[0].m, Model.maxCoord[0].m, int((Model.maxCoord[0].m/2)+1))
    y = np.linspace(Model.minCoord[1].m, Model.maxCoord[1].m, int((Model.maxCoord[1].m/2)+1)) 

    xi, yi = np.meshgrid(x, y)

    ### this strips units, so units have to be reasigned after creating the array
    coords = np.zeros(shape=(xi.flatten().shape[0], 3))
    coords[:,0] = xi.flatten()
    coords[:,1] = yi.flatten()
    coords[:,2] = np.zeros_like(coords[:,0])+air.bottom.m

uw.mpi.barrier()

coords = uw.mpi.comm.bcast(coords, root=0)

## add dim after bcast, fails if done before
coords = (coords * u.kilometer)

### set ve by the y  coord (lateral position)

ve_conditions = fn.branching.conditional([((Model.y >= GEO.nd(Model.maxCoord[1]/2.)) , GEO.nd(1. * u.millimeter/u.year)),
                                          (True, GEO.nd(0.01 * u.millimeter/u.year))])

# ### constant vs across the model
vs_conditions = fn.branching.conditional([(True, GEO.nd(0.05 * u.millimeter/u.year))])

Model.surfaceProcesses = velocitySurface3D(airIndex=air.index,
                        sedimentIndex= sediment.index,
                        surfaceArray = coords,                ### coords for the surface
                        vs_condition = vs_conditions,         ### branching condition for erosion rate
                        ve_condition = ve_conditions,         ### branching condition for sedimentation rate
                        surfaceElevation=air.bottom,
                        method='linear'
                    )

Let me know if you run into any issues.

Cheers

UW3DsurfaceProcesses.py.txt

PatriceFRey commented 1 year ago

Hello Ben, thanks a lot for your time doing this. I'll run a test over the next few days and keep you posted. Cheers

bknight1 commented 1 year ago

All good Patrice, looking forward to getting updates! Cheers

julesghub commented 1 year ago

@bknight1 - should I add the above method to UW2. Are there any tests for it?

bknight1 commented 1 year ago

Hey @PatriceFRey, how did the code work for you?