Slope stability iterative factor of safety calculation

[g8star]

I’ve been using STKO for slope stability analyses, but I could use some help automating some of the tasks. In slope stability analyses with finite elements, the factor of safety is typically computed via the strength reduction method (e.g. Griffith and Fenton, 2004). I’ve been doing the strength reduction manually by setting material properties, running the analyses, and then checking the convergence. I then adjust the material properties and rerun the analyses. See the attached file for an example. The factor of safety should be between 1.2 and 1.3. By reducing the cohesion and friction angles of the PIMY materials, eventually the model does not converge.

Is there a way to programmatically modify the material properties during the analyses? Updating the properties after the elastic gravity stage would be ideal, potentially using something like updateParameter. Looping over the entire model building and analysis stages would be acceptable too, since each set of analyses runs pretty quickly.

strength reduction example.zip

(778.81 KiB) Downloaded 220 times

[g8star]

Michael Scott talks a little bit about updated OpenSees parameters within the analyses loop in his blog: https://portwooddigital.com/2021/08/15/ ... -the-loop/

All of the examples seem to be related to updating elastic material stiffness or UP element permeability. What about nonlinear material properties like cohesion and friction angle? Or other element properties like horizontal acceleration for pseudostatic analyses? Can these be updated using a similar method?

[STKO Team]

Is there a way to programmatically modify the material properties during the analyses?

Yes, but it depends on what you want to achieve.
You can use the parameter command if the material supports that parameter.
Otherwise you can use the PythonAPI in stko to programmatically change a material parameter and run a new analysis.

But first you should explain what you want to achieve. Updating a parameter during an analysis? or running a new analysis with the updated parameter?

[g8star]

I would like to update the cohesion ($cohesi) and friction angle ($frictionAng) of the PressureIndependMultiYield nDmaterial. Ideally, I would like to update the parameter during the analysis. As far as I can tell, the "parameter" command doesn't support updates of those parameters. So, we will have to programmatically change the parameters and then rerun the analysis. How do I do that with the PythonAPI? I use Python extensively, but I haven't used it in conjunction with STKO PythonAPI.

[STKO Team]

Please provide a small example with the element/materials you want, and a list of parameters you would like to change.
We'll show you how to setup a python script that runs those parametric analyses

[g8star]

Sure. I previously attached an example 2D slope stability problem, which I’ve attached again for convenience. The model consists of two layers, each defined by the PressureIndependMultiYield nDmaterial. I would like to programmatically modify the parameters:
• Sand fill: $frictionAng,
• Recent alluvium: $cohesi, and
• Recent alluvium: $frictionAng
The objective is to run the analysis, determine if the model converges, reduce $cohesi and tan($frictionAng) by a factor of 0.01, and then run the analysis again. The iterative process would continue until $cohesion and $frictionAng are small enough that the model no longer converges. The reduction factor at the final iteration defines the factor of safety (FS = 1/reduction factor).

The analysis runs fairly quickly with the provided parameters. The factor of safety should be somewhere between 1.2 and 1.3.

Let me know if you have any additional questions.

strength reduction example.zip

(778.81 KiB) Downloaded 143 times

[g8star]

The Python API provided enough information for me to develop my own solution. I've pasted a code snippet of my solution, which is intended only to work with the example problem posted previously. Some of the material and analysis properties hard-coded into the script and need to be modified for other scenarios.

Code: Select all

### Compute a list of tolerance by iteratively reducing strength parameters
### The last instance of convergence defines the factor of safety

from PyMpc import *
import numpy
import multiprocessing, subprocess
import time

#define soil properties
c = numpy.array([0.,9.8])
phid = numpy.array([30.6, 20])
phi = (180./numpy.pi)*numpy.arcsin((3.*numpy.tan(phid*numpy.pi/180))/(2.*numpy.sqrt(3)+numpy.tan(phid*numpy.pi/180)))  #octahedtral
ns = len(c)

F = numpy.arange(1,1.5,0.05)
nf = len(F)

doc = App.caeDocument()

def runf(Fi,pi=0):
	#pi is subprocess folder
	for j in range(ns):
		doc.physicalProperties[j+1].XObject.attributes['cohesi'].quantityScalar.value = c[j]/Fi
		doc.physicalProperties[j+1].XObject.attributes['frictionAng'].real = (180./numpy.pi)*numpy.arctan(numpy.tan(phi[j]*numpy.pi/180.)/Fi)
	di = '.\FS_{:}'.format(pi)
	App.runCommand('WriteInputFile',di)
	p = subprocess.Popen(['C:/OpenSees-Solvers/opensees.bat', '.\main.tcl'],shell=False,cwd=di) #your opensees solver goes here
	
	return p
	
#run it once to get appropriate tolerance
doc.analysisSteps[8].XObject.attributes['tol/NormDispIncr'].value = 1e6  #set an initial high value
ps = runf(F[0],pi=0)
while (ps.poll() is None):
	time.sleep(1)
mon = numpy.loadtxt('FS_{:}\\STKO_monitor_statistics.stats'.format(0))
tolt = numpy.max(mon[:,4])
print('Allowable tolerance {:}'.format(10**(int(numpy.log10(tolt))+1)))
doc.analysisSteps[8].XObject.attributes['tol/NormDispIncr'].value = 10**(int(numpy.log10(tolt))+1)

#run it sequentially until non-convergence
for i in range(1,nf):
	ps = runf(F[i])
	while (ps.poll() is None):
		time.sleep(1)
	mon = numpy.loadtxt('FS_{:}\\STKO_monitor_statistics.stats'.format(0))
	if len(mon)<20:  #initiate break, it didnt go through all iterations
		print('FS = {:}'.format(F[i-1]))
		break

[STKO Team]

Sorry for the late reply. Anyway, you did a great job.

There are many things you can do to further improve your script:

• using the subprocess Popen you can also think about running some of those analyses in parallel (for example 8 at a time if you have 8 processors). Have a look as multiprocessing queue in Python
• Use PySide2 widgets to let the user interact with settings
• much more...

[g8star]

Thanks for the additional information and ideas. I'm actually running it in parallel using subprocess.Popen in my implementation. I thought it would be simpler to just provide the code for sequential operation. Interestingly, multiprocessing.Pool opens multiple instances of the STKO GUI, which is not ideal.

[STKO Team]

Yes, if you use a high-level interface to multithreading or multiprocessing you may have those issues. Those libraries think you are running your code into a regular Python interpreter, while here the interpreter is inside STKO!

But you can do something like this:

Code: Select all

import queue
import subprocess
import threading
import random

# maximum number of concurrent jobs
max_processes = 4

# total number of jobs
num_jobs = 10

# the code to run in each process.
# to make each process have a different random duration
# the number in range() is a parameter
code_template = '''
import time
for i in range({}):
	print('iter',i)
	time.sleep(1)
'''

# the function that process the queue
def worker(queue):
	while not queue.empty():
		job_id, template_param = queue.get()
		try:
			# run the current job
			print(f'running job {job_id} on thread {threading.get_ident()} with {template_param} iterations')
			job = ['python', '-c', code_template.format(template_param)]
			process = subprocess.Popen(job)
			process.wait()
		finally:
			queue.task_done()

# create the queue
the_queue = queue.Queue()

# fill it with the jobs to run.
# randomly assign 3 to 10 iterations for each job
for job_id in range(num_jobs):
	the_queue.put((job_id, random.randint(3, 10)))

# create and start the working threads
all_threads = []
for _ in range(max_processes):
	the_thread = threading.Thread(target=worker, args=(the_queue,))
	the_thread.start()
	all_threads.append(the_thread)

# wait for all jobs to complete
the_queue.join()
for the_thread in all_threads:
	the_thread.join()

This code uses num_jobs jobs in max_processes processes using queues. basically it runs max_processes jobs (of the num_jobs jobs) at a time. If some of them finish before the others, other jobs will run, so that you always keep your processes busy

Note that this code will freeze STKO until everything is finished, so you may want to apply here some of the concepts I used in the webinar on the PythonAPI to avoid freezing the STKO's event loop.