diff --git a/codes/MRI.py b/codes/MRI.py index 8ebba11..a1eab4b 100644 --- a/codes/MRI.py +++ b/codes/MRI.py @@ -18,13 +18,9 @@ from mpi4py import MPI from common import inout import time -if '2017' in dolfin.__version__: - #class MPI(MPI): - # comm_world = mpi_comm_world() - set_log_active(False) -else: - PRINT_BARS = False - parameters["std_out_all_processes"] = False + +PRINT_BARS = False +parameters["std_out_all_processes"] = False def Etcetera(mode): if mode=='happy': @@ -50,31 +46,31 @@ def Etcetera(mode): def dealiased(VENC,vx,vy,vz): - vx2 = np.zeros(vx.shape) - vy2 = np.zeros(vy.shape) - vz2 = np.zeros(vz.shape) - - badx_pos = np.where(vx>VENC) - badx_neg = np.where(vx<-VENC) - bady_pos = np.where(vy>VENC) - bady_neg = np.where(vy<-VENC) - badz_pos = np.where(vz>VENC) - badz_neg = np.where(vz<-VENC) + vx2 = np.zeros(vx.shape) + vy2 = np.zeros(vy.shape) + vz2 = np.zeros(vz.shape) + + badx_pos = np.where(vx > VENC) + badx_neg = np.where(vx < -VENC) + bady_pos = np.where(vy > VENC) + bady_neg = np.where(vy < -VENC) + badz_pos = np.where(vz > VENC) + badz_neg = np.where(vz < -VENC) # More than VENC - vx2[badx_pos] = -(2*VENC - vx[badx_pos]) - vy2[bady_pos] = -(2*VENC - vy[bady_pos]) - vz2[badz_pos] = -(2*VENC - vz[badz_pos]) + vx2[badx_pos] = -(2*VENC - vx[badx_pos]) + vy2[bady_pos] = -(2*VENC - vy[bady_pos]) + vz2[badz_pos] = -(2*VENC - vz[badz_pos]) # Less than VENC - vx2[badx_neg] = 2*VENC + vx[badx_neg] - vy2[bady_neg] = 2*VENC + vy[bady_neg] - vz2[badz_neg] = 2*VENC + vz[badz_neg] + vx2[badx_neg] = 2*VENC + vx[badx_neg] + vy2[bady_neg] = 2*VENC + vy[bady_neg] + vz2[badz_neg] = 2*VENC + vz[badz_neg] # The rest - otherx = np.where(vx2==0) - othery = np.where(vy2==0) - otherz = np.where(vz2==0) - vx2[otherx] = vx[otherx] - vy2[othery] = vy[othery] - vz2[otherz] = vz[otherz] + otherx = np.where(vx2 == 0) + othery = np.where(vy2 == 0) + otherz = np.where(vz2 == 0) + vx2[otherx] = vx[otherx] + vy2[othery] = vy[othery] + vz2[otherz] = vz[otherz] return [vx2,vy2,vz2] @@ -1326,35 +1322,35 @@ def SCANNER(options): # Creating M 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])) - 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])) + 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])) + Mz[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqz[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqz[:,:,:,t])) # 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] + 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] + 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] + 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 - 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 - 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) - Fz0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + 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 ) + 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 ) + 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 ) + 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 + 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) + 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) - 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]) - 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]) - 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]) + 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]) + 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]) + 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]) # Adding noise by cartesian components for t in range(Nt): - noise1 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) - noise2 = 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]) Sx20[:,:,:,t] += noise1 + 1j*noise2 Sy20[:,:,:,t] += noise1 + 1j*noise2 Sz20[:,:,:,t] += noise1 + 1j*noise2 diff --git a/codes/PostCheck.py b/codes/PostCheck.py index 1674730..95fbf73 100644 --- a/codes/PostCheck.py +++ b/codes/PostCheck.py @@ -175,6 +175,86 @@ def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, for bb in bnds: 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': 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']: xdmf_u.write(unew, (k-1)*dt_new) - if mode == 'average': N_av = options['Temporal-Average']['subsampling_rate'] dt = options['Temporal-Average']['dt'] @@ -1047,7 +1126,16 @@ def ROUTINE(options): ref_check = options['Temporal-Interpolation']['original_check'] + 'checkpoint/' out_check = options['Temporal-Interpolation']['out_check'] 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) diff --git a/codes/cfd_mono.py b/codes/cfd_mono.py deleted file mode 100644 index b6d17d4..0000000 --- a/codes/cfd_mono.py +++ /dev/null @@ -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) - diff --git a/codes/monolithic.py b/codes/monolithic.py new file mode 100644 index 0000000..193d96e --- /dev/null +++ b/codes/monolithic.py @@ -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) + diff --git a/corrector/vol/colormap1.mp4 b/corrector/vol/colormap1.mp4 deleted file mode 100644 index 0a7d0c7..0000000 Binary files a/corrector/vol/colormap1.mp4 and /dev/null differ diff --git a/corrector/vol/cor1.mp4 b/corrector/vol/cor1.mp4 deleted file mode 100644 index 5330cbe..0000000 Binary files a/corrector/vol/cor1.mp4 and /dev/null differ diff --git a/corrector/vol/max.png b/corrector/vol/max.png deleted file mode 100644 index 063d5b0..0000000 Binary files a/corrector/vol/max.png and /dev/null differ diff --git a/corrector/vol/mean.png b/corrector/vol/mean.png deleted file mode 100644 index d7b160f..0000000 Binary files a/corrector/vol/mean.png and /dev/null differ diff --git a/corrector/vol/obs1.mp4 b/corrector/vol/obs1.mp4 deleted file mode 100644 index a1f0d31..0000000 Binary files a/corrector/vol/obs1.mp4 and /dev/null differ