11#!/usr/bin/env python3
22
3- from numpy .core .fromnumeric import reshape
4- # TODO: uncomment for vtk output
5- # from evtk.hl import pointsToVTK
6- import numpy as np
3+ from nutils import mesh , function , solver , cli
4+ from nutils .expression_v2 import Namespace
5+ import numpy
6+ import treelog
7+ import matplotlib .pyplot as plt
78import precice
89
9- def main ():
10-
11- # number of nodes, length of domain and space interval
12-
13- length = 10
14- n = 20
15- h = 1 / n
16- t = 0
17- counter = 0
18-
19- # generate mesh
20-
21- x = np .zeros (n + 1 )
22- y = np .zeros (n + 1 )
23- z = np .linspace (0 ,length ,n + 1 )
24-
25- # initial data
26-
27- u = np .zeros ((n + 1 ,3 ))
28- p = np .zeros (n + 1 )
29- rhs = np .zeros (n + 1 )
3010
11+ def main (nelems : int , etype : str , degree : int , reynolds : float ):
12+ '''
13+ 1D channel flow problem.
14+ .. arguments::
15+ nelems [12]
16+ Number of elements along edge.
17+ etype [square]
18+ Element type (square/triangle/mixed).
19+ degree [2]
20+ Polynomial degree for velocity; the pressure space is one degree less.
21+ reynolds [1000]
22+ Reynolds number, taking the domain size as characteristic length.
23+ '''
3124 # preCICE setup
32-
3325 participant_name = "Fluid1D"
34- config_file_name = "./precice-config.xml"
26+ config_file_name = ".. /precice-config.xml"
3527 solver_process_index = 0
3628 solver_process_size = 1
3729 interface = precice .Interface (participant_name , config_file_name , solver_process_index , solver_process_size )
38-
3930 mesh_name = "Fluid1D-Mesh"
4031 mesh_id = interface .get_mesh_id (mesh_name )
41-
4232 velocity_name = "Velocity"
4333 velocity_id = interface .get_data_id (velocity_name , mesh_id )
44-
4534 pressure_name = "Pressure"
4635 pressure_id = interface .get_data_id (pressure_name , mesh_id )
47-
4836 positions = [0 , 0 , 0 ]
49-
5037 vertex_ids = interface .set_mesh_vertex (mesh_id , positions )
51-
5238 precice_dt = interface .initialize ()
5339
40+ # problem definition
41+ domain , geom = mesh .rectilinear ([numpy .linspace (0 , 1 , nelems + 1 )])
42+ ns = Namespace ()
43+ ns .δ = function .eye (domain .ndims )
44+ ns .Σ = function .ones ([domain .ndims ])
45+ ns .Re = reynolds
46+ ns .x = geom
47+ ns .define_for ('x' , gradient = '∇' , normal = 'n' , jacobians = ('dV' , 'dS' ))
48+ ns .ubasis = domain .basis ('std' , degree = 2 ).vector (domain .ndims )
49+ ns .pbasis = domain .basis ('std' , degree = 1 )
50+ ns .u = function .dotarg ('u' , ns .ubasis )
51+ ns .p = function .dotarg ('p' , ns .pbasis )
52+ ns .stress_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
53+ ures = domain .integral ('∇_j(ubasis_ni) stress_ij dV' @ ns , degree = 4 )
54+ pres = domain .integral ('pbasis_n ∇_k(u_k) dV' @ ns , degree = 4 )
5455
5556 while interface .is_coupling_ongoing ():
56- if interface .is_action_required (
57- precice .action_write_iteration_checkpoint ()):
57+ if interface .is_read_data_available (): # get dirichlet pressure outlet value from 3D solver
58+ p_read = interface .read_scalar_data (pressure_id , vertex_ids )
59+ p_read = numpy .maximum (0 , p_read ) # filter out unphysical negative pressure values
60+ else :
61+ p_read = 0
62+ usqr = domain .boundary ['left' ].integral ('(u_0 - 1)^2 dS' @ ns , degree = 2 )
63+ diricons = solver .optimize ('u' , usqr , droptol = 1e-15 )
64+ ucons = diricons
65+ stringintegral = '(p - ' + str (numpy .rint (p_read )) + ')^2 dS'
66+ psqr = domain .boundary ['right' ].integral (stringintegral @ ns , degree = 2 )
67+ pcons = solver .optimize ('p' , psqr , droptol = 1e-15 )
68+ cons = dict (u = ucons , p = pcons )
69+ with treelog .context ('stokes' ):
70+ state0 = solver .solve_linear (('u' , 'p' ), (ures , pres ), constrain = cons )
71+ x , u , p = postprocess (domain , ns , precice_dt , ** state0 )
5872
73+ if interface .is_action_required (
74+ precice .action_write_iteration_checkpoint ()):
5975 u_iter = u
6076 p_iter = p
61- t_iter = t
62- counter_iter = counter
63-
6477 interface .mark_action_fulfilled (
65- precice .action_write_iteration_checkpoint ())
66-
67- # determine time step size
68-
69- dt = 0.025
70- dt = np .minimum (dt ,precice_dt )
71-
72- # set boundary conditions
73-
74- u [0 ,2 ] = 1
75- # u[1,2] = 1 # dirichlet velocity inlet
76-
77- u [- 1 ,2 ] = u [- 2 ,2 ] # neumann velocity outlet
78-
79- p [0 ] = p [1 ] # neumann pressure inlet
80- # p[-1] = 0 # dirichlet pressure outlet
81- if interface .is_read_data_available (): # get dirichlet pressure outlet value from 3D solver
82- p_read_in = interface .read_scalar_data (pressure_id , vertex_ids )
83- p [- 1 ] = p_read_in
84- else :
85- p [- 1 ] = 0
86-
87- # compute right-hand side of 1D PPE
78+ precice .action_write_iteration_checkpoint ())
8879
89- for i in range (n - 1 ):
90-
91- rhs [i + 1 ] = (1 / dt ) * ((u [i + 2 ,2 ] - u [i ,2 ]) / 2 * h )
92-
93- # solve the PPE using a SOR solver
94-
95- tolerance = 0.001
96- error = 1
97- omega = 1.7
98- max_iter = 1000
99- iter = 0
100-
101- while error >= tolerance :
102-
103- p [0 ] = p [1 ] # renew neumann pressure inlet
104-
105- for i in range (n - 1 ):
106- p [i + 1 ] = (1 - omega ) * p [i + 1 ] + ((omega * h ** 2 ) / 2 ) * (((p [i ] + p [i + 2 ]) / h ** 2 ) - rhs [i + 1 ])
107-
108- sum = 0
109- for i in range (n - 1 ):
110- val = ((p [i ] - 2 * p [i + 1 ] + p [i + 2 ]) / h ** 2 ) - rhs [i + 1 ]
111- sum += val * val
112-
113- error = np .sqrt (sum / n )
114-
115- iter += 1
116- if iter >= max_iter :
117- print ("SOR solver did not converge.\n " )
118- break
119-
120-
121- # calculate new velocities
122-
123- for i in range (n - 1 ):
124- u [i + 1 ,2 ] = u [i + 1 ,2 ] - dt * ((p [i + 2 ] - p [i + 1 ]) / h )
125-
126-
127- # write velocity to 3D solver
128-
129- write_vel = u [- 2 ,:]
130-
131- if interface .is_write_data_required (dt ): # write new velocities to 3D solver
132- interface .write_vector_data (
133- velocity_id , vertex_ids , write_vel )
134-
135- # transform data and write output data to vtk files
136-
137- u_print = np .reshape (u [:,2 ], n + 1 )
138- p_print = np .reshape (p , n + 1 )
139-
140- u_print = np .ascontiguousarray (u_print , dtype = np .float32 )
141- p_print = np .ascontiguousarray (p_print , dtype = np .float32 )
142-
143- filename = "./results/Fluid1D_" + str (counter )
144-
145- # TODO: uncomment if pyEVTK is installed and vtk output of 1D participant is wanted
146- # pointsToVTK(filename, x, y, z, data = {"U" : u_print, "p" : p_print})
147-
148- # advance simulation time
149-
150- dt = interface .advance (dt )
151- t = t + dt
152- counter += 1
80+ write_vel = [0 , 0 , u [- 1 ]]
81+ if interface .is_write_data_required (precice_dt ): # write new velocities to 3D solver
82+ interface .write_vector_data (velocity_id , vertex_ids , write_vel )
83+ precice_dt = interface .advance (precice_dt )
15384
15485 if interface .is_action_required (
155- precice .action_read_iteration_checkpoint ()):
156-
86+ precice .action_read_iteration_checkpoint ()):
15787 u = u_iter
15888 p = p_iter
159- t = t_iter
160- counter = counter_iter
161-
16289 interface .mark_action_fulfilled (
163- precice .action_read_iteration_checkpoint ())
90+ precice .action_read_iteration_checkpoint ())
16491
16592 interface .finalize ()
93+ return state0
94+
95+
96+ def postprocess (domain , ns , dt , ** arguments ):
97+
98+ bezier = domain .sample ('bezier' , 9 )
99+ x , u , p = bezier .eval (['x_i' , 'u_i' , 'p' ] @ ns , ** arguments )
100+
101+ ax1 = plt .subplot (211 )
102+ ax1 .plot (x , u )
103+ ax1 .set_title ("velocity u" )
104+ ax2 = plt .subplot (212 )
105+ ax2 .plot (x , p )
106+ ax2 .set_title ("pressure p" )
107+ plt .savefig ("./results/Fluid1D_" + str (dt ) + ".png" )
108+ return x , u , p
166109
167110
168- if __name__ == " __main__"
169- main ( )
111+ if __name__ == ' __main__'
112+ cli . run ( main )
0 commit comments