Hi Romain,

We spoke about this last week, but restarting a model with a pressure BC has very odd effects. They seem to be independant of the model state too (restarting a model from timestep = 1 or timestep = 100 gives a similar outcome).

I have a simple model here: https://gist.githubusercontent.com/LukeMondy/8a5c0cb0d4953de1e05a6ea25e2d5085/raw/7072e5dbd93553505a88696c5fe138d67e24ef3f/stressBC_issue.py

You can see in this video a model that starts with a stress BC, runs for 5e6 years, is then restarted, and runs for another 5e6 years:
https://giphy.com/gifs/2mCSvMdixwUiE9Nb2x/html5
The left hand side shows the swarm, coloured by material index, and the right hand side shows the total pressure along the bottom wall.

The most obvious thing that is weird is that after the restart, the deformations changes a lot, and gets into a "lumpy" steady-state. The cause of this (I think) is that the pressure at the bottom of the model before the restart fluctuates, presumably in response to some dynamic pressures - but after the restart it stays fixed.

Any ideas?

(side note: I know my way of saving and reloading the pressure is weird, but this happens with or without that)

Hi Romain,

I am using UWGeodynamics on Magnus. I would like to use the inner solver 'superlusidt', however it appears that petsc was configured without the package superlu_dist and my model can't run. Could you please check if this is the case?

Thank you for your help!
Roberta

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for possible LU and Cholesky solvers
[0]PETSC ERROR: Could not locate solver package superlu_dist. Perhaps you must ./configure with --download-superlu_dist
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.8.3, Dec, 09, 2017
[0]PETSC ERROR: test_mymodel3d.py on a arch-linux2-c-opt named nid00218 by rcarluccio Wed Jun 20 11:32:54 2018
[0]PETSC ERROR: Configure options --download-hdf5=1 --with-debugging=0 --prefix=/opt/petsc --download-fblaslapack=1 --COPTFLAGS="-g -O2" --CXXOPTFLAGS="-g -O2" --FOPTFLAGS="-g -O2" --download-petsc4py=1 --download-mumps=1 --download-parmetis=1 --download-metis=1 --download-superlu=1 --download-hypre=1 --download-cmake=1 --download-scalapack=1
[0]PETSC ERROR: #1 MatGetFactor() line 4344 in /tmp/petsc-3.8.3/src/mat/interface/matrix.c
[0]PETSC ERROR: #2 PCSetUp_LU() line 93 in /tmp/petsc-3.8.3/src/ksp/pc/impls/factor/lu/lu.c
[0]PETSC ERROR: #3 PCSetUp() line 924 in /tmp/petsc-3.8.3/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #4 KSPSetUp() line 381 in /tmp/petsc-3.8.3/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #5 KSPSolve() line 612 in /tmp/petsc-3.8.3/src/ksp/ksp/interface/itfunc.c
srun: error: nid00218: tasks 1-7: Aborted
srun: Terminating job step 2935332.0
srun: Job step aborted: Waiting up to 32 seconds for job step to finish.
slurmstepd: error: *** STEP 2935332.0 ON nid00218 CANCELLED AT 2018-06-20T11:33:31 ***
srun: error: nid00218: task 0: Killed

The title is self explanatory, and enhancement for future version and for advanced use cases and for user who are at the Zen of PetSc

solver.options.scr
solver.options.scr.ksp_type

Hey,

The migration to python3 has found some small issues:

diff --git a/UWGeodynamics/Underworld_extended/_swarm.py b/UWGeodynamics/Underworld_extended/_swarm.py
index 3057377..3e76b69 100644
--- a/UWGeodynamics/Underworld_extended/_swarm.py
+++ b/UWGeodynamics/Underworld_extended/_swarm.py
@@ -185,7 +185,7 @@ class Swarm(uw.swarm.Swarm):
         if try_optimise:
             procCount = h5f.attrs.get('proc_offset')
             if procCount is not None and nProcs == len(procCount):
-                for p_i in xrange(rank):
+                for p_i in range(rank):
                     offset += procCount[p_i]
                 size = procCount[rank]
 
@@ -193,7 +193,7 @@ class Swarm(uw.swarm.Swarm):
         chunk=int(1e4) # read in this many points at a time
 
         (multiples, remainder) = divmod( size, chunk )
-        for ii in xrange(multiples+1):
+        for ii in range(multiples+1):
             # setup the points to begin and end reading in
             chunkStart = offset + ii*chunk
             if ii == multiples:
/Underworld_extended/_swarmvariable.py
index 3f6d196..58a0c2b 100644
--- a/UWGeodynamics/Underworld_extended/_swarmvariable.py
+++ b/UWGeodynamics/Underworld_extended/_swarmvariable.py
@@ -92,7 +92,7 @@ class SwarmVariable(uw.swarm.SwarmVariable):
 
         # calculate the hdf5 file offset
         offset=0
-        for i in xrange(rank):
+        for i in range(rank):
             offset += procCount[i]
 
         # open parallel hdf5 file

Python3 doesn't have xrange any more.

Hello Romain,

I plan to add a thermal anomaly defined by a rectangle in the asthenosphere. This rectangle is another material defined by a polygon shape.

I was wondering when we set the temperature boundary conditions, and reach the steady state temperature solution, does this temperature boundary conditions continue to be active during model run? Can I define temperature boundary conditions for the "model run" and initial temperature for the rectangle in my asthenosphere having 100 Kelvin higher than the background.

PoiseuilleFlowModel.set_temperatureBCs(top=293.15 * u.kelvin, 
                         bottom=1626.15 * u.kelvin,
                         materials=[(asthenosphere, 1626.15 * u.kelvin),(air, 293.15 * u.kelvin),(anomaly,1626.15 * u.kelvin)])

PoiseuilleFlowModel.init_model()

PoiseuilleFlowModel.temperature.data[PoiseuilleFlowModel.projMaterialField.data[:]==5]=GEO.nd(1726.0 * u.kelvin)

I was doing this routine in UW1 by calling the steady-state temperatures (TemperatureField.###.h5) in lmrInitials.xml and defining a hotter region on top that steady-state solution to model the thermal anomaly. At the same time, i could define other temperature BC as i like.

On magnus, I used shifter the pull the latest underworld_geodynamics:magnus, and tried to start a job.

It reported an issue in initialising the UWGeo module, saying http.client not found.

I will attach the error in a little bit - can't access Magnus where I am - but I think any model would trigger this issue.

So this line here:

Commenting it out makes things work fine.

It's a good addition though, I like what it does! Not sure why it causes issues.

Hello Romain,

when i was trying to track the velocities of each component (Vx and Vy) on the surface tracers using the following code:

x = np.linspace(GEO.nd(PoiseuilleFlowModel.minCoord[0]), GEO.nd(PoiseuilleFlowModel.maxCoord[0]), 10000)
y = GEO.nd(0. * u.kilometer)

surface_tracers = PoiseuilleFlowModel.add_passive_tracers(name="SurfaceTracer", vertices=[x, y])
surface_tracers.add_tracked_field(PoiseuilleFlowModel.velocityField, name="Velocities" , dataType="float", units=(u.centimeter / u.year), count=1, overwrite=True)

I got the following error:

File "/opt/UWGeodynamics/UWGeodynamics/_utils.py", line 157, in save
    obj.data[...] = field["value"].evaluate(self.swarm)
ValueError: could not broadcast input array from shape (1250,2) into shape (1250,1)
    tracers.save(self.outputDir, checkpointID, self.time)

So, how can we track vector fields on passive tracers?

svgwrite looks promising? ideas?

Hi Romain,

In https://uwgeodynamics.readthedocs.io/en/latest/User%20Guide.html#the-material-object/The Material object/CombinedShape/. you should change:
shape = GEO.shapes.CombinedShape([disk1, disk2]) FOR shape = disk1 & disk2

In the HalfSpace section below, to define the 3D shape (figure above), after the halfspaces definition:

halfspace1 = GEO.shapes.HalfSpace(normal=(-1.,0.,1.), origin=(4000. * u.km, 0. * u.km, -1000. * u.km))
halfspace2 = GEO.shapes.HalfSpace(normal=(0.,0.,1.), origin=(7000. * u.km, 1000. * u.km, 0. * u.km))
halfspace3 = GEO.shapes.HalfSpace(normal=(1.,0.,0.), origin=(9000. * u.km, 1000. * u.km, -500. * u.km))
halfspace4 = GEO.shapes.HalfSpace(normal=(0.,0.,-1.), origin=(6500. * u.km, 1000. * u.km, -1000. * u.km))

You should add:

compositeShape = halfspace1 & halfspace2 & halfspace3 & halfspace4
polygon= Model.add_material(name="polygon", shape=CompositeShape)

Cheers

Claire

the stdout message in the custom nonlinear solver branch has an error message, preventing the function to run,

sys.__stdout__.write(
   1408                         """Non linear solver - Residual {2:.8e}; Tolerance {3:.8e}""".format(
-> 1409                             residual, self._curTolerance))
   1410 
   1411             converged = False

IndexError: tuple index out of range

Hello Romain,

When I do the following:

GEO.rcParams["nonlinear.tolerance"]=0.2

The model still takes the tolerance value as 0.01 which is the default value.

Why is that so?
Omer

Hello Romain,

I am running a model with thermal anomaly defined as 100 K higher than the ambient asthenosphere. The residual for the non-linear solver is sometimes 500, 1.3, 0.4, so much varies. The tolerance is 0.01 for the model, so convergence takes a bit of time. In fact, it did not converge yet after 100 iteration with 80 CPUs, 2km resolution using 1920x384 mesh, 1.0 cm/yr maximum wall velocity, and penalty value of 1e6, on Raijin.

My question is, can penalty value force the calculation to be rough and give high residuals in each calculation?

Hello Guys,

I wanted to define conditional bottom bc using the following:

ConditionsBottom=fn.branching.conditional([(PoiseuilleFlowModel.x < GEO.nd(-1800.0 * u.kilometer), GEO.nd(25000.0 * u.megapascal)),
                                           (True, GEO.nd(0.0 * u.centimeter / u.year))]) #

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

Interestingly, I got the following error:

IndexErrorTraceback (most recent call last)
<ipython-input-26-618612d17b79> in <module>()
     23                                     top=[0.0 * u.centimeter / u.year, None],
     24                                     right=[0.0 * u.centimeter / u.year, None],
---> 25                                     left=[0.0 * u.centimeter / u.year, None])

/opt/UWGeodynamics/UWGeodynamics/_model.py in set_velocityBCs(self, left, right, top, bottom, front, back, indexSets)
    739                                        bottom=bottom, front=front,
    740                                        back=back, indexSets=indexSets)
--> 741         return self.velocityBCs.get_conditions()
    742 
    743     set_mechanicalBCs = set_velocityBCs

/opt/UWGeodynamics/UWGeodynamics/_velocity_boundaries.py in get_conditions(self)
    298         for set_ in self.order_wall_conditions:
    299             (condition, nodes) = self._wall_indexSets[set_]
--> 300             self.apply_condition_nodes(condition, nodes)
    301 
    302         if self.indexSets:

/opt/UWGeodynamics/UWGeodynamics/_velocity_boundaries.py in apply_condition_nodes(self, condition, nodes)
    201                     self.Model.velocityField.data[nodes.data, dim] = (
    202                         func.evaluate(
--> 203                             self.Model.mesh.data[nodes.data])[:, dim])
    204                     self.Model.boundariesField.data[nodes.data, dim] = (
    205                         func.evaluate(

IndexError: index 1 is out of bounds for axis 1 with size 1

Any idea what is going on?

Hi Romain,

I've noticed that after a restart, all the sediments loose their time of deposition. I think it's because of this line:

Cheers,
Luke

Hi, everyone, I tested the simple 3D programm,

the code as following attached
simple-rift-3D-2019-0227.ipynb.tar.gz

I don't know what is means that screen showed massage

Running with UWGeodynamics version 2.7.1 Options: -remove_constant_pressure_null_space False -ksp_k2_type NULL -change_backsolve False -pc_type none -force_correction True -Q22_pc_type gkgdiag -change_A11rhspresolve False -ksp_type bsscr -rescale_equations False -restore_K False -A11_pc_type lu -A11_pc_factor_mat_solver_package mumps -A11_ksp_type preonly -scr_ksp_type fgmres -scr_ksp_rtol 1e-05 -A11_mg_active False

Hello,

Just experiencing a bit of trouble with restarting a model. So I run this model: https://github.com/LukeMondy/Continental_Rifting/blob/master/lmondy_rift.py

for 300,000 years, which produces 3 checkpoints.

When I re-run the model, it writes this out:

Running with UWGeodynamics version 1.0.1
================================================================================

Restarting Model from Step 3 at Time = 0.0 year

(2018-12-14 13:49:17)
================================================================================

which is printed by:

So the bit that is weird is that it knows that Step 3 is the right one, but the Time = 0.0 year is not right. As the model progresses, it essentially overwrites the earlier model:
Step: 1 Model Time: 10000.00 year dt: 10000.00 year (2018-12-14 13:50:40)

Any ideas?

Exposing certain _variables would help out a lot, especially w.r.t stresses, i.e yield stress function (evaluating it over some passive markers in creating a profile etc, for example), and providing methods for fetching/getting shear and normal stresses on the swarm and or the mesh.

Linear and or exponential healing of accumulated plastic strain plastic over a user assigned time period.
Should be relatively simple/trivial?, could try implementing myself

I tried to run a model with no passive tracers, and it failed.

I think this block needs to be indented to fix:

Hello,

Is it possible that I define specific traction field magnitudes for each mesh node number in the left and right side walls?

I tried something like this,

appliedTractionField=PoiseuilleFlowModel.mesh.add_variable(nodeDofCount=2)

appliedTractionField.data[:,0]=0.00

.. Then I assign pressures to appliedTractionField.data array with some values having dimensions (Pascal) corresponding to the mesh node numbers of the side walls. Then I non-dimensionalise the array for further use, and do the following:

PoiseuilleFlowModel.set_velocityBCs(left=[appliedTractionField.data[PoiseuilleFlowModel.mesh.specialSets["MinI_VertexSet"].data[:],0] * GEO.Dimensionalize(1.0, u.pascal) ,0.0 * GEO.Dimensionalize(1.0, u.pascal)],
   							     right=[appliedTractionField.data[PoiseuilleFlowModel.mesh.specialSets["MaxI_VertexSet"].data[:],0] * GEO.Dimensionalize(1.0, u.pascal) ,0.0 * GEO.Dimensionalize(1.0, u.pascal)],
   							     top=[0.0 * u.centimeter / u.year, 0.0 * u.centimeter / u.year],
   							     bottom=[None, 0.0 * u.centimeter / u.year])

but seems that I am doing it wrong..

So I'm going through updating my model to match the changes in 1.0.0, and I've come across an issue where using nodeSets in boundary conditions is no longer parallel safe.

For example, I had this:

air_mask = air_shape.fn.evaluate(Model.mesh.data)
air_nodes = Model.mesh.data_nodegId[air_mask].ravel()
astheno_mask = astheno_shape.fn.evaluate(Model.mesh.data)
astheno_nodes = Model.mesh.data_nodegId[astheno_mask].ravel()


# Temp initial conditions
Model.set_temperatureBCs(top=293.15 * u.degK, 
                         bottom=1603.15 * u.degK, 
                         nodeSets = [(air_nodes, 293.15 * u.degK), (astheno_nodes, 1603.15 * u.degK)])

On a single cpu, it's OK.
On more than 1, it says:

jovyan@7e4d8acb9a80:/host$ mpirun -np 2 python -u lmondy_rift.py 
loaded rc file /opt/UWGeodynamics/UWGeodynamics/uwgeo-data/uwgeodynamicsrc
	Global element size: 304x95
	Local offset of rank 0: 0x0
	Local range of rank 0: 152x95
In func WeightsCalculator_CalculateAll(): for swarm "I6XW1XSN__swarm"
	done 33% (4814 cells)...
	done 67% (9627 cells)...
	done 100% (14440 cells)...
WeightsCalculator_CalculateAll(): finished update of weights for swarm "I6XW1XSN__swarm"
Traceback (most recent call last):
  File "lmondy_rift.py", line 310, in <module>
Traceback (most recent call last):
  File "lmondy_rift.py", line 310, in <module>
    nodeSets = [(air_mask, 293.15 * u.degK), (astheno_mask, 1603.15 * u.degK)])
  File "/opt/UWGeodynamics/UWGeodynamics/_model.py", line 611, in set_temperatureBCs
    nodeSets = [(air_mask, 293.15 * u.degK), (astheno_mask, 1603.15 * u.degK)])
  File "/opt/UWGeodynamics/UWGeodynamics/_model.py", line 611, in set_temperatureBCs
An uncaught exception was encountered on processor 0.
An uncaught exception was encountered on processor 1.
    return self._temperatureBCs.get_conditions()
  File "/opt/UWGeodynamics/UWGeodynamics/_boundary_conditions.py", line 200, in get_conditions
    self._apply_conditions_nodes(condition, nodes)
  File "/opt/UWGeodynamics/UWGeodynamics/_boundary_conditions.py", line 163, in _apply_conditions_nodes
    return self._temperatureBCs.get_conditions()
  File "/opt/UWGeodynamics/UWGeodynamics/_boundary_conditions.py", line 200, in get_conditions
    self._apply_conditions_nodes(condition, nodes)
  File "/opt/UWGeodynamics/UWGeodynamics/_boundary_conditions.py", line 163, in _apply_conditions_nodes
    fromObject=nodes)
  File "/opt/underworld2/underworld/mesh/_mesh.py", line 1224, in __init__
    fromObject=nodes)
  File "/opt/underworld2/underworld/mesh/_mesh.py", line 1224, in __init__
    super(FeMesh_IndexSet,self).__init__(object,*args,**kwargs)
  File "/opt/underworld2/underworld/container/_indexset.py", line 599, in __init__
    super(FeMesh_IndexSet,self).__init__(object,*args,**kwargs)
  File "/opt/underworld2/underworld/container/_indexset.py", line 599, in __init__
    super(ObjectifiedIndexSet,self).__init__(*args,**kwargs)
  File "/opt/underworld2/underworld/container/_indexset.py", line 84, in __init__
    super(ObjectifiedIndexSet,self).__init__(*args,**kwargs)
  File "/opt/underworld2/underworld/container/_indexset.py", line 84, in __init__
    self.add(fromObject)
  File "/opt/underworld2/underworld/container/_indexset.py", line 132, in add
    self.add(fromObject)
  File "/opt/underworld2/underworld/container/_indexset.py", line 132, in add
    self._addremove(indices,True);
  File "/opt/underworld2/underworld/container/_indexset.py", line 177, in _addremove
    self._addremove(indices,True);
  File "/opt/underworld2/underworld/container/_indexset.py", line 177, in _addremove
    self._AddOrRemoveWithNumpyArray(indices, isadding)
    self._AddOrRemoveWithNumpyArray(indices, isadding)
  File "/opt/underworld2/underworld/container/_indexset.py", line 275, in _AddOrRemoveWithNumpyArray
  File "/opt/underworld2/underworld/container/_indexset.py", line 275, in _AddOrRemoveWithNumpyArray
    raise TypeError("Incompatible array type ({}) provided. Must be of integer type.".format(ndarray.dtype))
TypeError: Incompatible array type (bool) provided. Must be of integer type.
    raise TypeError("Incompatible array type ({}) provided. Must be of integer type.".format(ndarray.dtype))
TypeError: Incompatible array type (bool) provided. Must be of integer type.
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 1 in communicator MPI COMMUNICATOR 3 DUP FROM 0
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
^Cmpirun: abort is already in progress...hit ctrl-c again to forcibly terminate

I then tried using the newer method, as shown in: https://github.com/underworldcode/UWGeodynamics/blob/master/examples/1_25_Hot_Canon_Ball.ipynb
where you just pass a shape to the nodeSets. This also didn't work.
I then ran the Hot Canon Ball example in parallel, and see the same issue. (I converted it to a python script, and then ran python -u canon.py, and mpirun -np 2 python -u canon.py. This is with UWGeo 1.0.1)

The original method used to be parallel safe, so I'm not sure what changed so that it is not.

Cheers,
Luke

Surface Swarms could track the inflow from boundary, for use cases in compression, and or convergence

python 1_08_ViscoElasticHalfSpace.ipynb.py
loaded rc file /Users/arijit/codes/UWGeodynamics/UWGeodynamics/uwgeo-data/uwgeodynamicsrc
Maxwell relaxation time =  316.887646408 year 10000000000.0 second
Observation time        =  119999.901227 year 3.786828e+12 second
effective viscosity     =  9.97366222542e+20 pascal * second
        Global element size: 128x32
        Local offset of rank 0: 0x0
        Local range of rank 0: 128x32
[[ 0.    0.05]]
Running with UWGeodynamics version 0.8.0-dev-120d091(HEAD)
Traceback (most recent call last):
  File "1_08_ViscoElasticHalfSpace.ipynb.py", line 244, in <module>
    Model.run_for(600000. * u.years)
  File "/Users/arijit/codes/UWGeodynamics/UWGeodynamics/_model.py", line 1546, in run_for
    self.restart(step=restartStep, restartDir=restartDir)
  File "/Users/arijit/codes/UWGeodynamics/UWGeodynamics/_model.py", line 379, in restart
    step = indices[step]
IndexError: list index out of range
-------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code.. Per user-direction, the job has been aborted.
-------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[41761,1],0]
  Exit code:    1
--------------------------------------------------------------------------

due to setting default restartStep=-1

workaround
Model.run_for(1. * u.years,restartStep=None)
Do want to modify all notebooks and examples?

When I try to run my UW2 model with the development docker image, it now gets stuck after doing the energy solve.

jovyan@510e3d184d4a:/host$ mpirun -np 6 python lmondy_rift.py l
	Global element size: 200x90
	Local offset of rank 0: 0x0
	Local range of rank 0: 67x45
In SystemLinearEquations_NonLinearExecute

Non linear solver - iteration 0
Linear solver (SRTKKVC5__system-execute) 
Linear solver (SRTKKVC5__system-execute), solution time 1.150250e-02 (secs)
Non linear solver - iteration 1
Linear solver (SRTKKVC5__system-execute) 
Linear solver (SRTKKVC5__system-execute), solution time 1.057148e-02 (secs)
In func SystemLinearEquations_NonLinearExecute: Iteration 1 of 500 - Residual 0.0015126 - Tolerance = 0.01
Non linear solver - Residual 1.51264297e-03; Tolerance 1.0000e-02 - Converged - 7.001925e-02 (secs)

In func SystemLinearEquations_NonLinearExecute: Converged after 1 iterations.

If you ctrl-c it, it says:

^C[mpiexec@510e3d184d4a] Sending Ctrl-C to processes as requested
[mpiexec@510e3d184d4a] Press Ctrl-C again to force abort


=====================================================================================
Error running Underworld - Signal 2 'SIGINT' (Termination Request).
Isn't it wonderbubble to have CTRL-C?


=====================================================================================
Error running Underworld - Signal 2 'SIGINT' (Termination Request).
Isn't it wonderbubble to have CTRL-C?


=====================================================================================
Error running Underworld - Signal 2 'SIGINT' (Termination Request).
Isn't it wonderbubble to have CTRL-C?


=====================================================================================
Error running Underworld - Signal 2 'SIGINT' (Termination Request).
Isn't it wonderbubble to have CTRL-C?


=====================================================================================
Error running Underworld - Signal 2 'SIGINT' (Termination Request).
Isn't it wonderbubble to have CTRL-C?


=====================================================================================
Error running Underworld - Signal 2 'SIGINT' (Termination Request).
Isn't it wonderbubble to have CTRL-C?
HDF5: infinite loop closing library
      D,G,A,S,T,F,F,AC,FD,P,PL,FD,P,FD,P,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD
HDF5: infinite loop closing library
      D,G,A,S,T,F,F,AC,FD,P,PL,FD,P,FD,P,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD
HDF5: infinite loop closing library
      D,G,A,S,T,F,F,AC,FD,P,PL,FD,P,FD,P,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD,FD
jovyan@510e3d184d4a:/host$

When I run without mpirun, it seems to work fine.

I'm sorry to ask a stupid question@@"

I download the uwgeodynamics image in kitematic

My host has ipyparallel function for jupyter already but don't have this kind of function when I used kitematic

Is it normal @@?

Hello Romain,

When I execute the pbs script you recently provided in 3D_rift, UWGeodynamics is not found in my directories.

import UWGeodynamics as GEO
    raise ImportError("Can not find Underworld, please check your installation")

When I check the pbs script, it says:

#UW_DIR=$HOME/programs/underworld2_master # master branch
UW_DIR=$HOME/programs/underworld2_development # Development branch

So, I don't have underworld2_development branch in my home directory. How can I install it?

By the way I already have UWGeodynamics in /home/562/ofb562/opt/ but maybe it is the old one..

So, could you please help me to update my directories and be able to run some UWGeodyanamics models on Raijin?
Thanks

used with the Thermomechanical Model
Traceback (most recent call last):
File "Tutorial_1_ThermoMechanical_Model.py", line 242, in
Model.save(filename='twol.json')
File "/media/student/Outputs/UWgeo/UWGeodynamics/UWGeodynamics/_model.py", line 1516, in save
json.dump(self, f, sort_keys=True, indent=4, cls=ObjectEncoder)
File "/usr/lib/python2.7/json/init.py", line 189, in dump
for chunk in iterable:
File "/usr/lib/python2.7/json/encoder.py", line 443, in _iterencode
for chunk in _iterencode(o, _current_indent_level):
File "/usr/lib/python2.7/json/encoder.py", line 434, in _iterencode
for chunk in _iterencode_dict(o, _current_indent_level):
File "/usr/lib/python2.7/json/encoder.py", line 408, in _iterencode_dict
for chunk in chunks:
File "/usr/lib/python2.7/json/encoder.py", line 332, in _iterencode_list
for chunk in chunks:
File "/usr/lib/python2.7/json/encoder.py", line 442, in _iterencode
o = _default(o)
File "/media/student/Outputs/UWgeo/UWGeodynamics/UWGeodynamics/json_encoder.py", line 8, in default
return obj.to_json()
File "/media/student/Outputs/UWgeo/UWGeodynamics/UWGeodynamics/_material.py", line 86, in to_json
val = self.dict[attribute]
KeyError: 'density'

Hello Romain,

In each time you make an update on the development version of the UWGeodynamics, does it mean that the Raijin is also updated? If not, should i update it by myself, maybe weekly ?

Cheers
Omer

Hi Romain,
It looks like the restart does not work. The model goes through the loading procedure but then stopped.
Cheers
Patrice

=========================================================================== Restarting Model from Step 93 at Time = 0.0 year

(2018-12-14 02:59:46)

Mesh loaded(2018-12-14 02:59:46)
Swarm loaded(2018-12-14 02:59:50)
pressureField loaded(2018-12-14 02:59:52)
velocityField loaded(2018-12-14 02:59:53)
materialField loaded(2018-12-14 02:59:56)
plasticStrain loaded(2018-12-14 02:59:58)
temperature loaded(2018-12-14 02:59:59)
FSE_Mantle loaded(2018-12-14 02:59:59)

I was just trying to reproduce the issue #65 , but now I face a different issue.

I take one of these models: https://github.com/EarthByte/UW2-tests-and-benchmarks/blob/master/isostasy/Isostasy%205%20-%20weak%20centre%20with%20sediments%20-%20PressureBC.ipynb and export it to a python file.

When I run it inside the latest underworld2_geodynamics (latest or dev), with this:
python isostasy.py
it works fine.

When I run it with:
mpirun -np 2 python isostasy.py # or any other number of processes
I get:

mpirun -np 2 python pt2.py 
	Global element size: 25x25
	Local offset of rank 0: 0x0
	Local range of rank 0: 25x13
Linear solver (J0OZ9Q5U__system-execute) 

BSSCR -- Block Stokes Schur Compliment Reduction Solver 
AUGMENTED LAGRANGIAN K2 METHOD - Penalty = 0.000000

  Setting schur_pc to "uw" 


SCR Solver Summary:

  RHS V Solve:            = 0.002512 secs / 1 its
  Pressure Solve:         = 0.1921 secs / 81 its
  Final V Solve:          = 0.002317 secs / 1 its

  Total BSSCR Linear solve time: 0.226781 seconds

Linear solver (J0OZ9Q5U__system-execute), solution time 2.283564e-01 (secs)
loaded rc file /opt/UWGeodynamics/UWGeodynamics/uwgeo-data/uwgeodynamicsrc
An uncaught exception was encountered on processor 0.
RuntimeError: Failed to execute the callback function, please check if it's valid

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "pt2.py", line 184, in <module>
    Model.solve()
  File "/opt/UWGeodynamics/UWGeodynamics/_model.py", line 1470, in solve
    nonLinearTolerance=self._curTolerance)
  File "/opt/underworld2/underworld/timing.py", line 323, in timed
    return routine(*args, **kwargs)
  File "/opt/underworld2/underworld/systems/_bsscr.py", line 451, in solve
    libUnderworld.StgFEM.SystemLinearEquations_UpdateSolutionOntoNodes(self._stokesSLE._cself, None)
SystemError: <built-in function SystemLinearEquations_UpdateSolutionOntoNodes> returned a result with an error set
application called MPI_Abort(comm=0x84000004, 1) - process 0

===================================================================================
=   BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
=   PID 130 RUNNING AT 7ff61ef05b11
=   EXIT CODE: 1
=   CLEANING UP REMAINING PROCESSES
=   YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================

If I change the model resolution to be higher, like 100x100, or 96x96, I see the same result.

Hi Romain,

Since starting with the new UWGeodynamics container, I have observed some notable differences between models that should be equivalent to those previously run with the underworld2_geodynamics container. The geometry that evolves is significantly different and subsidence is too large. I also notice that the vertical stress (calculated by integrating vertically density * gravity) deviates through time from the stress base boundary condition i have prescribed (from 3.61 GPa at the model start, which is equivalent to the basal boundary condition, to 3.33 GPa after 20 myr). This drift did not occur in models run with the old container. Is the implementation of the stress base boundary condition different in the new container? This behavior looks similar to behavior of the lecode isostasy boundary condition, which we have discussed previously.

Here is the boundary condition I have set:
Model.set_velocityBCs(left=[-0.55 * u.centimeter / u.year, 0.0 * u.centimeter / u.year],
right=[0.55 * u.centimeter / u.year, 0.0 * u.centimeter / u.year],
bottom=[0.0 * u.pascal, 3.607e9 * u.pascal])

Thanks,
Pete

Passive Tracers Grid: not working in parallel / nor restart

Hi, I wonder how to plot Variable field v.s depth fig, like temperature field, viscosity, yield stress...etc. v.s depth fig in GEO?

Figures explain the issue
Materials plus plastic strain(material 2 is transparent)

plastic strain

Colour scales

Just to enable mesh refinements, i guess swarm initialization could be deferred if a meshModification flag is passed on to the model initialization and called(e.g. Model.init_swarms() ) by user after he/she/they have made modification to the mesh.
Any Ideas?
if you favour this one i could send in a PR for review

Hello Romain,

I realized that when the convergence is not achieved after maximum number of iteration, the code does not stop, but tries again from 0. Why?

Hello Romain,

This afternoon, I had an issue that I could not fix. When I run the model using 8 cores that I was able to run, I could not go beyond some stage in calculation:

swig/python detected a memory leak of type 'Vec *', no destructor found.
swig/python detected a memory leak of type 'Vec *', no destructor found.
swig/python detected a memory leak of type 'Vec *', no destructor found.
In SystemLinearEquations_NonLinearExecute

Non linear solver - iteration 0
Linear solver (KEBIP2FO__system-execute) 

BSSCR -- Block Stokes Schur Compliment Reduction Solver 


-----  K2_GMG  ------

AUGMENTED LAGRANGIAN K2 METHOD - Penalty = 1000000.000000


	* K+p*K2 in time: 0.246727 seconds

  Setting schur_pc to "gkgdiag" 

Then, there is no response for hours... It was perfectly working before...

My model is:

PoiseuilleFlowModel = GEO.Model(minCoord=(-1920. * u.kilometer, -700. * u.kilometer),
                  maxCoord=(1920. * u.kilometer, 68. * u.kilometer),  
                  elementRes=(1920, 384),# 5 km resolution in both axes
                  gravity=(0.0, -9.8 * u.meter / u.second**2))

The boundary conditions are:

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

I use the version below:

Running with UWGeodynamics version 0.5.0-dev-64b24e0(development)

Most model run for a few time steps(about .5Ma) then fail with this

In func WeightsCalculator_CalculateAll(): for swarm "E3DSDPHZ__swarm" │··········
done 33% (10667 cells)... │··········
done 67% (21334 cells)... │··········
done 100% (32000 cells)... │··········
WeightsCalculator_CalculateAll(): finished update of weights for swarm "E3DSDPHZ__swarm" │··········
Processing surface with Badlands │··········
Traceback (most recent call last): │··········
File "two_weak_layers.py", line 122, in │··········
Model.run_for(8 * u.megayears, .01 * u.megayears) │··········
File "/home/student/ssd/uwgeo/UWGeodynamics/UWGeodynamics/_model.py", line 1321, in run_for │··········
self._update() │··········
File "/home/student/ssd/uwgeo/UWGeodynamics/UWGeodynamics/_model.py", line 1389, in _update │··········
self.surfaceProcesses.solve(dt) │··········
File "/home/student/ssd/uwgeo/UWGeodynamics/UWGeodynamics/surfaceProcesses.py", line 51, in solve │··········
self._BadlandsModel.solve(dt) │··········
File "/home/student/ssd/uwgeo/UWGeodynamics/UWGeodynamics/linkage/linkage.py", line 140, in solve │··········
self.badlands_model.run_to_time(self.time_years+dt_years) │··········
File "/home/student/pristine/pyBadlands_serial/pyBadlands/model.py", line 315, in run_to_time │··········
tflex=flexiso, scum=sload, Te=vTh,Ke=vKe, flexure=fflex, strat=fstrat, ero=fero) │··········
File "/home/student/pristine/pyBadlands_serial/pyBadlands/forcing/forceSim.py", line 752, in apply_XY_dispacements │··········
raise ValueError('Problem: IDs after merging is on previous vertex position.')

I had some basic ideas on this, derived from what is done in other codes.
creating getters or functions such as
getPlasticEta getViscousEta getElasticEta getMeltEta

followed by Different averaging schemes or selecting minimum of the three depending on a rcParams
something like effective.viscosity.scheme
@rbeucher do you agree, suggestions, or you ideas on the clean up?

left_InFlow = [(Model.y >= GEO.nd(-100.*u.kilometer), None),(True, 0.)]

Model.set_velocityBCs(left=[left_RidgeInFlow, None],
                      right=[0., None],
                      bottom=[None, 0.],
                      top=[None, 0.])

  File "/Users/arijit/codes/UWGeodynamics/UWGeodynamics/_model.py", line 837, in set_velocityBCs
    return self.velocityBCs.get_conditions()
  File "/Users/arijit/codes/UWGeodynamics/UWGeodynamics/_velocity_boundaries.py", line 312, in get_conditions
    self.apply_condition_nodes(condition, nodes)
  File "/Users/arijit/codes/UWGeodynamics/UWGeodynamics/_velocity_boundaries.py", line 224, in apply_condition_nodes
    func = fn.branching.conditional(condition[dim])
  File "/Users/arijit/underworld2/underworld/timing.py", line 320, in timed
    return routine(*args, **kwargs)
  File "/Users/arijit/underworld2/underworld/function/branching.py", line 225, in __init__
    self._fncself.insert( clause[0]._fncself, clause[1]._fncself )
AttributeError: 'NoneType' object has no attribute '_fncself'

ping me if you need an example, use case being keeping part of a wall free-slip and other part open.

Hello,

I was getting very unexpected behavior in my models, where at t=0 everything looked OK, but then at t=1, my entire pressure field would flip (so high pressure in the air, low in the mantle) - but my velocities were still pointing up (that being, material was flowing from low to high pressure)!

After much confusion, I discovered that by defining gravity as a vector, rather than a scalar, it has lost its sign somehow:

My fix was to define gravity in my models as such:

 Model = GEO.Model(elementRes=resolution,
                   minCoord=(-300 * u.kilometer, -180 * u.kilometer),
                   maxCoord=( 300 * u.kilometer,  20 * u.kilometer),
                   gravity = (0., -9.81 * u.m / u.s**2),                                                                                 
                   outputDir = output_dir)

Cheers,
Luke

after restarting a model, the swarm plots( defined in _plots.py ) don't work, as in they plot the initial model state, mesh plots do not show this behaviour or plots defined manually as glucifer figures in the script, affected ones are the material and viscosity fields. possible add a proper manual for restart???

Hello Romain,

The projStressField.h5 output has no units. How can I make it Pascal?
and, it would be cool if I had the 2nd invariant of the deviatoric stress tensor as another output :)
Cheers
Omer

Hello Romain,
I have been running some coupled UWGeo-Badlands pull-apart models. My issue is that the landscape doesn't move laterally as expected, interestingly there is no issue with the vertical motion of the surface. It seems that either the lateral motion of the UWGeo surface is not passed to Badlands, or it is passed with the wrong scale. The motion of the Badlands grid seems to confirm that there is a scaling issue.
Cheers
Patrice

Hi~
I'm a beginner @____@

I run the > Tutorial 3: The Numerical Sandbox, Extension Experiment (Mesh deform version)

with kitematic

I found the xdmf files produced from tutorial made paraview break down easily.
however they could open by visit

I thought if there were something wrong in the parameter "Model run" or not

I just guess, I'm looking for debugs and compare programs with underworld2

sorry to disturb

In many of the UW2 notebooks, they say something like this:

Note: the use of a pressure-sensitive rheology suggests that it is important to use a Q2/dQ1 element

(source:
https://github.com/underworldcode/underworld2/blob/debf718d3d3275ba3b6cdd8303b277a41463fe47/docs/development/models_broken/LemialeEtAl2008/data_comp/Lemiale-Parallel.py#L122 )

However, when you try to turn this on, the lithostatic pressure function fails, as it is unable to broadcast numpy arrays of different sizes.

In func SystemLinearEquations_NonLinearExecute: Converged after 1 iterations.
Traceback (most recent call last):
  File "lmondy_rift.py", line 142, in <module>
    Model.init_model()
  File "/opt/UWGeodynamics/UWGeodynamics/_model.py", line 1460, in init_model
    self.get_lithostatic_pressureField()
  File "/opt/UWGeodynamics/UWGeodynamics/_model.py", line 1361, in get_lithostatic_pressureField
    self.pressureField.data[:], LPressBot = lithoPress.solve()
ValueError: could not broadcast input array from shape (4500,1) into shape (18000,1)

Looking at the LithoPress code, I had a go at trying to interpolate the pressure result to the bigger dQ1 array, but got caught up in MPI stuff - but perhaps easier would be to project the density variable on to the subMesh (see here: https://github.com/rbeucher/UWGeodynamics/blob/master/UWGeodynamics/lithopress/lithopress.py#L34 ).

I had an attempt myself:

 lithoPress = LithostaticPressure(self.mesh.subMesh, self._densityFn, gravity)
 # from https://github.com/rbeucher/UWGeodynamics/blob/master/UWGeodynamics/_model.py#L1405

but there seems to be some strange differences between the self.mesh and the self.mesh.subMesh - mostly in that they don't have access to the same properties (maybe because the self.mesh is a UWGeodynamics.Underworld_extended._mesh.FeMesh_Cartesian object, while the subMesh is underworld.mesh._mesh.FeMesh)

This seems fairly important, as most of the stuff in UWGeodyamics does use pressure-sensitive rheologies, and in my models I'm seeing a lot of the checker-board pressure. Going to Q2/dQ0 made the model much more stable (in terms of convergence), but then fails because it's not an appropriate element for the _phi field.

Any help would be much appreciated!

worked before:
d8a2f61

Hello,

It would be very useful to be able to be able to pass options to the solver from UWGeodynamics.
Either getting access to the solver object like: Model.solver.options... or being able to set a method would be good:

def options(solver):
    solver.option.whatever = 'ok'
    return solver

Model.set_solver_callback(options)

Not quite sure on the details, but something like that would be great!

Hello Romain,

I learned that I can calculate it using the following:

import underworld.function as fn
fn.tensor.second_invariant(projStressField)

How can i make sure that I will have that output in all checkpoints?