From 39bd6eea22b32f3b718b7c5791163ed717c8eef7 Mon Sep 17 00:00:00 2001 From: brookeslawski Date: Tue, 20 Feb 2024 16:41:17 -0700 Subject: [PATCH] black reformatting --- examples/rayleigh-benard-chorin-dolfinx.py | 242 +++++++++++++-------- 1 file changed, 152 insertions(+), 90 deletions(-) diff --git a/examples/rayleigh-benard-chorin-dolfinx.py b/examples/rayleigh-benard-chorin-dolfinx.py index a9f842e7..6bcc95a3 100644 --- a/examples/rayleigh-benard-chorin-dolfinx.py +++ b/examples/rayleigh-benard-chorin-dolfinx.py @@ -12,13 +12,45 @@ import gmsh from dolfinx import io -from dolfinx.fem import Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical, locate_dofs_topological -from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc +from dolfinx.fem import ( + Constant, + Function, + FunctionSpace, + assemble_scalar, + dirichletbc, + form, + locate_dofs_geometrical, + locate_dofs_topological, +) +from dolfinx.fem.petsc import ( + assemble_matrix, + assemble_vector, + apply_lifting, + create_vector, + set_bc, +) + # from dolfinx.io import VTXWriter from dolfinx.mesh import create_rectangle, CellType, locate_entities_boundary + # from dolfinx.plot import vtk_mesh -from ufl import (FacetNormal, FiniteElement, Identity, TestFunction, TrialFunction, VectorElement, - div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym) +from ufl import ( + FacetNormal, + FiniteElement, + Identity, + TestFunction, + TrialFunction, + VectorElement, + div, + dot, + ds, + dx, + inner, + lhs, + nabla_grad, + rhs, + sym, +) # ================================================================ # Inputs @@ -31,16 +63,16 @@ y_max = 1.0 h = 0.05 -nx = 50 # 150 # int((x_max - x_min)/h) -ny = 10 # 50 # int((y_max - y_min)/h) +nx = 50 # 150 # int((x_max - x_min)/h) +ny = 10 # 50 # int((y_max - y_min)/h) -pv_panel_present = True # empty domain or with a pv panel in the center? +pv_panel_present = True # empty domain or with a pv panel in the center? # ================================================================ # Build Mesh # ================================================================ -if pv_panel_present: +if pv_panel_present: comm = MPI.COMM_WORLD gmsh.initialize() @@ -52,25 +84,28 @@ ndim = 2 - domain_width = 3.0 # box width - domain_height = 1.0 # # box height + domain_width = 3.0 # box width + domain_height = 1.0 # # box height - domain_id = gmsh_model.occ.addRectangle(0, 0, 0, domain_width, domain_height) # Notice this spans from [0, x_max], [0, y_max], your BCs may need adjustment + domain_id = gmsh_model.occ.addRectangle( + 0, 0, 0, domain_width, domain_height + ) # Notice this spans from [0, x_max], [0, y_max], your BCs may need adjustment domain_tag = (ndim, domain_id) - panel_width = 0.5 # Chord length, or width - panel_height = 0.05 # Sets the panel thickness, really - panel_angle = np.radians(30) # Sets the panel rotation (argument must be radians for gmsh) + panel_width = 0.5 # Chord length, or width + panel_height = 0.05 # Sets the panel thickness, really + panel_angle = np.radians( + 30 + ) # Sets the panel rotation (argument must be radians for gmsh) - panel_id = gmsh_model.occ.addRectangle(-0.5*panel_width, -0.5*panel_height, 0, panel_width, panel_height) + panel_id = gmsh_model.occ.addRectangle( + -0.5 * panel_width, -0.5 * panel_height, 0, panel_width, panel_height + ) panel_tag = (ndim, panel_id) # Rotate the panel and shift it into its correct position gmsh_model.occ.rotate([panel_tag], 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, panel_angle) - gmsh_model.occ.translate([panel_tag], - 0.5*domain_width, - 0.5*domain_height, - 0.0) + gmsh_model.occ.translate([panel_tag], 0.5 * domain_width, 0.5 * domain_height, 0.0) # Cookie cutter step, domain = domain - panel, is how to read this gmsh_model.occ.cut([domain_tag], [panel_tag]) @@ -79,7 +114,7 @@ all_pts = gmsh_model.occ.getEntities(0) - l_characteristic = 0.05 # Sets the characteristic size of the cells + l_characteristic = 0.05 # Sets the characteristic size of the cells gmsh_model.mesh.setSize(all_pts, l_characteristic) vol_tag_list = gmsh_model.occ.getEntities(ndim) @@ -95,8 +130,12 @@ else: # create an empty domain mesh - mesh = create_rectangle(MPI.COMM_WORLD, [np.array([x_min, y_min]), np.array([x_max, y_max])], - [nx, ny], CellType.triangle) + mesh = create_rectangle( + MPI.COMM_WORLD, + [np.array([x_min, y_min]), np.array([x_max, y_max])], + [nx, ny], + CellType.triangle, + ) # Two key physical parameters are the Rayleigh number (Ra), which @@ -129,9 +168,9 @@ v_cg2 = VectorElement("Lagrange", mesh.ufl_cell(), 2) q_cg1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1) -V = FunctionSpace(mesh, v_cg2) # velocity -Q = FunctionSpace(mesh, q_cg1) # pressure -S = FunctionSpace(mesh, q_cg1) # temperature +V = FunctionSpace(mesh, v_cg2) # velocity +Q = FunctionSpace(mesh, q_cg1) # pressure +S = FunctionSpace(mesh, q_cg1) # temperature # velocity u = TrialFunction(V) @@ -150,33 +189,39 @@ # temperature theta = TrialFunction(S) s = TestFunction(S) -T_n = Function(S) # for outputting T, calculated from theta for each timestep +T_n = Function(S) # for outputting T, calculated from theta for each timestep T_n.name = "T_n" theta_n = Function(S) -theta_ = Function(S) +theta_ = Function(S) -#%% ================================================================ +# %% ================================================================ # Build Boundary Conditions # ================================================================ -def left_wall( x): + +def left_wall(x): return np.isclose(x[0], x_min) -def right_wall( x): + +def right_wall(x): return np.isclose(x[0], x_max) -def bottom_wall( x): + +def bottom_wall(x): return np.isclose(x[1], y_min) -def top_wall( x): + +def top_wall(x): return np.isclose(x[1], y_max) -def internal_boundaries( x): + +def internal_boundaries(x): tol = 1e-3 x_test = np.logical_and(x_min + tol < x[0], x[0] < x_max - tol) y_test = np.logical_and(y_min + tol < x[1], x[1] < y_max - tol) return np.logical_and(x_test, y_test) + # Velocity Boundary Conditions left_wall_dofs = locate_dofs_geometrical(V, left_wall) u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType) @@ -193,9 +238,10 @@ def internal_boundaries( x): bcu = [bcu_left_wall, bcu_right_wall, bcu_bottom_wall, bcu_top_wall] -if pv_panel_present: - - boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, internal_boundaries) +if pv_panel_present: + boundary_facets = locate_entities_boundary( + mesh, mesh.topology.dim - 1, internal_boundaries + ) boundary_dofs = locate_dofs_topological(V, mesh.topology.dim - 1, boundary_facets) bcu_internal_walls = dirichletbc(u_noslip, boundary_dofs, V) @@ -213,46 +259,59 @@ def internal_boundaries( x): # temperature boundary conditions T0_top = 0 T0_bottom = 1 -T0_val = 0.5 # 0.5 #10 # TODO - change this to expression -DeltaT = T0_bottom-T0_top # ? should this be defined as Constant(mesh, PETSc.ScalarType(bottom-top)) ? +T0_val = 0.5 # 0.5 #10 # TODO - change this to expression +DeltaT = ( + T0_bottom - T0_top +) # ? should this be defined as Constant(mesh, PETSc.ScalarType(bottom-top)) ? # Interpolate initial condition # T0_int = Function(S) -T_n.interpolate(lambda x: (T0_bottom + (x[1]/y_max)*(T0_top-T0_bottom))) +T_n.interpolate(lambda x: (T0_bottom + (x[1] / y_max) * (T0_top - T0_bottom))) # theta0_int = Function(S) -theta_n.interpolate(lambda x: (-x[1]/y_max)) +theta_n.interpolate(lambda x: (-x[1] / y_max)) -#initialize T_n +# initialize T_n # T_n.x.array[:] = DeltaT*theta_n.x.array[:] + T0_bottom # is this necessary? # keeping these separate in case we want to specify different temperatures for when pv panels are present -if pv_panel_present==False: - - print('applying bottom wall temp = {}'.format((T0_bottom-T0_bottom)/DeltaT)) +if pv_panel_present == False: + print("applying bottom wall temp = {}".format((T0_bottom - T0_bottom) / DeltaT)) bottom_wall_dofs = locate_dofs_geometrical(S, bottom_wall) - bcT_bottom_wall = dirichletbc(PETSc.ScalarType((T0_bottom-T0_bottom)/DeltaT), bottom_wall_dofs, S) + bcT_bottom_wall = dirichletbc( + PETSc.ScalarType((T0_bottom - T0_bottom) / DeltaT), bottom_wall_dofs, S + ) - print('applying top wall temp = {}'.format((T0_top-T0_bottom)/DeltaT)) + print("applying top wall temp = {}".format((T0_top - T0_bottom) / DeltaT)) top_wall_dofs = locate_dofs_geometrical(S, top_wall) - bcT_top_wall = dirichletbc(PETSc.ScalarType((T0_top-T0_bottom)/DeltaT), top_wall_dofs, S) + bcT_top_wall = dirichletbc( + PETSc.ScalarType((T0_top - T0_bottom) / DeltaT), top_wall_dofs, S + ) bcT = [bcT_top_wall, bcT_bottom_wall] -if pv_panel_present: +if pv_panel_present: # apply temperature of 0 at all walls except PV panel which is 1 - print('applying bottom wall temp = {}'.format((T0_bottom-T0_bottom)/DeltaT)) + print("applying bottom wall temp = {}".format((T0_bottom - T0_bottom) / DeltaT)) bottom_wall_dofs = locate_dofs_geometrical(S, bottom_wall) - bcT_bottom_wall = dirichletbc(PETSc.ScalarType((T0_bottom-T0_bottom)/DeltaT), bottom_wall_dofs, S) + bcT_bottom_wall = dirichletbc( + PETSc.ScalarType((T0_bottom - T0_bottom) / DeltaT), bottom_wall_dofs, S + ) - print('applying top wall temp = {}'.format((T0_top-T0_bottom)/DeltaT)) + print("applying top wall temp = {}".format((T0_top - T0_bottom) / DeltaT)) top_wall_dofs = locate_dofs_geometrical(S, top_wall) - bcT_top_wall = dirichletbc(PETSc.ScalarType((T0_top-T0_bottom)/DeltaT), top_wall_dofs, S) - - print('applying internal boundary temp') - boundary_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, internal_boundaries) + bcT_top_wall = dirichletbc( + PETSc.ScalarType((T0_top - T0_bottom) / DeltaT), top_wall_dofs, S + ) + + print("applying internal boundary temp") + boundary_facets = locate_entities_boundary( + mesh, mesh.topology.dim - 1, internal_boundaries + ) boundary_dofs = locate_dofs_topological(S, mesh.topology.dim - 1, boundary_facets) - bcT_internal_walls = dirichletbc(PETSc.ScalarType((T0_bottom-T0_bottom)/DeltaT), boundary_dofs, S) + bcT_internal_walls = dirichletbc( + PETSc.ScalarType((T0_bottom - T0_bottom) / DeltaT), boundary_dofs, S + ) bcT = [bcT_top_wall, bcT_bottom_wall, bcT_internal_walls] @@ -265,7 +324,7 @@ def internal_boundaries( x): # bcp_left_wall = dirichletbc(PETSc.ScalarType(pressure_bc), left_wall_dofs, Q) # bcp_right_wall = dirichletbc(PETSc.ScalarType(pressure_bc), right_wall_dofs, Q) -bcp = [] # [bcp_left_wall, bcp_right_wall, bcp_bottom_wall, bcp_top_wall] +bcp = [] # [bcp_left_wall, bcp_right_wall, bcp_bottom_wall, bcp_top_wall] # ================================================================ # Build All Forms @@ -273,39 +332,42 @@ def internal_boundaries( x): # step 1: tentative velocity # chorin (removed the pressure term) -F1 = (1 / Pr) * ((1 / dt) * inner(u - u_n, v) * dx - + inner(nabla_grad(u_n) * u_n, v) * dx) # this might be dot not * ? -F1 += nu * inner(nabla_grad(u), nabla_grad(v)) * dx -F1 -= Ra*inner(theta_n*g,v)*dx +F1 = (1 / Pr) * ( + (1 / dt) * inner(u - u_n, v) * dx + inner(nabla_grad(u_n) * u_n, v) * dx +) # this might be dot not * ? +F1 += nu * inner(nabla_grad(u), nabla_grad(v)) * dx +F1 -= Ra * inner(theta_n * g, v) * dx -a1 = form(lhs(F1)) # dependent on u +a1 = form(lhs(F1)) # dependent on u L1 = form(rhs(F1)) # step 2: pressure correction -a2 = form(inner(nabla_grad(p), nabla_grad(q))*dx) -L2 = form(-(1.0/dt)*div(u_)*q*dx) # needs to be reassembled +a2 = form(inner(nabla_grad(p), nabla_grad(q)) * dx) +L2 = form(-(1.0 / dt) * div(u_) * q * dx) # needs to be reassembled # step 3: velocity update -a3 = form(inner(u, v)*dx) # doesn't need to be reassembled -L3 = form(inner(u_, v)*dx - dt*inner(nabla_grad(p_), v)*dx) +a3 = form(inner(u, v) * dx) # doesn't need to be reassembled +L3 = form(inner(u_, v) * dx - dt * inner(nabla_grad(p_), v) * dx) # step 4: temperature update? -a4 = form((1/dt)*inner((theta), s)*dx - + inner(nabla_grad(theta), nabla_grad(s))*dx - + inner(dot(u_, nabla_grad(theta)), s)*dx) # needs to be reassembled bc of u_ -L4 = form((1/dt)*inner(theta_n, s)*dx) # needs to be reassembled bc of theta_n +a4 = form( + (1 / dt) * inner((theta), s) * dx + + inner(nabla_grad(theta), nabla_grad(s)) * dx + + inner(dot(u_, nabla_grad(theta)), s) * dx +) # needs to be reassembled bc of u_ +L4 = form((1 / dt) * inner(theta_n, s) * dx) # needs to be reassembled bc of theta_n # Solver for step 1 solver1 = PETSc.KSP().create(mesh.comm) -solver1.setType(PETSc.KSP.Type.GMRES) # TODO - test solution with BCGS +solver1.setType(PETSc.KSP.Type.GMRES) # TODO - test solution with BCGS pc1 = solver1.getPC() pc1.setType(PETSc.PC.Type.HYPRE) pc1.setHYPREType("boomeramg") # Solver for step 2 solver2 = PETSc.KSP().create(mesh.comm) -solver2.setType(PETSc.KSP.Type.GMRES) # TODO - test solution with BCGS +solver2.setType(PETSc.KSP.Type.GMRES) # TODO - test solution with BCGS pc2 = solver2.getPC() pc2.setType(PETSc.PC.Type.HYPRE) # pc2.setHYPREType("boomeramg") # TODO - test solution with this instead (for speed?) @@ -314,7 +376,7 @@ def internal_boundaries( x): solver3 = PETSc.KSP().create(mesh.comm) solver3.setType(PETSc.KSP.Type.GMRES) pc3 = solver3.getPC() -pc3.setType(PETSc.PC.Type.JACOBI) # TODO - test solution with SOR +pc3.setType(PETSc.PC.Type.JACOBI) # TODO - test solution with SOR # Solver for step 2 solver4 = PETSc.KSP().create(mesh.comm) @@ -328,11 +390,11 @@ def internal_boundaries( x): # ================================================================ eps = 3.0e-16 -t = dt_num #dt # 0.0 -ct = 1 #0 -save_interval = 1 #50 +t = dt_num # dt # 0.0 +ct = 1 # 0 +save_interval = 1 # 50 -t_final = 0.1 #0.5 # 0.5 #0.1 # 0.000075 +t_final = 0.1 # 0.5 # 0.5 #0.1 # 0.000075 if pv_panel_present: save_fn = "rayleigh-benard-pv.xdmf" @@ -340,7 +402,6 @@ def internal_boundaries( x): save_fn = "rayleigh-benard.xdmf" with io.XDMFFile(mesh.comm, save_fn, "w") as xdmf: - xdmf.write_mesh(mesh) xdmf.write_function(u_n, 0) xdmf.write_function(p_n, 0) @@ -364,11 +425,9 @@ def internal_boundaries( x): b4 = assemble_vector(L4) while t < t_final + eps: - - T_n.x.array[:] = DeltaT*theta_n.x.array[:] + T0_bottom + T_n.x.array[:] = DeltaT * theta_n.x.array[:] + T0_bottom with io.XDMFFile(mesh.comm, save_fn, "a") as xdmf: - xdmf.write_function(u_n, t) xdmf.write_function(p_n, t) xdmf.write_function(T_n, t) @@ -378,8 +437,8 @@ def internal_boundaries( x): # Assemble and Build Solvers # ================================================================ - A1.zeroEntries() # resets the matrix - A1 = assemble_matrix(A1,a1,bcs=bcu) + A1.zeroEntries() # resets the matrix + A1 = assemble_matrix(A1, a1, bcs=bcu) A1.assemble() solver1.setOperators(A1) @@ -394,7 +453,7 @@ def internal_boundaries( x): # Step 1: Tentative velocity solve with b1.localForm() as loc_1: loc_1.set(0) - b1=assemble_vector(b1, L1) + b1 = assemble_vector(b1, L1) apply_lifting(b1, [a1], [bcu]) b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) set_bc(b1, bcu) @@ -404,7 +463,7 @@ def internal_boundaries( x): # Step 2: Pressure corrrection step with b2.localForm() as loc_2: loc_2.set(0) - b2=assemble_vector(b2, L2) + b2 = assemble_vector(b2, L2) apply_lifting(b2, [a2], [bcp]) b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) set_bc(b2, bcp) @@ -414,7 +473,7 @@ def internal_boundaries( x): # Step 3: Velocity correction step with b3.localForm() as loc_3: loc_3.set(0) - b3=assemble_vector(b3, L3) + b3 = assemble_vector(b3, L3) apply_lifting(b3, [a3], [bcu]) b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) solver3.solve(b3, u_.vector) @@ -423,13 +482,13 @@ def internal_boundaries( x): # Step 4: Temperature corrrection step with b4.localForm() as loc_4: loc_4.set(0) - b4=assemble_vector(b4, L4) + b4 = assemble_vector(b4, L4) apply_lifting(b4, [a4], [bcT]) b4.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) set_bc(b4, bcT) solver4.solve(b4, theta_.vector) theta_.x.scatter_forward() - + # Update variable with solution form this time step u_n.x.array[:] = u_.x.array[:] p_n.x.array[:] = p_.x.array[:] @@ -442,7 +501,10 @@ def internal_boundaries( x): T_n_sum = mesh.comm.allreduce(np.sum(T_n.vector.array), op=MPI.SUM) if ct % save_interval == 0: - print('Time = %.6f, u_max = %.6e, p_max = %.6e, T_max = %.6e, T_sum = %.6e' % (t, u_n_max, p_n_max, T_n_max, T_n_sum)) + print( + "Time = %.6f, u_max = %.6e, p_max = %.6e, T_max = %.6e, T_sum = %.6e" + % (t, u_n_max, p_n_max, T_n_max, T_n_sum) + ) # Move to next step t += float(dt) @@ -462,4 +524,4 @@ def internal_boundaries( x): # print('size1 = ',np.shape(coords_better[0, :])) # print('size2 = ',np.shape(coords_better[0, :])) # plt.scatter(coords_better[:, 0], coords_better[:, 1], c=np.sqrt(u_k.vector.array[0::2]**2 + u_k.vector.array[1::2]**2)) -# plt.show() \ No newline at end of file +# plt.show()