Hello Everyone
I am rather new to fenics and am trying to solve a problem, where I want to achieve a given heat-distribution (z) within my domain by optimising a boundary-term (g). My problem is that my boundary-term g only appears in the surface-integrals (ds) and in none of the dx. I cheated a bit and put an extra volume-integral with a tiny constant into the bilinear form, which set g (not on boundary) to basically 0 and made the code run. Otherwise the two print() in the end of the code printed “inf inf”.
So I’m wondering what the best way is to solve this problem, without artificially setting the inner g = 0. I have tried to create subdomains and “remove” the inner DoFs of g, but have not gotten it to work yet. Below is the code with the reg_term that works.
I’m using dolfinx & basix 0.10.0. Thank you for your help!
"""
1) Mesh and parameters
"""
domain = mesh.create_rectangle(
MPI.COMM_WORLD,
[(0.0, 0.0),
(2.0, 1.0)],
[40, 20],
cell_type=mesh.CellType.triangle,
)
# Parameters
rho = Constant(domain, PETSc.ScalarType(50))
alpha = Constant(domain, PETSc.ScalarType(1e-3))
# Desired Temperature Distribution
x = ufl.SpatialCoordinate(domain)
z = (x[0] - 1) ** 2
"""
2) Mixed functionspace
"""
# Base-element
E = element("Lagrange", domain.basix_cell(), 1)
# Mixed element space for (T, g, phi)
ME = mixed_element([E, E, E])
W = functionspace(domain, ME)
# Unknowns and Testfunctons
T, g, phi = TrialFunctions(W)
t, h, psi = TestFunctions(W)
"""
3) Bilinear & Linear form
"""
a = (
T * t * dx + inner(grad(t), grad(phi)) * dx + rho * phi * t * ds
+ alpha * g * h * ds - rho * phi * h * ds
+ inner(grad(T), grad(psi)) * dx + rho * psi * (T - g) * ds
)
L = z * t * dx
"""
4) Solve linear problem
"""
reg_term = Constant(domain, PETSc.ScalarType(1e-8))
a += reg_term * g * h * dx
problem = LinearProblem(a, L, petsc_options_prefix="heat_", petsc_options={"ksp_type": "gmres", "pc_type": "ilu"})
uh = problem.solve()
"""
5) Results
"""
V_T, map_T = W.sub(0).collapse()
V_g, map_g = W.sub(1).collapse()
T_function = Function(V_T, name="T")
g_function = Function(V_g, name="g")
T_function.x.array[:] = uh.x.array[map_T]
g_function.x.array[:] = uh.x.array[map_g]
print("uh min/max:", float(np.min(uh.x.array)), float(np.max(uh.x.array)))
print("T min/max:", float(np.min(T_function.x.array)), float(np.max(T_function.x.array)))