UW function that changes along Model.x as velocity BC

[JoeIbrahim]

Hello,

I'm trying to implement a variable boundary condition at the base of the model (that's supposed to simulate doming kinematically) that changes along the X-axis.

Conditionsbottom = fn.branching.conditional([(Model.x > GEO.nd(20. * u.kilometers), GEO.nd(0. * u.centimeter / u.year)),
              (True, GEO.nd(1.-(0.00308505*Model.x**2) + (1.58626*10.**-6. * Model.x**4.)))])

Model.set_velocityBCs(bottom=[None, Conditionsbottom],
                      left=[0.0 * u.centimeter / u.year, None],
                      right=[0.0 * u.centimeter / u.year, 0.0 * u.centimeter / u.year],
                      top=[None, None])

Ideally using this function the velocity field at the base along the Y-axis would look like this:

Though what I do end up getting is a Vy of 1cm/yr at Model.x < 20 km instead of the function:

Do you know where I might have gone wrong?

Thanks for your help!!

[JoeIbrahim]

I tried a simpler boundary condition along the base where after 20km along the X axis I impose positive Model.x BC and before 20km its -Model.x :

Conditionsbottom = fn.branching.conditional([(Model.x > GEO.nd(20. * u.kilometers), GEO.nd(Model.x)),
(True, GEO.nd(-Model.x))])

# fn_condition = fn.branching.conditional(conditions)

Model.set_velocityBCs(bottom=[None, Conditionsbottom],
                      left=[0.0 * u.centimeter / u.year, None],
                      right=[0.0 * u.centimeter / u.year, 0.0 * u.centimeter / u.year],
                      top=[None, None])

and I got a Vy profile along the base that looks like this:

I'm not sure how the scaling with Model.x works, is it being rescaled at some point?

I've also tried it without non dimensionalising Model.x, and got the same results

Conditionsbottom = fn.branching.conditional([(Model.x > GEO.nd(20. * u.kilometers), Model.x), (True,-Model.x)])

Thanks again!

[julesghub]

Hi Joe, Model.x is a uw.fn.

t1 = Model.x
t2 = GEO.nd(Model.x)
print(type(t1),type(t2))

As such it isn't "lazy executed" under GEO.nd() in the conditional for the velocityBCs.

Can you instead build a standalone numpy array using your empirical formula and use GEO.dimensionalise(Model.mesh.data[:,0], u.km).magnitude to build a vector of x values. Input this standalone numpy array into the conditional instead of the GEO.nd(Model.x)

[JoeIbrahim]

Hi Julian,

Thanks for the response!

So I've created a numpy array using my formula and the x values using:

mycosinfunction=np.array([1.-(0.00308505*i**2) + (1.58626*10.**-6. * i**4.) for i in GEO.dimensionalise(Model.mesh.data[:,0], u.km).magnitude ])

If I try to put that into the boundary condition like this:

Conditionsbottom = fn.branching.conditional([(Model.x > GEO.nd(20. * u.kilometers), GEO.nd(0. * u.centimeter / u.year)), (True, mycosinfunction)])

I get the error:

RuntimeError: Issue utilising function of class 'conditional' constructed at: Line 5 of notebook cell 63: (True, mycosinfunction)]) Error message: Issue with clause 0 of conditional function. Consequent function appears to return result of size 1. Previous function returned result of size 9849. All consequent function must return results of identical size

I might've misunderstood something, but would I have to include the first part of the conditional (Model.x > 20km, 0. Velocity) into the array as well? Also is it okay that the x values span the entire model, or should I only grab one row - so something like(Model.mesh.data[0:200,0], u.km)in my case?

[julesghub]

Ok the trick here is to use your existing formula with a substitution.

Make

real_x = GEO.dimensionalise(1, u.km).magnitude*Model.x

And swap all Model.x for real_x in the following

Conditionsbottom = fn.branching.conditional([(Model.x > GEO.nd(20. * u.kilometers), GEO.nd(0. * u.centimeter / u.year)),
              (True, GEO.nd(1.-(0.00308505*Model.x**2) + (1.58626*10.**-6. * Model.x**4.)))])

My previous answer about the standalone numpy array was wrong, thanks for trying it. Give the above a go and let me know it goes.

[JoeIbrahim]

Perfect that did the trick!

Thanks for your help Julian!

[julesghub]

@JoeIbrahim thinking about this some more the use of the "trick" will contain real world values in the numerics (ie, unscaled). This is bad. A better way forward is NonDimensionalise all the coefficients in the forumla rather than Dimensionalise Model.x.

So (0.00308505*Model.x**2) would become GEO.nd(0.00308505 * u.km*-2)Model.x**2.

Also there's no need to wrap the entire formula in GEO.nd().

Let me know how this goes.

[JoeIbrahim]

I gave this way a go and it works too. I'll keep non-dimensionalising the coefficients moving forward. Thanks for all the help!