Skip to content

Commit

Permalink
fixes for 1D0D Reymond validation example: turn off gravity, fix issu…
Browse files Browse the repository at this point in the history 
…es with time scale and reference areas, improve result analysis script
  • Loading branch information
@dladd
dladd committed Aug 18, 2015
1 parent 555452e commit aa0bb31
Show file tree
Hide file tree
Showing 7 changed files with 789 additions and 600 deletions.
162 changes: 101 additions & 61 deletions FluidMechanics/NavierStokes/Coupled1DCellML/Python/Reymond/Post.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,26 +2,32 @@
# Post Process #
##################

import os
import os,sys
import numpy as np
import subprocess
import pylab
from matplotlib import pyplot as plt
from subprocess import Popen, PIPE
import FluidExamples1DUtilities as Utilities

def Post(nodes):
# Set the reference values
Ts = 1.0#0.001 # Time (s)
Qs = 100.0e-8 # Flow (m3/s)
nodeCounter = 0

# Set the time parameters
solveTime = 1345.#790.0
timeIncrement = 0.05
cycleTime = 790.0
numberOfCycles = 4.0
solveTime = cycleTime*numberOfCycles
timeIncrement = 0.2
outputFrequency = 10

#Choose type of plot(s)
plotFlow = True
plotPressure = False

# Set up numpy arrays
numberOfTimesteps = (solveTime/timeIncrement)/outputFrequency
numberOfTimesteps = int((solveTime/timeIncrement)/outputFrequency)
timeResult = np.zeros((numberOfTimesteps))
numberOfNodes = len(nodes)
flowResult = np.zeros((numberOfNodes,numberOfTimesteps))
Expand All @@ -39,70 +45,104 @@ def Post(nodes):
Utilities.CsvNodeReader(filename,inputNodeNumbers,bifurcationNodeNumbers,trifurcationNodeNumbers,coupledNodeNumbers,nodeCoordinates,arteryLabels)

# Set up reference types
refs = ['Reymond2009','Reymond2011']
refTypes = ['Experimental','Model']
refs = ['Reymond2009'] # ,'Reymond2011'
refTypes = ['Model','Experimental']
colours = ['r','c','b','g','k']

# Node number for data extraction
nodeLabels = []
# TODO: should swap these loops around to speed things up for multiple nodes (currently re-reads files for each node)
for node in nodes:
# Get artery names
arteryLabel = arteryLabels[node-1]
print('Reading results for ' + arteryLabel)
nodeLabels.append(arteryLabel)
# Create the Result file to store the flow and area versus time for each node
createFile = open("results/node_"+str(node),'w+')
# Loop through the output files
timestep = 0
for x in range(0,int(solveTime/timeIncrement),outputFrequency):
outputFile = open("output/MainTime_"+str(x)+".part0.exnode")
outputLINE = outputFile.readlines()
for i in range(0,len(outputLINE)):
if outputLINE[i].split() == ['Node:', str(node)]:
# Extract the variables from output files
Flow = float(''.join(outputLINE[i+4].split()))
Pressure = float(''.join(outputLINE[i+12].split()))
Time = x*Ts*timeIncrement
Flow = Flow*Qs*1000000.0
# Write in the Result file
#Result = open("results/node_"+str(nodeCounter),'a+')
#Result = open("results/node_"+str(node),'a+')
#print >> Result,"%.4f"%Time,Flow,Pressure

# Store results
timeResult[timestep] = Time
flowResult[nodeCounter,timestep] = Flow
pressureResult[nodeCounter,timestep] = Pressure

timestep += 1


pylab.plot(timeResult,flowResult[nodeCounter,:],'r',alpha=0.5,label=arteryLabel)

for ref in refs:
for refType in refTypes:
filename = './reference/'+ref+'/'+refType+'/node_'+str(node)
if os.path.exists(filename):
refData = np.genfromtxt(filename,delimiter='\t')
refLabel = ref + ' ' + refType + ' ' + arteryLabel
print('Found reference '+refLabel)
pylab.plot(refData[:,0],refData[:,1],'--r',alpha=0.5,label=refLabel)

# Plot this node
pylab.show()



nodeCounter = nodeCounter+1

#Read in all node data
print("Reading output file data (this can take a while...)")
flowNodeData = np.zeros((totalNumberOfNodes,numberOfTimesteps))
pressureNodeData = np.zeros((totalNumberOfNodes,numberOfTimesteps))
#for timestep in range(0,int(solveTime/timeIncrement),outputFrequency):
for timestep in range(numberOfTimesteps):
filename = "output/MainTime_"+str(timestep*outputFrequency)+".part0.exnode"
#print(filename)

outputFile = open(filename)
lines = outputFile.readlines()
for i in range(0,len(lines)):
if 'Node:' in lines[i]:
nodeNumber = [int(s) for s in lines[i].split() if s.isdigit()][0]
#print(nodeNumber)
# Extract the variables from output files
Flow = float(lines[i+4].split()[0])
Flow = Flow*Qs*1000000.0
Pressure = float(lines[i+12].split()[0])

# Store results
flowNodeData[nodeNumber-1,timestep] = Flow
pressureNodeData[nodeNumber-1,timestep] = Pressure
outputFile.close()
timeResult = np.arange(0.0,solveTime,Ts*timeIncrement*outputFrequency)

if plotFlow:
nodeCounter = 0
for node in nodes:
# Get artery names
arteryLabel = arteryLabels[node-1]
print('Generating plot for ' + arteryLabel)
nodeLabels.append(arteryLabel)

# Do a subplot if looking at several nodes at once
if numberOfNodes > 1:
plt.subplot(int((numberOfNodes)/2)+numberOfNodes%2,2,nodeCounter+1)

# Plot this node
plt.plot(timeResult,flowNodeData[node-1,:],'k-',label='OpenCMISS Model')
plt.title(arteryLabel+' (Node '+str(node)+')')
colourNum = 0
for ref in refs:
for refType in refTypes:
filename = './reference/'+ref+'/'+refType+'/node_'+str(node)
if os.path.exists(filename):
refData = np.genfromtxt(filename,delimiter='\t')
numberOfRefTimesteps = len(refData[:,0])
refLabel = ref + ' ' + refType# + ' ' + arteryLabel
print('Found reference from '+refLabel)
refDataCycles=refData
addTime = np.zeros(numberOfRefTimesteps*numberOfCycles)
for cycle in range(1,int(numberOfCycles)):
refDataCycles = np.vstack((refDataCycles,refData))
addTime[numberOfRefTimesteps*cycle:]+=cycleTime
refDataCycles[:,0]=np.add(refDataCycles[:,0],addTime)
plt.plot(refDataCycles[:,0],refDataCycles[:,1],colours[colourNum],alpha=1.0,label=refLabel)
colourNum+=1

nodeCounter = nodeCounter+1
# Plot all nodes
plt.legend(loc = (1.2, 0.5))
plt.show()

if plotPressure:
nodeCounter = 0
for node in nodes:
# Get artery names
arteryLabel = arteryLabels[node-1]
print('Generating plot for ' + arteryLabel)
nodeLabels.append(arteryLabel)

# Do a subplot if looking at several nodes at once
if numberOfNodes > 1:
plt.subplot(int((numberOfNodes)/2)+numberOfNodes%2,2,nodeCounter+1)

# Plot this node
plt.plot(timeResult,pressureNodeData[node-1,:],'b-')
plt.title(arteryLabel+' (Node '+str(node)+')')
nodeCounter = nodeCounter+1
# Plot all nodes
plt.show()

return

nodes = input('Nodes: ')
if len(sys.argv) > 1:
nodes = int(sys.argv[1:])
else:
nodes = [1,93,120,127,170]
Post(nodes)

print "."
print "."
print "."
print "Processing Completed!"

#Popen(['gnuplot','PlotNode.p'])
Expand Down
52 changes: 29 additions & 23 deletions FluidMechanics/NavierStokes/Coupled1DCellML/Python/Reymond/ReymondExample.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,8 @@
# Initialise OpenCMISS and any other needed libraries
#================================================================================================================================
import numpy as np
import math,cmath,csv,time,sys,os,glob
import FluidExamples1DUtilities as Utilities
import math,csv,time,sys,os,glob,shutil
import FluidExamples1DUtilities as Utilities1D

sys.path.append(os.sep.join((os.environ['OPENCMISS_ROOT'],'cm','bindings','python')))
from opencmiss import CMISS
Expand Down Expand Up @@ -99,7 +99,7 @@
# Set the flags
RCRBoundaries = True # Set to use coupled 0D Windkessel models (from CellML) at model outlet boundaries
nonReflecting = False # Set to use non-reflecting outlet boundaries
CheckTimestepStability = True # Set to do a basic check of the stability of the hyperbolic problem based on the timestep size
CheckTimestepStability = False # Set to do a basic check of the stability of the hyperbolic problem based on the timestep size
initialiseFromFile = False # Set to initialise values
ProgressDiagnostics = True # Set to diagnostics

Expand All @@ -120,9 +120,9 @@
coupledNodeNumbers = []
arteryLabels = []
filename = 'input/Node.csv'
numberOfNodes = Utilities.GetNumberOfNodes(filename)
numberOfNodes = Utilities1D.GetNumberOfNodes(filename)
nodeCoordinates = np.zeros([numberOfNodes,4,3])
Utilities.CsvNodeReader(filename,inputNodeNumbers,bifurcationNodeNumbers,trifurcationNodeNumbers,coupledNodeNumbers,
Utilities1D.CsvNodeReader(filename,inputNodeNumbers,bifurcationNodeNumbers,trifurcationNodeNumbers,coupledNodeNumbers,
nodeCoordinates,arteryLabels)
numberOfInputNodes = len(inputNodeNumbers)
numberOfBifurcations = len(bifurcationNodeNumbers)
Expand All @@ -134,7 +134,7 @@
elementNodes.append([0,0,0])
bifurcationElements = (numberOfBifurcations+1)*[3*[0]]
trifurcationElements = (numberOfTrifurcations+1)*[4*[0]]
Utilities.CsvElementReader('input/Element.csv',elementNodes,bifurcationElements,trifurcationElements,numberOfBifurcations,numberOfTrifurcations)
Utilities1D.CsvElementReader('input/Element.csv',elementNodes,bifurcationElements,trifurcationElements,numberOfBifurcations,numberOfTrifurcations)
numberOfElements = len(elementNodes)-1

if (ProgressDiagnostics):
Expand All @@ -153,11 +153,9 @@
# Set the material parameters
Rho = 1050.0 # Density (kg/m3)
Mu = 0.004 # Viscosity (Pa.s)
D = 100.0 # Diffusivity (m2/s)
G0 = 9.81 # Gravitational acceleration (m/s2)
G0 = 0.0 # Gravitational acceleration (m/s2)
Pext = 0.0 # External pressure (Pa)
Pv = 2667.0 # Venous pressure = 20.0 mmHg (Pa)
Alpha = 1.3 # Flow profile type
Alpha = 1.0 # Flow profile type

# Material parameter scaling factors
Ls = 1000.0 # Length (m -> mm)
Expand All @@ -170,16 +168,14 @@
Rhos = Ms/(Ls**3.0) # Density (kg/m3)
Mus = Ms/(Ls*Ts) # Viscosity (kg/(ms))
Ps = Ms/(Ls*Ts**2.0) # Pressure (kg/(ms2))
Ds = (Ls**2.0)/Ts # Diffusivity (m2/s)
Zs = Ps/Qs # Impedance (pa/(m3/s))
Gs = Ls/(Ts**2.0) # Acceleration (m/s2)

# Initialise the node-based parameters
A0 = np.zeros((numberOfNodes+1,4)) # Area (m2)
H = np.zeros((numberOfNodes+1,4)) # Thickness (m)
E = np.zeros((numberOfNodes+1,4)) # Elasticity (Pa)
# Read the MATERIAL csv file
Utilities.CsvMaterialReader('input/Material.csv',A0,E,H)
Utilities1D.CsvMaterialReader('input/Material.csv',A0,E,H)

# Apply scale factors
Rho = Rho*Rhos
Expand All @@ -188,7 +184,6 @@
A0 = A0*As
E = E*Es
H = H*Hs
D = D*Ds
G0 = G0*Gs

Q = np.zeros((numberOfNodes+1,4))
Expand Down Expand Up @@ -235,7 +230,7 @@
# Set the time parameters
numberOfPeriods = 4.0
timePeriod = 790.
timeIncrement = 0.05
timeIncrement = 0.2
startTime = 0.0
stopTime = numberOfPeriods*timePeriod
dynamicSolverNavierStokesTheta = [1.0]
Expand All @@ -257,15 +252,21 @@
RESTART_VALUE = 3000 # default: 30

# N-S/C coupling tolerance
couplingTolerance1D = 1.0E+6
couplingTolerance1D = 1.0E+10
# 1D-0D coupling tolerance
couplingTolerance1D0D = 0.001

# Navier-Stokes solver
EquationsSetSubtype = CMISS.EquationsSetSubtypes.COUPLED1D0D_NAVIER_STOKES
# Characteristic solver
EquationsSetCharacteristicSubtype = CMISS.EquationsSetSubtypes.CHARACTERISTIC
ProblemSubtype = CMISS.ProblemSubTypes.COUPLED1D0D_NAVIER_STOKES
if(RCRBoundaries):
EquationsSetSubtype = CMISS.EquationsSetSubtypes.COUPLED1D0D_NAVIER_STOKES
# Characteristic solver
EquationsSetCharacteristicSubtype = CMISS.EquationsSetSubtypes.CHARACTERISTIC
ProblemSubtype = CMISS.ProblemSubTypes.COUPLED1D0D_NAVIER_STOKES
else:
EquationsSetSubtype = CMISS.EquationsSetSubtypes.TRANSIENT1D_NAVIER_STOKES
# Characteristic solver
EquationsSetCharacteristicSubtype = CMISS.EquationsSetSubtypes.CHARACTERISTIC
ProblemSubtype = CMISS.ProblemSubTypes.TRANSIENT1D_NAVIER_STOKES

#================================================================================================================================
# Coordinate System
Expand Down Expand Up @@ -648,7 +649,8 @@
nodeIdx = coupledNodeNumbers[terminalIdx-1]
nodeDomain = Decomposition.NodeDomainGet(nodeIdx,meshComponentNumber)
if (nodeDomain == computationalNodeNumber):
# Incoming
# Incoming (parent) - outgoing component to come from 0D
versionIdx = 1
IndependentFieldNavierStokes.ParameterSetUpdateNodeDP(CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.VALUES,
versionIdx,1,nodeIdx,1,1.0)

Expand Down Expand Up @@ -745,6 +747,7 @@
nodeDomain = Decomposition.NodeDomainGet(nodeIdx,meshComponentNumber)
if (nodeDomain == computationalNodeNumber):
#print("Terminal node: " + str(nodeIdx))
versionIdx = 1
CellMLModelsField.ParameterSetUpdateNode(CMISS.FieldVariableTypes.U,CMISS.FieldParameterSetTypes.VALUES,
versionIdx,1,nodeIdx,1,CellMLModelIndex[terminalIdx])

Expand Down Expand Up @@ -1018,6 +1021,7 @@
SolverEquationsCharacteristic.BoundaryConditionsCreateStart(BoundaryConditionsCharacteristic)

# Area-outlet
versionIdx = 1
for terminalIdx in range (1,numberOfTerminalNodes+1):
nodeNumber = coupledNodeNumbers[terminalIdx-1]
nodeDomain = Decomposition.NodeDomainGet(nodeNumber,meshComponentNumber)
Expand All @@ -1041,13 +1045,15 @@
SolverEquationsNavierStokes.BoundaryConditionsCreateStart(BoundaryConditionsNavierStokes)

# Inlet (Flow)
versionIdx = 1
for inputIdx in range (1,numberOfInputNodes+1):
nodeNumber = inputNodeNumbers[inputIdx-1]
nodeDomain = Decomposition.NodeDomainGet(nodeNumber,meshComponentNumber)
if (nodeDomain == computationalNodeNumber):
BoundaryConditionsNavierStokes.SetNode(DependentFieldNavierStokes,CMISS.FieldVariableTypes.U,
versionIdx,1,nodeNumber,1,CMISS.BoundaryConditionsTypes.FIXED_FITTED,Q[inputIdx][0])
# Area-outlet
versionIdx = 1
for terminalIdx in range (1,numberOfTerminalNodes+1):
nodeNumber = coupledNodeNumbers[terminalIdx-1]
nodeDomain = Decomposition.NodeDomainGet(nodeNumber,meshComponentNumber)
Expand All @@ -1067,7 +1073,7 @@

if (CheckTimestepStability):
QMax = 430.0
maxTimestep = Utilities.GetMaxStableTimestep(elementNodes,QMax,nodeCoordinates,H,E,A0,Rho)
maxTimestep = Utilities1D.GetMaxStableTimestep(elementNodes,QMax,nodeCoordinates,H,E,A0,Rho)
if (timeIncrement > maxTimestep):
sys.exit('Timestep size '+str(timeIncrement)+' above maximum allowable size of '+str(maxTimestep)+'. Please reduce step size and re-run')

Expand All @@ -1087,7 +1093,7 @@

# Remove CellML tmp files
for filename in glob.glob("./tmp.cellml2code*"):
os.remove(filename)
shutil.rmtree(filename)

#================================================================================================================================
# Finish Program
Expand Down
Loading

0 comments on commit aa0bb31

Please sign in to comment.