sonicLiquid.py 5.18 KB
Newer Older
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)