Dear all,
I am trying to solve a transient heat equation. I am completely newbie in FEniCS.
For inverse part, i wanted to find the value of the coefficients of thermal conductivity.
However, at first iteration |proj g|= 1.98120D-14 and program stoped.
Please guide me, how can i solve the issue.
My complete code is below:
#“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”“”"
from dolfin import *
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from fenics_adjoint import *
rho = Constant(1.46555e6)
C= Constant(20.0)
dt= 1
num_steps=10
mesh = RectangleMesh(Point(0,0),Point(1,1), 10, 10)
V = FunctionSpace(mesh, "Lagrange", 1)
x=SpatialCoordinate(mesh)
M = FunctionSpace(mesh,"CG",1)
Q_heat = 20000.0
Q = Constant(0.0)
g = Expression('Q_heat', degree=1, Q_heat=Constant(Q_heat))
#Define boudary conditions
u_0 = Constant(0.0) #Initial temperature over the boundary is zero.
u0=interpolate(u_0, V)
#Define Left Boundary
def left_boundary( x, on_boundary):
return on_boundary and (abs(x[0])< DOLFIN_EPS)
bc = DirichletBC(V, Constant(0.0), left_boundary)
class right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and ((abs(x[0] - 1.) < DOLFIN_EPS))
vertical = right()
tdim = mesh.topology().dim()
boundaries = MeshFunction('size_t', mesh, tdim-1)
boundaries.set_all(0)
vertical.mark(boundaries, 1)
dsB = Measure("ds", subdomain_id=1, subdomain_data=boundaries)
d_target=interpolate(Constant(20.0), M)
d_guess =interpolate(Constant(40.0), M)
e_target=interpolate(Constant(100.0), M)
e_guess=interpolate(Constant(100.0), M)
def forward(d, e):
u = TrialFunction(V)
v = TestFunction(V)
u0 = Function(V)
#Define Thermal conductivity
def k():
return d + e*x[1]
a= rho*C*u*v*dx+ dt*inner(grad(v),k()*grad(u))*dx
L= dt*g*v*dsB+rho*C*inner(u0,v)*dx+dt*Q*v*dx
u = Function(V)
#vtkfile = File('Heat_Transfer/solution.pvd')
t = 0
for n in range (num_steps):
t += dt
#Compute solution
solve(a==L, u, bc)
#Save solution
#vtkfile << (u, t)
plot(u)
#update solution
u0.assign(u)
return u,d,e
[u_target, d_target_1, e_target_1] = forward(d_target, e_target)
File("forward_modified/targetTempereature_modified.pvd") << u_target
File("forward_modified/target_d_modified.pvd") << d_target_1
File("forward_modified/target_e_up.pvd") << e_target_1
[u,d,e] = forward(d_guess, e_guess)
output = File("inverse_Up/inverse_result_up.pvd")
d_viz = Function(M, name="d")
def eval_cb(j, d):
d_viz.assign(d)
output << d_viz
alpha = Constant(0.0)
beta = Constant(1e-2)
#objective function
J = assemble((0.5*inner((u-u_target),u-u_target))*dx+0.5*alpha*sqrt(beta*beta+dot(grad(d),grad(d)))*dx)
#Optimize variables
d_c = Control(d)
J_rf = ReducedFunctional(J, d_c, eval_cb_post=eval_cb)
# scipy optimizer
d_opt = minimize(J_rf,method='L-BFGS-B', bounds=(5.0, 200.0), options={"gtol":1e-15,"ftol":1e-15})