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)