Analytical conduction adjoint
The cost function is described as:
where,
are the unknowns, and in our case are prior estimates of unknowns are the uncertainty of these prior estimates is the model we are solving, in our case the heat equation are observables on the solution, surface heat flow in our case is the uncertainty of our observables
We are solving the steady state heat equation which requires input variables
where all of our prior estimates and uncertainties are added to the cost function.
The temperature,
we can formulate this into a forward model where the cost function is returned:
def f(x):
''' forward model, returns the cost function J(x) '''
k, H, Q0 = tuple(x)
T = T0 + Q0*z/k + H*dz**2/k * (1-exp(-z/dz))
q = k * (T[1]-T[0]) / (z[1]-z[0])
return (q-q_obs)**2/sigma_q**2 + (k-k_prior)**2/sigma_k**2 + (H-H_prior)**2/sigma_H**2 + (Q0-Q0_prior)**2/sigma_Q0**2
q_obs
and sigma_q
are surface heat flow observations and uncertainty respectively.
The right hand side of this equation is all our prior estimates of k
, H
, and Q0
are the only parameters that change.
Now the tangent linear model uses the derivatives of each line in the forward model with respect to k
, H
, and Q0
.
def f_tl(x, dx):
''' tangent linear model '''
k, H, Q0 = tuple(x)
dkdx, dHdx, dQ0dx = 1, 1, 1
dk, dH, dQ0 = tuple(dx)
# T = f(k, H, Q0)
T = T0 + Q0*z/k + H*dz**2/k * (1-exp(-z/dz))
dTdk = -Q0*np.array(z)/k**2 - H*dz**2/k**2 * (1-exp(z/dz))
dTdH = dz**2/k * (1-exp(z/dz))
dTdQ0 = z/k
dT = dTdk*dk + dTdH*dH + dTdQ0*dQ0 # similar to the dot product
# q_model = f(T0, T1, k)
q_model = k * (T[1]-T[0])/(z[1]-z[0])
dqdT1 = k / (z[1]-z[0])
dqdT0 = -k / (z[1]-z[0])
dqdk = (T[1]-T[0]) / (z[1]-z[0])
dq_model = dot([dqdT0, dqdT1, dqdk], [dT[0], dT[1], dk])
# return f(q_model, k, H, Q0)
d2q = (2*q_model-2*q_obs)/sigma_q**2
d2k = (2*k-2*k_prior)/sigma_k**2
d2H = (2*H-2*H_prior)/sigma_H**2
d2Q0 = (2*Q0-2*Q0_prior)/sigma_Q0**2
d2x = dot([d2q, d2k, d2H, d2Q0], [dq_model, dk, dH, dQ0])
return d2x
This returns reasonable values of
We backwards-propagate a small change in
def f_ad(x, df):
''' adjoint model '''
k, H, Q0 = tuple(x)
T = T0 + Q0*z/k + H*dz**2/k * (1-exp(z/dz))
q_model = k * (T[1]-T[0])/(z[1]-z[0])
d2q = (2*q_model-2*q_obs)/sigma_q**2
d2k = (2*k-2*k_prior)/sigma_k**2
d2H = (2*H-2*H_prior)/sigma_H**2
d2Q0 = (2*Q0-2*Q0_prior)/sigma_Q0**2
dq_ad = d2q*df
dk_ad = d2k*df
dH_ad = d2H*df
dQ0_ad = d2Q0*df
dqdT0 = -k / (z[1]-z[0])
dqdT1 = k / (z[1]-z[0])
dqdk = (T[1]-T[0]) / (z[1]-z[0])
dT0_ad = dqdT0*dq_ad
dT1_ad = dqdT1*dq_ad
dk_ad += dqdk*dq_ad
dTdk = -Q0*np.array(z)/k**2 - H*dz**2/k**2 * (1-exp(z/dz))
dTdH = dz**2/k * (1-exp(z/dz))
dTdQ0 = z/k
dk_ad += dTdk[0]*dT0 + dTdk[1]*dT1
dH_ad += dTdH[0]*dT0 + dTdH[1]*dT1
dQ0_ad += dTdQ0[0]*dT0 + dTdQ0[1]*dT1
return dk_ad, dH_ad, dQ0_ad
The adjoint model returns
Apart from being able to calculate T[0]
, T[1]
, T[2]
… T[n]
, but the adjoint only requires computation on elements you actually need - in this case only T[0]
and T[1]
.
When writing the adjoint, only enough of the forward model needs to be computed so that you have the appropriate variables. For most cases this is a complete forward run.