P283370/NuMRI
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

asdf

Browse Source
This commit is contained in:
jeremias 2020-02-17 11:00:29 +01:00
parent 49e6de207b
commit 66b0e7a8e6
9 changed files with 354 additions and 381 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

100
codes/MRI.py
View File

@ -18,13 +18,9 @@ from mpi4py import MPI
from common import inout from common import inout
import time import time
if '2017' in dolfin.__version__:
#class MPI(MPI): PRINT_BARS = False
# comm_world = mpi_comm_world() parameters["std_out_all_processes"] = False
set_log_active(False)
else:
PRINT_BARS = False
parameters["std_out_all_processes"] = False
def Etcetera(mode): def Etcetera(mode):
if mode=='happy': if mode=='happy':
@ -50,31 +46,31 @@ def Etcetera(mode):
def dealiased(VENC,vx,vy,vz): def dealiased(VENC,vx,vy,vz):
vx2 = np.zeros(vx.shape) vx2 = np.zeros(vx.shape)
vy2 = np.zeros(vy.shape) vy2 = np.zeros(vy.shape)
vz2 = np.zeros(vz.shape) vz2 = np.zeros(vz.shape)
badx_pos = np.where(vx>VENC) badx_pos = np.where(vx > VENC)
badx_neg = np.where(vx<-VENC) badx_neg = np.where(vx < -VENC)
bady_pos = np.where(vy>VENC) bady_pos = np.where(vy > VENC)
bady_neg = np.where(vy<-VENC) bady_neg = np.where(vy < -VENC)
badz_pos = np.where(vz>VENC) badz_pos = np.where(vz > VENC)
badz_neg = np.where(vz<-VENC) badz_neg = np.where(vz < -VENC)
# More than VENC # More than VENC
vx2[badx_pos] = -(2*VENC - vx[badx_pos]) vx2[badx_pos] = -(2*VENC - vx[badx_pos])
vy2[bady_pos] = -(2*VENC - vy[bady_pos]) vy2[bady_pos] = -(2*VENC - vy[bady_pos])
vz2[badz_pos] = -(2*VENC - vz[badz_pos]) vz2[badz_pos] = -(2*VENC - vz[badz_pos])
# Less than VENC # Less than VENC
vx2[badx_neg] = 2*VENC + vx[badx_neg] vx2[badx_neg] = 2*VENC + vx[badx_neg]
vy2[bady_neg] = 2*VENC + vy[bady_neg] vy2[bady_neg] = 2*VENC + vy[bady_neg]
vz2[badz_neg] = 2*VENC + vz[badz_neg] vz2[badz_neg] = 2*VENC + vz[badz_neg]
# The rest # The rest
otherx = np.where(vx2==0) otherx = np.where(vx2 == 0)
othery = np.where(vy2==0) othery = np.where(vy2 == 0)
otherz = np.where(vz2==0) otherz = np.where(vz2 == 0)
vx2[otherx] = vx[otherx] vx2[otherx] = vx[otherx]
vy2[othery] = vy[othery] vy2[othery] = vy[othery]
vz2[otherz] = vz[otherz] vz2[otherz] = vz[otherz]
return [vx2,vy2,vz2] return [vx2,vy2,vz2]
@ -1326,35 +1322,35 @@ def SCANNER(options):
# Creating M # Creating M
for t in range(Nt): for t in range(Nt):
Mx[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqx[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqx[:,:,:,t])) Mx[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqx[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqx[:,:,:,t]))
My[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqy[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqy[:,:,:,t])) My[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqy[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqy[:,:,:,t]))
Mz[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqz[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqz[:,:,:,t])) Mz[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqz[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqz[:,:,:,t]))
# Organs # Organs
Mx[:,:,:,t] = Mx[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) Mx[:,:,:,t] = Mx[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 )
Mx[:,:,:,t] = Mx[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) Mx[:,:,:,t] = Mx[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 )
My[:,:,:,t] = My[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) My[:,:,:,t] = My[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 )
My[:,:,:,t] = My[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) My[:,:,:,t] = My[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 )
Mz[:,:,:,t] = Mz[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) Mz[:,:,:,t] = Mz[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 )
Mz[:,:,:,t] = Mz[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) Mz[:,:,:,t] = Mz[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 )
FxG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqx[:,:,:,t]/VENC FxG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqx[:,:,:,t]/VENC
FyG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqy[:,:,:,t]/VENC FyG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqy[:,:,:,t]/VENC
FzG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqz[:,:,:,t]/VENC FzG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqz[:,:,:,t]/VENC
Fx0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) Fx0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4)
Fy0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) Fy0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4)
Fz0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) Fz0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4)
Sx20[:,:,:,t] = Mx[:,:,:,t]*np.cos(Fx0[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(Fx0[:,:,:,t]) Sx20[:,:,:,t] = Mx[:,:,:,t]*np.cos(Fx0[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(Fx0[:,:,:,t])
Sy20[:,:,:,t] = My[:,:,:,t]*np.cos(Fy0[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(Fy0[:,:,:,t]) Sy20[:,:,:,t] = My[:,:,:,t]*np.cos(Fy0[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(Fy0[:,:,:,t])
Sz20[:,:,:,t] = Mz[:,:,:,t]*np.cos(Fz0[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(Fz0[:,:,:,t]) Sz20[:,:,:,t] = Mz[:,:,:,t]*np.cos(Fz0[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(Fz0[:,:,:,t])
Sx2G[:,:,:,t] = Mx[:,:,:,t]*np.cos(FxG[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(FxG[:,:,:,t]) Sx2G[:,:,:,t] = Mx[:,:,:,t]*np.cos(FxG[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(FxG[:,:,:,t])
Sy2G[:,:,:,t] = My[:,:,:,t]*np.cos(FyG[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(FyG[:,:,:,t]) Sy2G[:,:,:,t] = My[:,:,:,t]*np.cos(FyG[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(FyG[:,:,:,t])
Sz2G[:,:,:,t] = Mz[:,:,:,t]*np.cos(FzG[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(FzG[:,:,:,t]) Sz2G[:,:,:,t] = Mz[:,:,:,t]*np.cos(FzG[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(FzG[:,:,:,t])
# Adding noise by cartesian components # Adding noise by cartesian components
for t in range(Nt): for t in range(Nt):
noise1 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) noise1 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz])
noise2 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) noise2 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz])
Sx20[:,:,:,t] += noise1 + 1j*noise2 Sx20[:,:,:,t] += noise1 + 1j*noise2
Sy20[:,:,:,t] += noise1 + 1j*noise2 Sy20[:,:,:,t] += noise1 + 1j*noise2
Sz20[:,:,:,t] += noise1 + 1j*noise2 Sz20[:,:,:,t] += noise1 + 1j*noise2

92
codes/PostCheck.py
View File

@ -175,6 +175,86 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname,
for bb in bnds: for bb in bnds:
QQ[bb] = [] QQ[bb] = []
if mode == 'perturbation':
u = Function(W)
unew = Function(W)
u.rename('velocity', outname)
unew.rename('velocity', outname)
if options['Perturbation']['xdmf']:
xdmf_u = XDMFFile(output_path+'u.xdmf')
if not options['Perturbation']['type']['SNR']=='inf':
Noise = True
def Add_Noise(signal,SNR):
Psignal = signal**2
Psignal_av = np.mean(Psignal)
Psignal_av_db = 10*np.log10(Psignal_av)
Pnoise_av_db = Psignal_av_db - SNR
Pnoise_av = 10**(Pnoise_av_db/10)
noise_std = np.sqrt(Pnoise_av)
noise = np.random.normal(0,noise_std,len(signal))
return signal + noise
else:
Noise = False
if not options['Perturbation']['type']['phase_contrast']==0:
Phase_Contrast = True
else:
Phase_Contrast = False
noise_in_coil = options['Perturbation']['type']['coil']
for k in indexes:
path = checkpoint_path + str(k) + '/'+filename+'.h5'
hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r')
hdf.read(u, 'u/vector_0')
time = hdf.attributes('u/vector_0').to_dict()['timestamp']
hdf.close()
uvec = u.vector().get_local()
if Phase_Contrast:
ufactor = options['Perturbation']['type']['phase_contrast']/100
VENC = np.max(np.abs(uvec))*ufactor
gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei
#B0 = 1.5 # Tesla Magnetic Field Strenght
TE = 5e-3 # Echo-time
Phi1 = gamma*10*TE + 0*uvec
Phi2 = gamma*10*TE + np.pi*uvec/VENC
M1 = np.cos(Phi1) + 1j*np.sin(Phi1)
M2 = np.cos(Phi2) + 1j*np.sin(Phi2)
if noise_in_coil:
SNR = options['Perturbation']['type']['SNR']
m1x_new = Add_Noise(np.real(M1),SNR)
m1y_new = Add_Noise(np.imag(M1),SNR)
m2x_new = Add_Noise(np.real(M2),SNR)
m2y_new = Add_Noise(np.imag(M2),SNR)
M1_new = m1x_new + 1j*m1y_new
M2_new = m2x_new + 1j*m2y_new
uvec = (np.angle(M2_new) - np.angle(M1_new))*VENC/np.pi
unew.vector()[:] = uvec
else:
uvec = (np.angle(M2) - np.angle(M1))*VENC/np.pi
else:
if noise_in_coil:
raise Exception('In order to perturb in coils some PC should be selected')
if not noise_in_coil:
if Noise:
SNR = options['Perturbation']['type']['SNR']
unew.vector()[:] = Add_Noise(uvec,SNR)
else:
unew.vector()[:] = uvec
print('Writing checkpoint number ',k)
write_path = output_path + 'checkpoint/{i}/'.format(i=k)
hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w')
hdf2.write(unew, '/u', time)
hdf2.close()
if options['Perturbation']['xdmf']:
xdmf_u.write(unew, time)
if mode == 'interpolation': if mode == 'interpolation':
dt_new = options['Temporal-Interpolation']['dt_new'] dt_new = options['Temporal-Interpolation']['dt_new']
@ -239,7 +319,6 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname,
if options['Temporal-Interpolation']['xdmf']: if options['Temporal-Interpolation']['xdmf']:
xdmf_u.write(unew, (k-1)*dt_new) xdmf_u.write(unew, (k-1)*dt_new)
if mode == 'average': if mode == 'average':
N_av = options['Temporal-Average']['subsampling_rate'] N_av = options['Temporal-Average']['subsampling_rate']
dt = options['Temporal-Average']['dt'] dt = options['Temporal-Average']['dt']
@ -1047,7 +1126,16 @@ def ROUTINE(options):
ref_check = options['Temporal-Interpolation']['original_check'] + 'checkpoint/' ref_check = options['Temporal-Interpolation']['original_check'] + 'checkpoint/'
out_check = options['Temporal-Interpolation']['out_check'] out_check = options['Temporal-Interpolation']['out_check']
READcheckpoint(MESH,'interpolation', out_check,ref_check,'u','u',options) READcheckpoint(MESH,'interpolation', out_check,ref_check,'u','u',options)
if 'Perturbation' in options:
if options['Perturbation']['apply']:
if rank==0:
print('--- Perturbation in measurements ---')
MESH = LOADmesh(options['Perturbation']['meshpath'])
checkpath = options['Perturbation']['checkpath'] + 'checkpoint/'
out_check = options['Perturbation']['checkpath'] + 'Perturbation/'
READcheckpoint(MESH,'perturbation', out_check,checkpath,'u','u',options)

327
codes/cfd_mono.py
View File

@ -1,327 +0,0 @@
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import dolfin
from common import inout
from mpi4py import MPI
import sys
import os
#
# NAVIER STOKES PROBLEM IN THE AORTA with a MONOLITHIC SOLVER
# THIS SCRIPT INCLUDE THE 0-WINDKESSEL BOUNDARY CONDITION
#
# Written by Jeremias Garay L: j.e.garay.labra@rug.nl
#
if '2017' in dolfin.__version__:
class MPI(MPI):
comm_world = mpi_comm_world()
set_log_active(False)
else:
parameters["std_out_all_processes"] = False
def BACH():
import os
tempof = 0.7
semicorchea = 0.1*tempof
corchea = 0.2*tempof
negra = 0.4*tempof
blanca = 0.8*tempof
LA3 = 220
SIb = 233.08
SI = 246.94
DO = 261
REb = 277.18
RE = 293.66
FA = 349.23
MIb = 311.13
MI = 329.63
SOL = 391
LA = 440
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, DO))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, LA))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, SOL))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, FA))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (semicorchea, MI))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (semicorchea, FA))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (corchea, SOL))
os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (blanca, FA))
def flux(u,n,ds,idd):
Q = 0
for k in idd:
Q += assemble(dot(u,n)*ds(k))
return Q
def save_outlet_data(options,Qin,Qout3,Qout4,Qout5,Qout6,Pin,Pout3,Pout4,Pout5,Pout6,dt):
# saving the fluxes
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q_in_dt' + str(dt) + '.txt', Qin)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q3_dt' + str(dt) + '.txt', Qout3)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q4_dt' + str(dt) + '.txt', Qout4)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q5_dt' + str(dt) + '.txt', Qout5)
np.savetxt(options['outlets_path'] + 'fluxes/' + 'q6_dt' + str(dt) + '.txt', Qout6)
# saving the pressures
np.savetxt(options['outlets_path'] + 'pressures/' + 'p_in_dt' + str(dt) + '.txt', Pin)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p3_dt' + str(dt) + '.txt', Pout3)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p4_dt' + str(dt) + '.txt', Pout4)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p5_dt' + str(dt) + '.txt', Pout5)
np.savetxt(options['outlets_path'] + 'pressures/' + 'p6_dt' + str(dt) + '.txt', Pout6)
def windkessel_update(u0,n,ds,fluxes,press,pii,pii0,windkessel):
# Updating the time dependent windkessel parameters
if windkessel['C']:
for nk in windkessel['id']:
k = str(nk)
fluxes[k].assign(assemble(dot(u0,n)*ds(nk)))
pii0[k].assign(pii[k])
pii[k].assign(windkessel['alpha'][k]*pii0[k] + windkessel['beta'][k]*fluxes[k])
press[k].assign(windkessel['gamma'][k]*fluxes[k] + windkessel['alpha'][k]*pii0[k])
else:
for nk in windkessel['id']:
k = str(nk)
fluxes[k] = assemble(dot(u0,n)*ds(nk))
press[k].assign( windkessel['R_p'][k]*assemble(dot(u0,n)*ds(nk)))
def solv_NavierStokes(options):
# Assign physical parameters
rho = Constant(options['density'])
mu = Constant(options['dynamic_viscosity'])
otheta = Constant(1-options['param']['theta'])
theta = Constant(options['param']['theta'])
Tf = options['Tf']
dt = options['dt']
dt_w = options['dt_write']
Qin = np.zeros([int(Tf/dt)])
Qout3 = np.zeros([int(Tf/dt)])
Qout4 = np.zeros([int(Tf/dt)])
Qout5 = np.zeros([int(Tf/dt)])
Qout6 = np.zeros([int(Tf/dt)])
Pin = np.zeros([int(Tf/dt)])
Pout3 = np.zeros([int(Tf/dt)])
Pout4 = np.zeros([int(Tf/dt)])
Pout5 = np.zeros([int(Tf/dt)])
Pout6 = np.zeros([int(Tf/dt)])
barye2mmHg = 1/1333.22387415
# CREATING THE FILES
xdmf_u = XDMFFile(options['savepath']+'u.xdmf')
xdmf_p = XDMFFile(options['savepath']+'p.xdmf')
xdmf_u.parameters['rewrite_function_mesh'] = False
xdmf_p.parameters['rewrite_function_mesh'] = False
# LOADING THE MESH
#mesh = Mesh(options['mesh_path'])
#boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
#boundaries = MarkBoundaries(boundaries)
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), options['mesh_path'] , 'r')
hdf.read(mesh, '/mesh', False)
boundaries = MeshFunction('size_t', mesh , mesh.topology().dim() - 1)
hdf.read(boundaries, '/boundaries')
# To save the boundaries information
#path2 = '/home/p283370/Desktop/marked_fine/boundaries.xdmf'
#XDMFFile(path2).write(boundaries)
# DEFINE FUNCTION SPACES
V = VectorElement('Lagrange', mesh.ufl_cell(), options['param']['Nvel'])
Q = FiniteElement('Lagrange', mesh.ufl_cell(), options['param']['Npress'])
TH = V * Q
W = FunctionSpace(mesh, TH)
# No-slip boundary condition for velocity on walls
noslip = Constant((0, 0, 0))
inflow = Expression(('0','0','(-0.5*U*fabs(sin(w*t)) - 0.5*U*sin(w*t))*(t<0.7)'), degree=3,t=0,U=options['param']['U'], w=2*DOLFIN_PI/options['param']['period'])
bc_inflow = DirichletBC(W.sub(0), inflow, boundaries, 2)
bc_walls = DirichletBC(W.sub(0), noslip, boundaries, 1)
# COLLECTING THE BOUNDARY CONDITIONS
BCS = [bc_inflow,bc_walls]
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
w = Function(W)
n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data=boundaries)
h = CellDiameter(mesh)
f = Constant((0,0,0))
u0, p0 = w.split()
u_ = theta*u + (1 - theta)*u0
theta_p = theta
p_ = theta_p*p + (1 - theta_p)*p0
# The variational formulation
# Mass matrix
F = (
(rho/dt)*dot(u - u0, v)*dx
+ mu*inner(grad(u_), grad(v))*dx
- p_*div(v)*dx + q*div(u)*dx
+ rho*dot(grad(u_)*u0, v)*dx
)
if options['stab']['temam']:
F += 0.5*rho*div(u0)*dot(u_, v)*dx
if options['stab']['pspg']:
eps = Constant(options['stab']['epsilon'])
F += eps/mu*h**2*inner(grad(p_), grad(q))*dx
if options['stab']['backflow']:
def abs_n(x):
return 0.5*(x - abs(x))
for nk in options['stab']['back_id']:
if rank==0:
print('adding backflow stabilization in border number:' + str(nk))
F -= 0.5*rho*abs_n(dot(u0, n))*dot(u_, v)*ds(nk)
a = lhs(F)
L = rhs(F)
# Initialization of Windkessel Boundaries
if options['windkessel']['id']:
windkessel = options['windkessel']
# Coeficients
fluxes = {str(k):[] for k in windkessel['id']}
pii0 = {str(k):[] for k in windkessel['id']}
pii = {str(k):[] for k in windkessel['id']}
press = {str(k):[] for k in windkessel['id']}
if windkessel['C']:
for nk in windkessel['id']:
k = str(nk)
if rank==0:
print('Capacitance of windkessel model is: ', windkessel['C'][k])
# computing coeficients
windkessel['alpha'][k] = windkessel['R_d'][k]*windkessel['C'][k]/(windkessel['R_d'][k]*windkessel['C'][k] + options['dt'])
windkessel['beta'][k] = windkessel['R_d'][k]*(1-windkessel['alpha'][k])
windkessel['gamma'][k] = windkessel['R_p'][k] + windkessel['beta'][k]
if rank==0:
print('Using 0-Windkessel complete at outlet number: ' + str(k))
# setting initial values for flux, distal pressure and proximal pressure
fluxes[k] = Constant(0)
pii0[k] = Constant(47/barye2mmHg)
pii[k] = Constant(windkessel['alpha'][k]*pii0[k] + windkessel['beta'][k]*fluxes[k])
press[k] = Constant(windkessel['gamma'][k]*fluxes[k] + windkessel['alpha'][k]*pii0[k])
# Adding to RHS
L = L - dt*press[k]*dot(v,n)*ds(nk)
else:
for nk in windkessel['id']:
k = str(nk)
if rank==0:
print('Using 0-Windkessel reduced at outlet number: ' + str(nk))
fluxes[k] = Constant(0)
press[k] = Constant(windkessel['R_p'][k]*0)
# Adding to RHS
L = L - dt*press[k]*dot(v,n)*ds(nk)
# The static part of the matrix
A = assemble(a)
#b = assemble(L)
[bc.apply(A) for bc in BCS]
#[bc.apply(b) for bc in BCS]
#solv = LUSolver()
#solv.set_operator(A)
#solv.parameters['linear_solver'] = 'mumps'
#solv.parameters['reuse_factorization'] = True
u, p = w.split()
u.rename('velocity', 'u')
p.rename('pressure', 'p')
ind = 0
t = dt
ones = interpolate(Constant(1),W.sub(1).collapse())
A2 = assemble(ones*ds(2))
A3 = assemble(ones*ds(3))
A4 = assemble(ones*ds(4))
A5 = assemble(ones*ds(5))
A6 = assemble(ones*ds(6))
checkcicle = int(options['checkpoint_dt']/options['dt'])
writecicle = int(options['checkpoint_dt']/options['dt_write'])
while t<=Tf+dt:
if options['windkessel']['id']:
windkessel_update(u,n,ds,fluxes,press,pii,pii0,windkessel)
# To solve
assemble(a, tensor=A)
b = assemble(L)
[bc.apply(A, b) for bc in BCS]
solve(A, w.vector(), b)
Qin[ind] = flux(u,n,ds,[2])
Qout3[ind] = flux(u,n,ds,[3])
Qout4[ind] = flux(u,n,ds,[4])
Qout5[ind] = flux(u,n,ds,[5])
Qout6[ind] = flux(u,n,ds,[6])
Pin[ind] = barye2mmHg*assemble(p*ds(2))/A2
Pout3[ind] = barye2mmHg*assemble(p*ds(3))/A3
Pout4[ind] = barye2mmHg*assemble(p*ds(4))/A4
Pout5[ind] = barye2mmHg*assemble(p*ds(5))/A5
Pout6[ind] = barye2mmHg*assemble(p*ds(6))/A6
if rank==0:
print('t = ',t)
# print('|u|:', norm(u0))
# print('|p|:', norm(p0))
# print('div(u):', assemble(div(u0)*dx))
print('Dp = ',np.round(Pin[ind]-Pout3[ind],3),'mmHg')
ind += 1
if options['write_xdmf']:
if np.mod(ind,writecicle)<0.1 or ind==1:
xdmf_u.write(u, t)
xdmf_p.write(p, t)
if np.mod(ind,checkcicle)<0.1 or ind==1:
if options['write_checkpoint']:
checkpath = options['savepath'] +'checkpoint/{i}/'.format(i=ind)
comm = u.function_space().mesh().mpi_comm()
inout.write_HDF5_data(comm, checkpath + '/u.h5', u, '/u', t=t)
inout.write_HDF5_data(comm, checkpath + '/p.h5', p, '/p', t=t)
inflow.t = t
t += dt
# saving the data at outlets: fluxes and pressures
if options['write_outlets']:
if rank==0:
save_outlet_data(options,Qin,Qout3,Qout4,Qout5,Qout6,Pin,Pout3,Pout4,Pout5,Pout6,options['dt'])
if __name__ == '__main__':
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
if len(sys.argv) > 1:
if os.path.exists(sys.argv[1]):
inputfile = sys.argv[1]
if rank==0:
print('Found input file ' + inputfile)
else:
raise Exception('Command line arg given but input file does not exist:'
' {}'.format(sys.argv[1]))
else:
raise Exception('An input file is required as argument!')
if 'Zion' in os.popen('hostname').read():
user = 'yeye'
np.set_printoptions(threshold=5)
if 'fwn-bborg-5-166' in os.popen('hostname').read():
user = 'p283370'
if rank==0:
print('Welcome user {uss}'.format(uss=user))
options = inout.read_parameters(inputfile)
solv_NavierStokes(options)

216
codes/monolithic.py Normal file
View File

@ -0,0 +1,216 @@
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
from common import inout
from mpi4py import MPI
import sys
import os
#
# NAVIER STOKES PROBLEM IN THE AORTA with a MONOLITHIC SOLVER
# THIS SCRIPT INCLUDE THE 0-WINDKESSEL BOUNDARY CONDITION
#
# Written by Jeremias Garay L: j.e.garay.labra@rug.nl
#
parameters["std_out_all_processes"] = False
def solv_NavierStokes(options):
# Assign physical parameters
rho = Constant(options['density'])
mu = Constant(options['dynamic_viscosity'])
otheta = Constant(1-options['theta'])
theta = Constant(options['theta'])
Tf = options['Tf']
dt = options['dt']
dt_w = options['dt_write']
# CREATING THE FILES
xdmf_u = XDMFFile(options['savepath']+'u.xdmf')
xdmf_p = XDMFFile(options['savepath']+'p.xdmf')
# LOADING THE MESH
mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(), options['mesh_path'] , 'r')
hdf.read(mesh, '/mesh', False)
bnds = MeshFunction('size_t', mesh , mesh.topology().dim() - 1)
hdf.read(bnds, '/boundaries')
# DEFINE FUNCTION SPACES
if options['fem_space'] == 'p2p1':
V = VectorElement('Lagrange', mesh.ufl_cell(), 2)
Q = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
pspg = False
elif options['fem_space'] == 'p1p1':
V = VectorElement('Lagrange', mesh.ufl_cell(), 1)
Q = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
pspg = True
TH = V * Q
W = FunctionSpace(mesh, TH)
n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data=bnds)
v, q = TestFunctions(W)
# Boundary Conditions
BCS = []
bc = options['boundary_conditions']
# For Windkessel implementation
flows = {}
pii0 = {}
pii = {}
press = {}
alpha = {}
beta = {}
gamma = {}
Windkvar = {}
windkessel = False
for nbc in range(len(bc)):
nid = bc[nbc]['id']
if bc[nbc]['type'] == 'dirichlet':
if rank==0:
print('Adding Dirichlet BC at boundary ', nid)
val = bc[nbc]['value']
if isinstance(val, (int, float)):
val = Constant(val)
BCS.append(DirichletBC(W.sub(0), val, bnds, nid))
else:
params = bc[nbc]['parameters'] if 'parameters' in bc[nbc] else dict()
inflow = Expression(val, degree=3, **params)
BCS.append(DirichletBC(W.sub(0), inflow, bnds, nid))
if bc[nbc]['type'] == 'windkessel':
windkessel = True
if rank==0:
print('Adding Windkessel BC at boundary ',nid)
[R_p,C,R_d] = bc[nbc]['value']
# coeficients
alpha[nid] = R_d*C/(R_d*C + dt)
beta[nid] = R_d*(1-alpha[nid])
gamma[nid] = R_p + beta[nid]
# dynamical terms
if C ==0:
flows[nid] = Constant(0)
press[nid] = Constant(R_p*0)
else:
flows[nid] = Constant(0)
pii0[nid] = Constant(bc[nbc]['p0'][0]*bc[nbc]['p0'][1])
pii[nid] = Constant(alpha[nid]*pii0[nid] + beta[nid]*flows[nid])
press[nid] = Constant(gamma[nid]*flows[nid] + alpha[nid]*pii0[nid])
Windkvar[nid] = dt*press[nid]*dot(v,n)*ds(nid)
u, p = TrialFunctions(W)
w = Function(W)
h = CellDiameter(mesh)
u0, p0 = w.split()
u_ = theta*u + otheta*u0
p_ = theta*p + otheta*p0
# The variational formulation
# Mass matrix
F = (
(rho/dt)*dot(u - u0, v)*dx
+ mu*inner(grad(u_), grad(v))*dx
- p_*div(v)*dx + q*div(u)*dx
+ rho*dot(grad(u_)*u0, v)*dx
)
# Stabilization Terms
if options['stabilization']['temam']:
if rank==0:
print('Addint Temam stabilization term')
F += 0.5*rho*div(u0)*dot(u_, v)*dx
if pspg:
if rank==0:
print('Adding PSPG stabilization term')
eps = Constant(options['stabilization']['eps'])
F += eps/mu*h**2*inner(grad(p_), grad(q))*dx
if options['stabilization']['forced_normal']['enabled']:
gparam = options['stabilization']['forced_normal']['gamma']
for nid in options['stabilization']['forced_normal']['boundaries']:
if rank==0:
print('Forcing normal velocity in border ', nid)
ut = u - n*dot(u,n)
vt = v - n*dot(v,n)
F += gparam*dot(ut,vt)*ds(nid)
if len(options['stabilization']['backflow_boundaries'])>0:
def abs_n(x):
return 0.5*(x - abs(x))
for nk in options['stabilization']['backflow_boundaries']:
if rank==0:
print('adding backflow stabilization in border number:' + str(nk))
F -= 0.5*rho*abs_n(dot(u0, n))*dot(u_, v)*ds(nk)
a = lhs(F)
L = rhs(F)
if windkessel:
for nid in Windkvar.keys():
L = L - Windkvar[nid]
# The static part of the matrix
A = assemble(a)
u, p = w.split()
u.rename('velocity', 'u')
p.rename('pressure', 'p')
ind = 0
t = dt
checkcicle = int(options['checkpoint_dt']/options['dt'])
writecicle = int(options['checkpoint_dt']/options['dt_write'])
while t<=Tf+dt:
if windkessel:
for k in flows.keys():
flows[k].assign(assemble(dot(u0,n)*ds(k)))
pii0[k].assign(pii[k])
pii[k].assign(alpha[k]*pii0[k] + beta[k]*flows[k])
press[k].assign(gamma[k]*flows[k] + alpha[k]*pii0[k])
# To solve
assemble(a, tensor=A)
b = assemble(L)
[bcs.apply(A, b) for bcs in BCS]
solve(A, w.vector(), b , 'lu')
ind += 1
if options['write_xdmf']:
if np.mod(ind,writecicle)<0.1 or ind==1:
xdmf_u.write(u, t)
xdmf_p.write(p, t)
if np.mod(ind,checkcicle)<0.1 or ind==1:
if options['write_checkpoint']:
checkpath = options['savepath'] +'checkpoint/{i}/'.format(i=ind)
comm = u.function_space().mesh().mpi_comm()
inout.write_HDF5_data(comm, checkpath + '/u.h5', u, '/u', t=t)
inout.write_HDF5_data(comm, checkpath + '/p.h5', p, '/p', t=t)
inflow.t = t
t += dt
assign(u0, w.sub(0))
if __name__ == '__main__':
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
if len(sys.argv) > 1:
if os.path.exists(sys.argv[1]):
inputfile = sys.argv[1]
if rank==0:
print('Found input file ' + inputfile)
else:
raise Exception('Command line arg given but input file does not exist:'
' {}'.format(sys.argv[1]))
else:
raise Exception('An input file is required as argument!')
options = inout.read_parameters(inputfile)
solv_NavierStokes(options)

BIN
corrector/vol/colormap1.mp4
View File

Binary file not shown.

BIN
corrector/vol/cor1.mp4
View File

Binary file not shown.

BIN
corrector/vol/max.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 310 KiB

BIN
corrector/vol/mean.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 290 KiB

BIN
corrector/vol/obs1.mp4
View File

Binary file not shown.