febio_model.py

###################################################################################
#   FEBio-UncertainSCI
#=======================
# 
# This code applies the UncertainSCI library to FEBio. 
# The code allows users to evaluate uncertainty and sensitivity of FEBio model
# parameters on model output variables. 
#
# This code is distributed under the MIT license. See LICENSE file for details.
# Copyright 2021 - All rights reserved.
###################################################################################
import subprocess, time, platform
import numpy as np
#------------------------------------------------------------------------------
# This function creates the "run" file, which will run an FEBio model with
# the given input parameters and will write the values of the output parameters
# to an output file.
#
# PARAMETERS:
#   runFileName : the name of the FEBio run file. 
#   outFileName : the name of the output file that FEBio will write results to
#   p           : the values of the input parameters to use in this run
#   inparams    : the names of the input parameters
#   outparmas   : the names of the requested output parameters
#
def createRunFile(runFileName, outFileName, p, inparams, outparams):
    fp = open(runFileName, 'w')
    fp.write('<?xml version="1.0"?>\n')
    fp.write('<febio_run version="1.0">\n')
    fp.write('    <Parameters>\n')
    for i in range(p.shape[0]):
        fp.write('          <param name="');fp.write(inparams[i]);fp.write('">'); fp.write(str(p[i])); fp.write('</param>\n')
    fp.write('    </Parameters>\n')
    fp.write('    <Output>\n')
    fp.write('      <file>');fp.write(outFileName);fp.write('</file>\n')
    for p in outparams:
        fp.write(f'      <param name="{p}"/>\n')
    fp.write('    </Output>\n')
    fp.write('</febio_run>\n')
    fp.close()

#------------------------------------------------------------------------------
# This function processes the FEBio output file that contains the values of the 
# requested output parameters and returns the values from the file. 
#
# PARAMETERS:
#   outFileName : the name of the FEBio output file name. 
def get_febio_output(outFileName):
    of = open(outFileName, 'r')
    lines = of.readlines()
    v = np.empty([len(lines)])
    for i in range(len(lines)):
        v[i] = float(lines[i])
    return v

#------------------------------------------------------------------------------
# This function calls FEBio with the given values for the parameters p
#
# PARAMETERS:
#   p           : the values of the parameters at current sample point
#   febioFile   : name of FEBio input file. 
#   runFileName : name of the run file
#   outFileName : name of output file that stores values of output parameters
#   inparams    : list of names of all input parameters
#   outparams   : list of names of all output parameters
def febio_function(p, febioFile, runFileName, outFileName, inparams, outparams):

    # prepare the parameter run file
    createRunFile(runFileName, outFileName, p, inparams, outparams)

    # run febio
    print('calling febio ...', str(p))
    subprocess.run(['febio4', '-i', febioFile, '-task=param_run', runFileName, '-silent'])

    # read the output file
    v = get_febio_output(outFileName)
    print(v, 'done\n')

    return v

#------------------------------------------------------------------------------
# write the model output to a file. This file will be read during a restart
#
# PARAMETERS:
#   model_output : output generated by FEBio
#   fileName     : the name of the file to write data to
#
def write_model_output(model_output, fileName):
    N = model_output.shape[0]
    M = model_output.shape[1]
    fp = open('pceresults.txt', 'w')
    for j in range(M):
        for i in range(N):
            if i!=0:
                fp.write(' ')
            fp.write(str(model_output[i, j, 0]))
        fp.write('\n')
    fp.close()

#------------------------------------------------------------------------------
# Generate the FEBio output from the provided samples
#
# PARAMETERS:
#   samples   : list of all samples to evaluate
#   febioFile : name of FEBio input file
#   inparams  : list of input parameter names
#   outparams : list of output parameter names
def febio_output(samples, febioFile, inparams, outparams):
    # get sample size and dimensions
    M = samples.shape[0]

    # nr of output variables
    N = len(outparams)

    # create empty output array
    model_output = np.empty([N, M, 1])
    
    for ind in range(M):
        v = febio_function(samples[ind, :], febioFile, 'run.feb', 'out.txt', inparams, outparams)
        for i in range(N):
            model_output[i, ind, :] = v[i]

    # store all results to a file
    write_model_output(model_output, 'pceresults.txt')

    return model_output

#------------------------------------------------------------------------------
# This function reads the FEBio model output, that is, the file written by the 
# write_model_output function
#
# PARAMETERS:
#   outFileName : name of output file. 
#
def read_model_output(outFilename):
    of = open(outFilename, 'r')
    lines = of.readlines()
    of.close()
    M = len(lines)
    N = len(lines[0].split(' '))
    v = np.empty([N, M, 1])
    for j in range(M):
        d = lines[j].split(' ')
        for i in range(N):
            v[i, j, 0] = float(d[i])
    return v

#------------------------------------------------------------------------------
# Creates model output by calling FEBio in parallel on local computer.
#
# PARAMETERS:
#   samples          : the list of samples to evaluate
#   febioFile        : name of FEBio input file
#   inparams         : list of input parameter names
#   outparams        : list of output parameter names
#   numParallelJobs  : number of parallel FEBio jobs to execute
#   numThreadsPerJob : number of (OpenMP) threads to use for each job
#
def febio_output_parallel(samples, febioFile, inparams, outparams, numParallelJobs, numThreadsPerJob):
    
    # get sample size
    totalJobs = samples.shape[0]

    # nr of output variables
    N = len(outparams)

    # create empty output array
    model_output = np.empty([N, totalJobs, 1])
    
    # find febio3 full path for use in Popen
    febio3 = "febio3"
    if platform.system() != "Windows":
        out = subprocess.run(['which', 'febio3'], capture_output=True)
        febio3 = out.stdout.decode('utf8').strip()

    # This will store the processes
    jobs = {}

    # start all jobs
    current = 0
    while True:
        while (current < totalJobs) and (len(jobs) < numParallelJobs):

            currentString = str(current)
            runFileName = 'run' + currentString + '.feb'
            outFileName = 'out' + currentString + '.txt'

            # create the run file
            createRunFile(runFileName, outFileName, samples[current, :], inparams, outparams)

            # start FEBio
            command = [febio3, '-i', febioFile, '-task=param_run', runFileName, '-silent']
            jobs[currentString] = subprocess.Popen(command, env={"OMP_NUM_THREADS": str(numThreadsPerJob)}, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
            current += 1

        # if there are no more jobs running, we're done
        if (len(jobs) == 0):
            break

        # find all finished jobs
        finished = []
        for jobId in jobs:
            if jobs[jobId].poll() != None:
                finished.append(jobId)
        
        # process all finished jobs
        for jobId in finished:
            outFileName = 'out' + jobId + '.txt'
            v = get_febio_output(outFileName)

            id = int(jobId)
            for i in range(N):
                model_output[i, id, :] = v[i]

            print(jobId, 'done')

            # cleanup
            del jobs[jobId]

        # take a little nap
        time.sleep(0.001)

    # write model output to file
    write_model_output(model_output, 'pceresults.txt')

    # All done, so return 
    return model_output