sonicLiquid.py 5.18 KB
 Benjamin Uekermann committed Jul 05, 2019 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 ``````#! /usr/bin/python3 from nutils import * import numpy, itertools, treelog import precice_future as precice #import precice from mpi4py import MPI @function.replace def subs0(f): if isinstance(f, function.Argument) and f._name == 'lhs': return function.Argument(name='lhs0', shape=f.shape, nderiv=f._nderiv) def main(nelems=200, dt=.005, refdensity=1e3, refpressure=101325., psi=1e-6, viscosity=1e-3, theta=.55, side = 'Left'): #domain, geom = mesh.rectilinear([numpy.linspace(0, 1000, nelems+1), numpy.linspace(0, 1, 6)]) domain, geom = mesh.rectilinear([numpy.linspace(0, 1000, nelems+1)]) if side == 'Left': domain = domain[:nelems//2] elif side == 'Right': domain = domain[nelems//2:] else: raise Exception('invalid side {!r}'.format(side)) bezier = domain.sample('bezier', 2) ns = function.Namespace() ns._functions['t'] = lambda f: theta * f + (1-theta) * subs0(f) ns._functions_nargs['t'] = 1 ns._functions['δt'] = lambda f: (f - subs0(f)) / dt ns._functions_nargs['δt'] = 1 ns.x = geom ns.ρref = refdensity ns.pref = refpressure ns.pin = 98100 ns.μ = viscosity ns.ψ = psi ns.ρbasis, ns.ubasis = function.chain([domain.basis('std', degree=1), domain.basis('std', degree=2).vector(domain.ndims)]) ns.ρ = 'ρref + ρbasis_n ?lhs_n' # density is shifted by ref density ns.u_i = 'ubasis_ni ?lhs_n' ns.p = 'pref + (ρ - ρref) / ψ' # pressure-density connection, see equ (8) ns.σ_ij = 'μ (u_i,j + u_j,i) - p δ_ij' # diffusive term and pressure gradient ns.h = 1 / nelems ns.k = 'ρ h / μ' # needs work, stabilization coeff res = domain.integral('ρbasis_n (δt(ρ) + t((ρ u_k)_,k)) d:x' @ ns, degree=4) # mass balance res += domain.integral('(ubasis_ni (δt(ρ u_i) + t((ρ u_i u_j)_,j)) + ubasis_ni,j t(σ_ij)) d:x' @ ns, degree=4) # momentum balance if side == "Left": res += domain.boundary['left'].integral('pin ubasis_ni n_i d:x' @ ns, degree=4) # pressure set at inlet #res += domain.integral('k ubasis_ni,j (u_j / sqrt(u_k u_k)) (δρu_i / δt + (ρ u_i u_j - σ_ij)_,j) d:x' @ ns, degree=4) # SUPG stabilization configFileName = "precice-config.xml" participantName = "Nutils" + side meshName = participantName + "-Mesh" solverProcessIndex = 0 solverProcessSize = 1 interface = precice.Interface(participantName, solverProcessIndex, solverProcessSize) interface.configure(configFileName) meshID = interface.get_mesh_id(meshName) writeData = "Pressure" if side == "Left" else "Velocity" readData = "Velocity" if side == "Left" else "Pressure" writedataID = interface.get_data_id(writeData, meshID) readdataID = interface.get_data_id(readData, meshID) couplinginterface = domain.boundary['right' if side == 'Left' else 'left'] vertex = numpy.array([0.0, 500.0, 0.0]) dataIndex = interface.set_mesh_vertex(meshID, vertex) precice_dt = interface.initialize() lhs = numpy.zeros(res.shape) t = 0 n = 0 f = open("watchpoint" + side + ".txt", "w") res0 = res lhs0 = numpy.zeros(res.shape) while interface.is_coupling_ongoing(): with log.context('dt', n): #export.triplot('output/pressure.png', x, p, tri=bezier.tri) #export.triplot('output/denity.png', x, ρ, tri=bezier.tri) #export.triplot('output/velocity.png', x, u[:,0], tri=bezier.tri) #interface.read_scalar_data(readdataID, dataIndex[0], readdata) readdata = interface.read_vector_data(readdataID, dataIndex) if side == "Left": uOut2 = readdata[1] else: uOut2 = min(.2 * t, 1) ns.pin = readdata sqr = domain.boundary['right'].integral('(u_0 - ?uOut)^2' @ ns, degree=4) cons = solver.optimize('lhs', sqr, arguments=dict(uOut=uOut2), droptol=1e-14) if side == "Right": #sqr = domain.boundary['left'].integral('(p - ?pIN)^2' @ ns, degree=4) #cons += solver.optimize('lhs', sqr, arguments=dict(pIN=readdata[0]), droptol=1e-14) res = res0 + domain.boundary['left'].integral('pin ubasis_ni n_i d:x' @ ns, degree=4) # save checkpoint if interface.is_action_required(precice.action_write_iteration_checkpoint()): lhscheckpoint = lhs0 interface.fulfilled_action(precice.action_write_iteration_checkpoint()) lhs = solver.newton('lhs', res, constrain=cons, arguments=dict(lhs0=lhs0), lhs0=lhs0).solve(1e-8) x, p, ρ, u = bezier.eval(['x_i', 'p', 'ρ', 'u_i'] @ ns, lhs=lhs) writedata = u[0] if side == "Right" else p[199] interface.write_scalar_data(writedataID, dataIndex, writedata) precice_dt = interface.advance(dt) # read checkpoint if required if interface.is_action_required(precice.action_read_iteration_checkpoint()): interface.fulfilled_action(precice.action_read_iteration_checkpoint()) lhs0 = lhscheckpoint else: x, p, ρ, u = bezier.eval(['x_i', 'p', 'ρ', 'u_i'] @ ns, lhs=lhs) f.write("%e; %e; %e; %e; %e; %e; %e\n" % (t, p[0], u[0], p[199], u[199], p[99], u[99])) f.flush() t += dt n += 1 lhs0 = lhs interface.finalize() f.close() if __name__ == '__main__': cli.run(main)``````