Created
July 30, 2025 15:01
-
-
Save RemDelaporteMathurin/c4a2bc4794ca85b0391cdbb865b9b7ac to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from mpi4py import MPI | |
| import dolfinx | |
| from dolfinx import log | |
| from dolfinx.fem.petsc import NonlinearProblem | |
| from dolfinx.nls.petsc import NewtonSolver | |
| import numpy as np | |
| from scifem import create_real_functionspace | |
| import ufl | |
| M = 20 | |
| mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, M, dtype=np.float64) | |
| V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1)) | |
| R = create_real_functionspace(mesh) | |
| W = ufl.MixedFunctionSpace(V, R) | |
| tdim = mesh.topology.dim | |
| fdim = tdim - 1 | |
| mesh.topology.create_connectivity(fdim, tdim) | |
| def boundary_right(x): | |
| return np.isclose(x[0], 1) | |
| def other_boundaries(x): | |
| return np.logical_or(np.isclose(x[0], 0), np.isclose(x[1], 0), np.isclose(x[1], 1)) | |
| dofs_right = dolfinx.fem.locate_dofs_geometrical(V, boundary_right) | |
| dofs_other = dolfinx.fem.locate_dofs_geometrical(V, other_boundaries) | |
| u = dolfinx.fem.Function(V) | |
| u_n = dolfinx.fem.Function(V) | |
| dt = dolfinx.fem.Constant(mesh, 1.0) | |
| K_S = dolfinx.fem.Constant(mesh, 2.0) | |
| e = 2.0 | |
| volume = 4.0 | |
| v, pressure_v = ufl.TestFunctions(W) | |
| pressure = dolfinx.fem.Function(R) | |
| pressure_n = dolfinx.fem.Function(R) | |
| pressure_ini_expr = dolfinx.fem.Expression( | |
| dolfinx.fem.Constant(mesh, 3.5), R.element.interpolation_points() | |
| ) | |
| pressure_n.interpolate(pressure_ini_expr) | |
| pressure.x.array[:] = pressure_n.x.array[:] | |
| bc_right_expr = dolfinx.fem.Expression( | |
| K_S * pressure_n**0.5, V.element.interpolation_points() | |
| ) | |
| u_bc_right = dolfinx.fem.Function(V) | |
| u_bc_right.interpolate(bc_right_expr) | |
| bc_right = dolfinx.fem.dirichletbc(u_bc_right, dofs_right) | |
| bc_other = dolfinx.fem.dirichletbc(dolfinx.fem.Constant(mesh, 0.0), dofs_other, V) | |
| n = ufl.FacetNormal(mesh) | |
| flux = ufl.inner(ufl.grad(u), n) * ufl.ds() | |
| F = -(u - u_n) / dt * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx | |
| F += (pressure - pressure_n) / dt * pressure_v * ufl.dx - e / volume * flux * ufl.dx | |
| problem = NonlinearProblem(F, u, bcs=[bc_right, bc_other]) | |
| solver = NewtonSolver(MPI.COMM_WORLD, problem) | |
| log.set_log_level(log.LogLevel.INFO) | |
| t = 0 | |
| while t < 10: | |
| t += dt.value | |
| print(f"Solving at time {t:.2f}") | |
| # Solve the problem | |
| n, converged = solver.solve(u) | |
| u_n.x.array[:] = u.x.array[:] |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment