error in my BCs

[elisafierro]

Hello Underworld2 team,

I had the following error after running my configuration file and I can't go on with the modeling. Could you help me? Could I send you by e-mail my configuration file?

Thanks!! Elisa

/root/underworld2/underworld/systems/_bsscr.py:449: UserWarning: A petsc error has been detected during the solve.
The resultant solution fields are most likely erroneous, check them thoroughly.
This is possibly due to many things, for example solar flares or insufficient boundary conditions.
The resultant KSPConvergedReasons are (f_hat, outer, backsolve) (0,-3,0).
  warnings.warn(estring)

ValueErrorTraceback (most recent call last)
<ipython-input-29-1b2a209e8394> in <module>()
    122                     v = velocityField.data[index]
    123                 mesh.data[count][0] = mesh.data[count][0] + dt * drivingV(mesh.data[count][0], v, 0., v * time + 4 * res)
--> 124                 count += 1
    125 
    126 

/usr/lib/python2.7/contextlib.pyc in __exit__(self, type, value, traceback)
     33                 value = type()
     34             try:
---> 35                 self.gen.throw(type, value, traceback)
     36                 raise RuntimeError("generator didn't stop after throw()")
     37             except StopIteration, exc:

/root/underworld2/underworld/mesh/_mesh.pyc in deform_mesh(self, isRegular, remainsRegular)
    254                                                      +"Note that any submesh should not be modified directly, "
    255                                                      +"but will instead be updated automatically on return from "
--> 256                                                      +"the 'deform_mesh' context manager. \nEncountered exception message:\n")
    257         finally:
    258             self.data.flags.writeable = False

/root/underworld2/underworld/mesh/_mesh.pyc in deform_mesh(self, isRegular, remainsRegular)
    247         self.data.flags.writeable = True
    248         try:
--> 249             yield
    250         except Exception as e:
    251             raise uw._prepend_message_to_exception(e, "An exception was raised during mesh deformation. "

<ipython-input-29-1b2a209e8394> in <module>()
    121                 else:
    122                     v = velocityField.data[index]
--> 123                 mesh.data[count][0] = mesh.data[count][0] + dt * drivingV(mesh.data[count][0], v, 0., v * time + 4 * res)
    124                 count += 1
    125 

ValueError: An exception was raised during mesh deformation. Your mesh may only be partially deformed. You can reset your mesh using the 'reset' method. Note that any submesh should not be modified directly, but will instead be updated automatically on return from the 'deform_mesh' context manager. 
Encountered exception message:

setting an array element with a sequence.

[rbeucher]

Hi Elisa,

Looks like the pb comes from the drivingV function which I suppose is returning a list or a tuple.Can you post it so we can have a look? v is a tuple of dimension mesh.dim and is used as an argument to your drivingV function. Could that be that you forgot to take a specific dimension of v in your calculation?

[elisafierro]

Hi, I posted 1) my velocity BCs, 2) the drivingV function and 3) the mesh deformation 'function' inside the 'timestepping loop'. Thanks so much!!! Elisa

1)

#VelocityBCs

vxWalls = mesh.specialSets["MinI_VertexSet"] + \
          mesh.specialSets["MaxI_VertexSet"] + \
          mesh.specialSets["MinJ_VertexSet"]

vyWalls = mesh.specialSets["MinJ_VertexSet"] + \
          mesh.specialSets["MaxJ_VertexSet"] + \
          mesh.specialSets["MaxI_VertexSet"]

max_i = 0

V_ff = 1.5 / velocity_scaling_cmyr
vxLeftWall  = 0.
vyTopWall = 0.
for i in range (0, res):
    index = 4 * res * (i + 1) - 1
    if temperatureField.data[index] < 1.:
        max_i = i - 1
        break
delta = V_ff / max_i
v = 0.

for index in mesh.specialSets["MinI_VertexSet"]:
    if temperatureField.data[index] < 1.:
        velocityField.data[index] = [0.,0.]
    else:
        velocityField.data[index,0] = vxLeftWall

for index in mesh.specialSets["MinJ_VertexSet"]:
    velocityField.data[index] = [0.,0.]
for index in mesh.specialSets["MaxJ_VertexSet"]:
    velocityField.data[index,1] = vyTopWall
for index in mesh.specialSets["MaxI_VertexSet"]:
    if temperatureField.data[index] < 1.:
        velocityField.data[index] = [V_ff,0.]
    else:
        velocityField.data[index] = [v,0.]
        v += delta

velocityBC = uw.conditions.DirichletCondition( variable = velocityField, 
                                               indexSetsPerDof = (vxWalls, vyWalls) )`

2)

#driving V function

def drivingV(x, v, minX, maxX):
    m = (v - vxLeftWall) / (maxX-minX) 
    return m * (x - minX) + vxLeftWall  `

3)

#perform timestepping (just a part of it)

 while step < steps_end:
    solver.solve(nonLinearIterate=True, nonLinearMaxIterations = 20)
    dt=advDiff.get_max_dt()

    # increment strain. this is your numerical integration.
    swarm_strain.data[:] += strainRate.evaluate(swarm)*(strain_rate_scaling*dt*time_scaling)
    # project onto meshvariable
    strain_proj.solve()

    # Advect using this timestep size.
    advDiff.integrate(dt)
    # advect particles too. note we set the update_owners flag to false, 
    # because some particles end up outside the mesh (which hasn't been updated yet)
    # and we don't want them deleted.  the mesh will call update_owners when `deform_mesh`
    # is called.
    advector.integrate(dt, update_owners=False)  
    with mesh.deform_mesh():
        count = 0
        for i in range (0, res):
            for j in range (0, 4*res):
                index = 4 * res * (i + 1) - 1
                if temperatureField.data[count] <= 1.:
                    v = V_ff
                else:
                    v = velocityField.data[index]
                mesh.data[count][0] = mesh.data[count][0] + dt * drivingV(mesh.data[count][0], v, 0., v * time + 4 * res)
                count += 1

[rbeucher]

OK Try maybe:

#driving V function

def drivingV(x, v, minX, maxX):
    m = (v[0] - vxLeftWall) / (maxX-minX) 
    return m * (x - minX) + vxLeftWall  `

[jmansour]

Yep @rbeucher has pointed out at least one issue:

if temperatureField.data[count] <= 1.:
    v = V_ff   # SCALAR
else:
    v = velocityField.data[index]  # VECTOR

so it might even be better to use

if temperatureField.data[count] <= 1.:
    v = V_ff   # SCALAR
else:
    v = velocityField.data[index][0]   # NOW JUST GETTING V_X

However, I'm a bit concerned that you are moving mesh nodes only if they meet a certain temperature criterion as you will very likely tangle up your mesh. What are you actually attempting to model with this type of deformation?

[elisafierro]

Hi,

I tried to modify the drivingV function and also the variable 'v' to get V_X but I got a similar error (that I wrote down here).

What I'm attempting to model is the motion of the right plate relative to the left one (that is fixed) assuming variable velocity BCs at the right boundary of the model: constant velocity (V_ff = 1.5 / velocity_scaling_cmyr) from the surface (T= 0) to the bottom of the Thermal Boundary Layer (T=1), decreasing at the sub-lithosphere mantle as v += delta (where delta=V_ff / max_i). I also assume free slip BCs at the top and left boundaries of the model and no slip BCs on the bottom of the model.

Thanks for your help!! Elisa

/root/underworld2/underworld/systems/_bsscr.py:449: UserWarning: A petsc error has been detected during the solve.
The resultant solution fields are most likely erroneous, check them thoroughly.
This is possibly due to many things, for example solar flares or insufficient boundary conditions.
The resultant KSPConvergedReasons are (f_hat, outer, backsolve) (0,-3,0).
  warnings.warn(estring)

TypeErrorTraceback (most recent call last)
<ipython-input-29-6c68e1b4b7d6> in <module>()
    123 
    124                 mesh.data[count][0] = mesh.data[count][0] + dt * drivingV(mesh.data[count][0], 0., v * time + 4 * res)
--> 125                 count += 1
    126 
    127 

/usr/lib/python2.7/contextlib.pyc in __exit__(self, type, value, traceback)
     33                 value = type()
     34             try:
---> 35                 self.gen.throw(type, value, traceback)
     36                 raise RuntimeError("generator didn't stop after throw()")
     37             except StopIteration, exc:

/root/underworld2/underworld/mesh/_mesh.pyc in deform_mesh(self, isRegular, remainsRegular)
    254                                                      +"Note that any submesh should not be modified directly, "
    255                                                      +"but will instead be updated automatically on return from "
--> 256                                                      +"the 'deform_mesh' context manager. \nEncountered exception message:\n")
    257         finally:
    258             self.data.flags.writeable = False

/root/underworld2/underworld/mesh/_mesh.pyc in deform_mesh(self, isRegular, remainsRegular)
    247         self.data.flags.writeable = True
    248         try:
--> 249             yield
    250         except Exception as e:
    251             raise uw._prepend_message_to_exception(e, "An exception was raised during mesh deformation. "

<ipython-input-29-6c68e1b4b7d6> in <module>()
    122                     v = velocityField.data[index][0]   # NOW JUST GETTING V_X
    123 
--> 124                 mesh.data[count][0] = mesh.data[count][0] + dt * drivingV(mesh.data[count][0], 0., v * time + 4 * res)
    125                 count += 1
    126 

<ipython-input-23-8e3823a073e7> in drivingV(x, minX, maxX)
      1 def drivingV(x, minX, maxX):
----> 2     m = (v[0] - vxLeftWall) / (maxX-minX) # normalizzazione della velocità
      3     return m * (x - minX) + vxLeftWall # velocità normalizzata * posizione + velocità

TypeError: An exception was raised during mesh deformation. Your mesh may only be partially deformed. You can reset your mesh using the 'reset' method. Note that any submesh should not be modified directly, but will instead be updated automatically on return from the 'deform_mesh' context manager. 
Encountered exception message:

'float' object has no attribute '__getitem__'

[jmansour]

You need to ensure you pass drivingV(x, v, minX, maxX) a consistent v. Either always pass in v as a scalar (ie pass in V_x), or always pass in the vector (in which case use v[0] inside the function to extract V_x). Check my previous post again.

[elisafierro]

Hi,

as you suggested me, I tried to pass drivingV(x, v, minX, maxX):

1) a vector v (v[0]), but I got an error (I attached it here as a .png file called 'error')

2) then I tried to pass in v as a scalar (v=velocityField.data[index][0]) and:

   - After 1 hour I got outputs (but I just set 3 timesteps!!) 
   - The outputs are weird and not correct
   - I got a UserWarning (I attached it here as a .png file called 'UserWarning') after the model 
     finished to run and also some notes in the container logs (I attached it here as a .png file called 
     'Container_logs')

Then, I realized that I set 2 velocityBCs at the right boundary of the model, not 1:

 1. v_x=0, v_y=0 inside velocityBC instructions 

 2. a constant velocity (V_ff = 1.5 / velocity_scaling_cmyr) from the surface to the bottom of the Thermal Boundary Layer, decreasing at the sub-lithosphere mantle as v += delta:

Could that be the problem?

Thanks so much!!!!!

# VelocityBCs

vxWalls = mesh.specialSets["MinI_VertexSet"] + \
          mesh.specialSets["MinJ_VertexSet"] + \
          mesh.specialSets["MaxI_VertexSet"] 

vyWalls = mesh.specialSets["MinJ_VertexSet"] + \
          mesh.specialSets["MaxJ_VertexSet"] #+ \
          mesh.specialSets["MaxI_VertexSet"]

max_i = 0

V_ff = 1.5 / velocity_scaling_cmyr
vxLeftWall  = 0.
vyTopWall = 0.
for i in range (0, res):
    index = 4 * res * (i + 1) - 1
    if temperatureField.data[index] < 1.:
        max_i = i - 1
        break
delta = V_ff / max_i
v = 0.

for index in mesh.specialSets["MinI_VertexSet"]:
    if temperatureField.data[index] < 1.:
        velocityField.data[index] = [0.,0.]
    else:
        velocityField.data[index,0] = vxLeftWall

for index in mesh.specialSets["MinJ_VertexSet"]:
    velocityField.data[index] = [0.,0.]
for index in mesh.specialSets["MaxJ_VertexSet"]:
    velocityField.data[index,1] = vyTopWall
for index in mesh.specialSets["MaxI_VertexSet"]:
    if temperatureField.data[index] < 1.:
        velocityField.data[index][0] = V_ff
    else:
        velocityField.data[index][0] = v
        v += delta

velocityBC = uw.conditions.DirichletCondition( variable = velocityField, 
                                               indexSetsPerDof = (vxWalls, vyWalls) )`

[jmansour]

Have you visualised the mesh @elisafierro? You can use the glucifer.objects.Mesh() drawing object.

I highly suspect your mesh deformations are problematic. I would highly suggest you move the mesh uniformly first (pure extension), and only after you have this working satisfactorily move to more complex deformations.

Also, you will need to be a bit careful with how you advect particles when you are moving the mesh. Generally you should advect particles first, setting the update_owners=False. This prevents issues with the particles moving outside of the mesh (because they are advected first). Check the 2_09_ShearBandsPureShear example for details.

[jmansour]

Closing due to inactivity. Reopen if necessary.