underworldcode / underworld2

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

Issue Associated with (underworldcode/UW2/cylindrical:cylindrical) image annulus models #379

Closed Prasanna2989 closed 5 years ago

Prasanna2989 commented 5 years ago

Associated with (underworldcode/UW2/cylindrical:cylindrical) image annulus models

I was trying to define various functions for different layers of annulus model using the swarm particles. MantleIndex, lowerCrustIndex...represent the indexes given for swarm particles. I want to apply different buoyancy functions and viscosity functions for different layers of the mantle convection model of Earth. The following method work in the cartesian models. However, when I apply the same mechanism in annulus model the attached error is received. Any suggestions to fix this issue will be appreciated.

Screen Shot 2019-05-12 at 4 21 31 pm

MantleIndex = 0 lowerCrustIndex = 1 upperCrustIndex = 2 lowerMantleIndex = 3 stickyAirIndex = 4 lithoIndex = 5 continentalCrustIndex = 6

####################### E = 23.03 ArrFunction = viscosity_M uw.function.math.exp((E/(tField +1.))-(E/(delta_T_M +1.))) d_eta_UM = 1. creep_UM = d_eta_UMArrFunction

Similar method for creep_LC,.....

#######################

viscosityMap = { MantleIndex : creep_UM,
lowerCrustIndex : creep_LC, upperCrustIndex : creep_UC, lithoIndex : creep_Lth, stickyAirIndex : 1e4,#creep_diffusion/10. lowerMantleIndex : creep_LM, continentalCrustIndex : creep_CC} viscosityFn = fn.branching.map( fn_key = materialVariable, mapping = viscosityMap )

bodyForceFn = g tField Ra / (annulus.radialLengths[1]-annulus.radialLengths[0])

densityMap = { MantleIndex : bodyForceFn, lowerCrustIndex : bodyForceFn, upperCrustIndex : bodyForceFn, lithoIndex : bodyForceFn, stickyAirIndex : bodyForceFn,#Need to review lowerMantleIndex : bodyForceFn, continentalCrustIndex : bodyForceFn }

g = 1.0*annulus.fn_unitvec_radial() densityFn = fn.branching.map( fn_key = materialVariable, mapping = densityMap )

stokesSLE = uw.systems.Stokes( vField, pField, fn_viscosity=viscosityFn, fnbodyforce=bodyForceFn, conditions=vBC, removeBCs=False)

Thanks Pras

jmansour commented 5 years ago

Hi @Prasanna2989

This is either a bug in our code, or your materialVariable swarm variable has not been configured correctly.

So with the fn.branching.map functions, your key function will be evaluated first, and its results will determine which function to evaluate to get the required result (ie, the mapping). In your case, the key function is returning a strange number (2342143069), whereas going by what you've said above it should take one of [0,1,2,3,4,5,6]. So it's possibly using an uninitialised value.

I'd suggest the following:

  1. Init your materialVariable to some known value.
    materialVariable.data[:] = 99999
  2. Then go ahead and configure materialVariable as normal.
  3. Run the rest of the script and check the error message.

If it reports key value (99999) without a mapped function, then you know you need to reconsider step 2 and check how you're missing particles. If on the other hand it reports another random number, it may be a bug in our API, in which case I'd suggest you post your model (but only after getting rid of everything unnecessary).

Prasanna2989 commented 5 years ago

Thanks jmansour.

I did not add the line of initialising. The initialising lines as follows.

MantleIndex = 0 materialVariable.data[:] = MantleIndex

still I get the same error.

I tried with following simple mapping and still have the same error.

g = 1.0 bodyForceFn = g tField Ra / (annulus.radialLengths[1]-annulus.radialLengths[0]) densityMap = {0:bodyForceFn} densityFn = fn.branching.map( fn_key = materialVariable, mapping = densityMap ) buoyancyFn = densityFn * annulus.fn_unitvec_radial()

viscosityMap = {0:ArrFunction} viscosityFn = fn.branching.map( fn_key = materialVariable, mapping = viscosityMap )

stokesSLE = uw.systems.Stokes(vField, pField, fn_viscosity=viscosityFn, fn_bodyforce=buoyancyFn, conditions=vBC, _removeBCs=False)

jmansour commented 5 years ago

Ok, can you post a simplified script? Please remove everything that is unnecessary as this will allow us to debug the problem much quicker.

Prasanna2989 commented 5 years ago

When I save materialIndex variable, it save the swarm with required material indexes. It is only density and viscosity mapping seems the issue.

jmansour commented 5 years ago

Do any map functions work?

Prasanna2989 commented 5 years ago

My bad. I have not initialised the material before send to the solver. Now it works. Appreciate your suggestions jmansour.

jmansour commented 5 years ago

Ok, easy fix, great!