import numpy, treelog, nutils import matplotlib.pyplot as plt from scipy import special from nutils import mesh, function, solver, cli, export import math import precice from mpi4py import MPI def main(degree: int, nelems_z: int, Omega: float, Ro: float, gap: float, reynolds: float): ''' .. arguments:: degree [2] nelems_z [15] Omega [1.0] Ro [1.0] gap [0.5] reynolds [1.0] dt = 1 ''' # Define/allocate constants degree = 2 nelems_z = 25 nelems_r = 25 Omega = 1 gap = 0.005 R0 = 1e-10 # Ideally equal to 0 -> causes instability w_rot = 0.04 w_bi = 0.007 w_m = 0.008 Ro_c = 0.158 Ri_c = 0.088 Rout = 0.17 pi = numpy.pi ns = function.Namespace() ns.pi = pi dt = 1 r_gap = 0.005 L_rotor = 0.17 W_rotor = 0.04 verts_r = numpy.linspace(R0, Rout, nelems_r + 1) verts_z = numpy.linspace(r_gap, r_gap + w_rot, nelems_z + 1) h = gap / nelems_z # Characteristic element length domain, (r, θ, z) = mesh.rectilinear([verts_r, [-2.5*numpy.pi/180, 2.5*numpy.pi/180], verts_z], periodic=[1]) rdir = function.stack([ function.cos(θ), 0,function.sin(θ)]) # outward radial zdir = function.stack([0, 1, 0]) # axial direction geom = r * rdir + z * zdir ns.x = geom degree = 2 # Specify degree of bsis function ns.basis = domain.basis('spline',degree=degree) ns.u = 'basis_n ?lhs_n' ns.q = -100 ns.flux = 'basis_n ?fluxdofs_n' ns.k = 20.0 ns.km = 7.6 ns.kbi = 40.0 ns.TK = 273.15 ns.Tin = '350' #h_fraction_Ro = Ro_c / Rout #h_fraction_Ri = Ri_c / Rout #w_fraction_bi = w_bi / w_rot #w_fraction_m = w_m / w_rot #N_fraction_Ro = math.ceil(h_fraction_Ro * nelems_r) #N_fraction_Ri = math.ceil(h_fraction_Ri * nelems_r) #N_fraction_bi = math.ceil((w_fraction_bi + w_fraction_m) * nelems_z) #N_fraction_m = math.ceil(w_fraction_m * nelems_z) #NameSquare = 'Magnet' #idz1 = nelems_z + 1 - N_fraction_m # Nelemsr #idz2 = nelems_z + 1 # Nelemsr #idy1 = 0 #idy2 = 1 #idx1 = N_fraction_Ri # Nelemsz #idx2 = N_fraction_Ro # Nelemsz #domain = domain.withgroups(vgroups={NameSquare: domain[idx1:idx2, idy1:idy2, idz1:idz2]}) #NameSquare2 = 'Backiron' #idz1 = nelems_z + 1 - N_fraction_bi # Nelemsr #idz2 = nelems_z + 1 - N_fraction_m # Nelemsr #idy1 = 0 #idy2 = 1 #idx1 = N_fraction_Ri # Nelemsz #idx2 = N_fraction_Ro # Nelemsz #domain = domain.withgroups(vgroups={NameSquare2: domain[idx1:idx2, idy1:idy2, idz1:idz2]}) #the weak form res = domain.integral('( k basis_n,i u_,i ) d:x' @ ns, degree=degree * 2) #res += domain[NameSquare].integral('( basis_n Q) d:x' @ ns, degree=2 * degree) #res += domain[NameSquare2].integral('( kbi basis_n,i u_,i ) d:x' @ ns, degree=2 * degree) # Dirichlet boundary condition sqr = domain.boundary['left'].integral('(u - Tin)^2 d:x' @ ns, degree=2 * degree) cons = solver.optimize('lhs', sqr, droptol=1e-15) # preCICE setup configFileName = "../precice-config.xml" participantName = "Nutils" solverProcessIndex = 0 solverProcessSize = 1 interface = precice.Interface(participantName, configFileName, solverProcessIndex, solverProcessSize) # define coupling meshes meshNameGP = "Nutils-Mesh-GP" # Gauss points meshNameCC = "Nutils-Mesh-CC" # cell centers (potentially sub-sampled) meshIDGP = interface.get_mesh_id(meshNameGP) meshIDCC = interface.get_mesh_id(meshNameCC) Bnames = 'back,front,right' domain = domain.withboundary(MyBoundary=Bnames) couplinginterface = domain.boundary['MyBoundary'] couplingsampleGP = couplinginterface.sample('gauss', degree=degree * 2) couplingsampleCC = couplinginterface.sample('uniform', 4) # number of sub-samples for better mapping verticesGP = couplingsampleGP.eval(ns.x) verticesCC = couplingsampleCC.eval(ns.x) dataIndicesGP = interface.set_mesh_vertices(meshIDGP, verticesGP) dataIndicesCC = interface.set_mesh_vertices(meshIDCC, verticesCC) # coupling data writeData = "Heat-Flux" readData = "Temperature" writedataID = interface.get_data_id(writeData, meshIDCC) readdataID = interface.get_data_id(readData, meshIDGP) # heat flux computation projectionmatrix = couplinginterface.integrate(ns.eval_nm('basis_n basis_m d:x'), degree=degree * 2) projectioncons = numpy.zeros(res.shape) projectioncons[projectionmatrix.rowsupp(1e-15)] = numpy.nan fluxdofs = lambda v: projectionmatrix.solve(v, constrain=projectioncons) precice_dt = interface.initialize() cons0 = cons # to not lose the Dirichlet BC at the bottom timestep = 1 while interface.is_coupling_ongoing(): # read temperature from interface if interface.is_read_data_available(): readdata = interface.read_block_scalar_data(readdataID, dataIndicesGP) coupledata = couplingsampleGP.asfunction(readdata) sqr = couplingsampleGP.integral((ns.u - coupledata)**2) cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons0) # save checkpoint if interface.is_action_required(precice.action_write_iteration_checkpoint()): interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint()) # potentially adjust non-matching timestep sizes dt = min(dt, precice_dt) # solve nutils timestep lhs = nutils.solver.solve_linear('lhs', res, constrain=cons) # write heat fluxes to interface if interface.is_write_data_required(dt): fluxvalues = res.eval(lhs=lhs) writedata = couplingsampleCC.eval('-flux' @ ns, fluxdofs=fluxdofs(fluxvalues)) interface.write_block_scalar_data(writedataID, dataIndicesCC, writedata) # do the coupling precice_dt = interface.advance(dt) # read checkpoint if required if interface.is_action_required(precice.action_read_iteration_checkpoint()): interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint()) else: # go to next timestep and visualize bezier = domain.sample('bezier', 2) x, u = bezier.eval(['x_i', 'u'] @ ns, lhs=lhs) print(u) with treelog.add(treelog.DataLog()): if timestep % 1 == 0: nutils.export.vtk('Solid_' + str(timestep), bezier.tri, x, T=u) timestep += 1 interface.finalize() if __name__ == '__main__': nutils.cli.run(main)