Skip to content

Commit 727ae8b

Browse files
RemDelaporteMathurinRemDelaporteMathurin
committed
move OF reader to its own file + MWE for multi-material
1 parent 7b7dabd commit 727ae8b

File tree

5 files changed

+221
-91
lines changed

5 files changed

+221
-91
lines changed

mwe_OF_reader.py renamed to OF_reader.py

Lines changed: 0 additions & 77 deletions
Original file line numberDiff line numberDiff line change
@@ -116,80 +116,3 @@ def find_closest_value(values: list[float], target: float) -> float:
116116
"""
117117
values = np.asarray(values) # Ensure input is a NumPy array
118118
return values[np.abs(values - target).argmin()]
119-
120-
121-
my_of_reader = OpenFOAMReader(filename=examples.download_cavity(load=False))
122-
vel = my_of_reader.create_dolfinx_function(t=2.5)
123-
124-
writer = dolfinx.io.VTXWriter(
125-
MPI.COMM_WORLD, "velocity_points_alt_test.bp", vel, engine="BP5"
126-
)
127-
writer.write(t=0)
128-
129-
130-
class TopSurface(F.SurfaceSubdomain):
131-
def locate_boundary_facet_indices(self, mesh):
132-
fdim = mesh.topology.dim - 1
133-
indices = locate_entities(
134-
mesh, fdim, lambda x: np.logical_and(np.isclose(x[1], 0.1), x[0] < 0.05)
135-
)
136-
return indices
137-
138-
139-
class BottomSurface(F.SurfaceSubdomain):
140-
def locate_boundary_facet_indices(self, mesh):
141-
fdim = mesh.topology.dim - 1
142-
indices = locate_entities(mesh, fdim, lambda x: np.isclose(x[1], 0))
143-
return indices
144-
145-
146-
my_model = F.HydrogenTransportProblem()
147-
148-
new_mesh = dolfinx.mesh.create_box(
149-
MPI.COMM_WORLD, [[0, 0, 0], [0.1, 0.1, 0.01]], n=[10, 10, 3]
150-
)
151-
152-
153-
# interpolate velocity on the new mesh
154-
element = basix.ufl.element("Lagrange", new_mesh.ufl_cell().cellname(), 1, shape=(3,))
155-
V = dolfinx.fem.functionspace(new_mesh, element)
156-
new_vel = dolfinx.fem.Function(V)
157-
F.helpers.nmm_interpolate(new_vel, vel)
158-
159-
160-
# my_model.mesh = F.Mesh(mesh=my_of_reader.dolfinx_mesh)
161-
my_model.mesh = F.Mesh(mesh=new_mesh)
162-
163-
my_mat = F.Material(name="mat", D_0=0.1, E_D=0)
164-
vol = F.VolumeSubdomain(id=1, material=my_mat)
165-
top = TopSurface(id=2)
166-
bot = BottomSurface(id=3)
167-
168-
my_model.subdomains = [vol, top, bot]
169-
170-
H = F.Species("H", mobile=True, subdomains=[vol])
171-
my_model.species = [H]
172-
173-
174-
my_model.temperature = 500
175-
176-
my_model.boundary_conditions = [
177-
F.FixedConcentrationBC(subdomain=top, value=1, species=H),
178-
F.FixedConcentrationBC(subdomain=bot, value=0, species=H),
179-
]
180-
181-
my_model.exports = [
182-
F.VTXSpeciesExport(filename="test_with_coupling_OF.bp", field=H),
183-
]
184-
185-
my_model.settings = F.Settings(
186-
atol=1e-10, rtol=1e-10, transient=True, stepsize=0.0001, final_time=0.03
187-
)
188-
189-
my_model.advection_terms = [
190-
F.AdvectionTerm(velocity=50 * new_vel, subdomain=vol, species=[H]),
191-
]
192-
193-
194-
my_model.initialise()
195-
my_model.run()

‎examples/advection_diffusion_multi_material.py‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,12 @@ def velocity_func(x):
9393

9494

9595
u = create_velocity_field()
96+
element = basix.ufl.element(
97+
"Lagrange", mesh.ufl_cell().cellname(), 1, shape=(mesh.geometry.dim,)
98+
)
99+
V = dolfinx.fem.functionspace(mesh, element)
100+
u_new = dolfinx.fem.Function(V)
101+
F.helpers.nmm_interpolate(u_new, u)
96102

97103
writer = dolfinx.io.VTXWriter(MPI.COMM_WORLD, "velocity.bp", u, engine="BP5")
98104
writer.write(t=0)
@@ -147,7 +153,7 @@ def velocity_func(x):
147153
]
148154

149155
my_model.advection_terms = [
150-
F.AdvectionTerm(velocity=u, species=[H], subdomain=bottom_domain)
156+
F.AdvectionTerm(velocity=u_new, species=[H], subdomain=bottom_domain)
151157
]
152158

153159

‎mwe_lid_driven_cavity.py‎

Lines changed: 91 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,91 @@
1+
from mpi4py import MPI
2+
3+
import basix
4+
import dolfinx
5+
import numpy as np
6+
import pyvista
7+
import ufl
8+
from dolfinx import mesh as dmesh
9+
from dolfinx.mesh import locate_entities
10+
from pyvista import examples
11+
12+
import festim as F
13+
14+
from OF_reader import OpenFOAMReader
15+
16+
17+
my_of_reader = OpenFOAMReader(filename=examples.download_cavity(load=False))
18+
vel = my_of_reader.create_dolfinx_function(t=2.5)
19+
20+
writer = dolfinx.io.VTXWriter(
21+
MPI.COMM_WORLD, "velocity_points_alt_test.bp", vel, engine="BP5"
22+
)
23+
writer.write(t=0)
24+
25+
26+
class TopSurface(F.SurfaceSubdomain):
27+
def locate_boundary_facet_indices(self, mesh):
28+
fdim = mesh.topology.dim - 1
29+
indices = locate_entities(
30+
mesh, fdim, lambda x: np.logical_and(np.isclose(x[1], 0.1), x[0] < 0.05)
31+
)
32+
return indices
33+
34+
35+
class BottomSurface(F.SurfaceSubdomain):
36+
def locate_boundary_facet_indices(self, mesh):
37+
fdim = mesh.topology.dim - 1
38+
indices = locate_entities(mesh, fdim, lambda x: np.isclose(x[1], 0))
39+
return indices
40+
41+
42+
my_model = F.HydrogenTransportProblem()
43+
44+
new_mesh = dolfinx.mesh.create_box(
45+
MPI.COMM_WORLD, [[0, 0, 0], [0.1, 0.1, 0.01]], n=[10, 10, 3]
46+
)
47+
48+
49+
# interpolate velocity on the new mesh
50+
element = basix.ufl.element("Lagrange", new_mesh.ufl_cell().cellname(), 1, shape=(3,))
51+
V = dolfinx.fem.functionspace(new_mesh, element)
52+
new_vel = dolfinx.fem.Function(V)
53+
F.helpers.nmm_interpolate(new_vel, vel)
54+
55+
56+
# my_model.mesh = F.Mesh(mesh=my_of_reader.dolfinx_mesh)
57+
my_model.mesh = F.Mesh(mesh=new_mesh)
58+
59+
my_mat = F.Material(name="mat", D_0=0.1, E_D=0)
60+
vol = F.VolumeSubdomain(id=1, material=my_mat)
61+
top = TopSurface(id=2)
62+
bot = BottomSurface(id=3)
63+
64+
my_model.subdomains = [vol, top, bot]
65+
66+
H = F.Species("H", mobile=True, subdomains=[vol])
67+
my_model.species = [H]
68+
69+
70+
my_model.temperature = 500
71+
72+
my_model.boundary_conditions = [
73+
F.FixedConcentrationBC(subdomain=top, value=1, species=H),
74+
F.FixedConcentrationBC(subdomain=bot, value=0, species=H),
75+
]
76+
77+
my_model.exports = [
78+
F.VTXSpeciesExport(filename="test_with_coupling_OF.bp", field=H),
79+
]
80+
81+
my_model.settings = F.Settings(
82+
atol=1e-10, rtol=1e-10, transient=True, stepsize=0.0001, final_time=0.03
83+
)
84+
85+
my_model.advection_terms = [
86+
F.AdvectionTerm(velocity=50 * new_vel, subdomain=vol, species=[H]),
87+
]
88+
89+
90+
my_model.initialise()
91+
my_model.run()

‎mwe_lid_driven_cavity_fluid_solid.py‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
from mpi4py import MPI
2+
3+
import basix
4+
import dolfinx
5+
import numpy as np
6+
from dolfinx.mesh import locate_entities
7+
from pyvista import examples
8+
9+
import festim as F
10+
11+
from OF_reader import OpenFOAMReader
12+
13+
14+
my_of_reader = OpenFOAMReader(filename=examples.download_cavity(load=False))
15+
vel = my_of_reader.create_dolfinx_function(t=2.5)
16+
17+
writer = dolfinx.io.VTXWriter(
18+
MPI.COMM_WORLD, "velocity_points_alt_test.bp", vel, engine="BP5"
19+
)
20+
writer.write(t=0)
21+
22+
23+
class TopSurface(F.SurfaceSubdomain):
24+
def locate_boundary_facet_indices(self, mesh):
25+
fdim = mesh.topology.dim - 1
26+
indices = locate_entities(
27+
mesh, fdim, lambda x: np.logical_and(np.isclose(x[1], 0.1), x[0] < 0.05)
28+
)
29+
return indices
30+
31+
32+
class BottomSurface(F.SurfaceSubdomain):
33+
def locate_boundary_facet_indices(self, mesh):
34+
fdim = mesh.topology.dim - 1
35+
indices = locate_entities(mesh, fdim, lambda x: np.isclose(x[1], 0))
36+
return indices
37+
38+
39+
class FluiDomain(F.VolumeSubdomain):
40+
def locate_subdomain_entities(self, mesh):
41+
indices = locate_entities(
42+
mesh,
43+
3,
44+
lambda x: x[0] <= of_x + 1e-16,
45+
)
46+
return indices
47+
48+
49+
class SolidDomain(F.VolumeSubdomain):
50+
def locate_subdomain_entities(self, mesh):
51+
indices = locate_entities(
52+
mesh,
53+
3,
54+
lambda x: x[0] > of_x - 1e-16,
55+
)
56+
return indices
57+
58+
59+
my_model = F.HTransportProblemDiscontinuous()
60+
61+
of_x, of_y, of_z = 0.1, 0.1, 0.01
62+
solid_x = 0.01
63+
64+
new_mesh = dolfinx.mesh.create_box(
65+
MPI.COMM_WORLD, [[0, 0, 0], [of_x + solid_x, of_y, of_z]], n=[11, 10, 3]
66+
)
67+
68+
69+
# interpolate velocity on the new mesh
70+
element = basix.ufl.element("Lagrange", new_mesh.ufl_cell().cellname(), 1, shape=(3,))
71+
V = dolfinx.fem.functionspace(new_mesh, element)
72+
new_vel = dolfinx.fem.Function(V)
73+
F.helpers.nmm_interpolate(new_vel, vel)
74+
75+
76+
# I want to time the values by 50
77+
new_vel.x.array[:] *= 50
78+
79+
# my_model.mesh = F.Mesh(mesh=my_of_reader.dolfinx_mesh)
80+
my_model.mesh = F.Mesh(mesh=new_mesh)
81+
82+
mat_fluid = F.Material(D_0=0.1, E_D=0, K_S_0=1, E_K_S=0)
83+
mat_solid = F.Material(D_0=0.1, E_D=0, K_S_0=3, E_K_S=0)
84+
vol_fluid = FluiDomain(id=1, material=mat_fluid)
85+
vol_solid = SolidDomain(id=2, material=mat_solid)
86+
top = TopSurface(id=2)
87+
bot = BottomSurface(id=3)
88+
89+
my_model.subdomains = [vol_fluid, vol_solid, top, bot]
90+
my_model.interfaces = [F.Interface(4, (vol_fluid, vol_solid), penalty_term=30)]
91+
my_model.surface_to_volume = {top: vol_fluid, bot: vol_fluid}
92+
H = F.Species("H", mobile=True, subdomains=[vol_fluid, vol_solid])
93+
my_model.species = [H]
94+
95+
96+
my_model.temperature = 500
97+
98+
my_model.boundary_conditions = [
99+
F.FixedConcentrationBC(subdomain=top, value=1, species=H),
100+
F.FixedConcentrationBC(subdomain=bot, value=0, species=H),
101+
]
102+
103+
my_model.exports = [
104+
F.VTXSpeciesExport(
105+
filename="test_with_coupling_OF_fluid.bp", field=H, subdomain=vol_fluid
106+
),
107+
F.VTXSpeciesExport(
108+
filename="test_with_coupling_OF_solid.bp", field=H, subdomain=vol_solid
109+
),
110+
]
111+
112+
my_model.settings = F.Settings(
113+
atol=1e-10, rtol=1e-10, transient=True, stepsize=0.0001, final_time=0.03
114+
)
115+
116+
my_model.advection_terms = [
117+
F.AdvectionTerm(velocity=new_vel, subdomain=vol_fluid, species=[H]),
118+
]
119+
120+
121+
my_model.initialise()
122+
my_model.run()

‎src/festim/hydrogen_transport_problem.py‎

Lines changed: 1 addition & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1169,19 +1169,7 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
11691169
conc = spe.subdomain_to_solution[subdomain]
11701170
vel = adv_term.velocity
11711171

1172-
# since the meshes are different, we need to interpolate onto a new functionspace
1173-
submesh = subdomain.submesh
1174-
v_cg = basix.ufl.element(
1175-
"Lagrange",
1176-
submesh.topology.cell_name(),
1177-
2,
1178-
shape=(submesh.geometry.dim,),
1179-
)
1180-
V_velocity = dolfinx.fem.functionspace(submesh, v_cg)
1181-
vel_2 = dolfinx.fem.Function(V_velocity)
1182-
nmm_interpolate(vel_2, vel)
1183-
1184-
form += ufl.inner(ufl.dot(ufl.grad(conc), vel_2), v) * self.dx(
1172+
form += ufl.inner(ufl.dot(ufl.grad(conc), vel), v) * self.dx(
11851173
subdomain.id
11861174
)
11871175

0 commit comments

Comments
 (0)