0

My code tries to solve an optimization problem constrained to hyperbolic PDE by matching a numerical solution to a desired one. The basic idea is that: I've initial condition "U" as input which is being updated by this line: U = Initialize(Update_Design_Parameter(U)). I expect that updated U will affect all intermediate results throughout the code; when I print "Evolve_in_One_Timestep(U)" under the for loop (for i in range(5)), results show that Evolve_in_One_Timestep(U) is changing with updated U, but when I print "U_data", results show that it remains constant which also affects the gradient making it constant throughout! Can anyone help how to fix this problem?

```
""" Adjoint Programming """
from numpy import *
from math import *
from scitools.std import *
import glob,os
for filename in glob.glob('Optimizer*.png'):
os.remove(filename)
N = 1001; M = 75; Dx = 1.5/M ; Dt = 0.0005; a = 1.0
U = zeros(M)
def desired_U(z):
return (1-e**(z-0.5))/(1 + e**(z-0.5))
z = linspace(-20,20, M)
U_d = desired_U(z)
def U_initial_function (x):
value = cos(pi*(x/M))
return value
U = fromfunction(U_initial_function,(M,))
def Initialize(U):
return U
def Evolve_in_One_Timestep (U):
#create an array to store the new values of the solution
U_new = zeros((M,), float32)
#loop over every gridpoint, except for k=0 and k=M-1
for k in range (1, M-1):
derivative = -0.25*(U[k+1]**2 - U[k-1]**2)*(Dt/Dx) + 0.5*sqrt(a)*(U[k+1] - 2*U[k] + U[k-1])*(Dt/Dx)
U_new[k] = U[k] + derivative
#Apply terminal conditions
U_new[0] = U_new[1]
U_new[M-1] = U_new[M-2]
return U_new
#def Iterate_With_Time(U,N):
U_data = Initialize(U)
for n in range(0,N):
U_data = Evolve_in_One_Timestep(U_data)
#Terminal conditions
def Initial_Values():
L = U_data - desired_U(z)
return L
def Evolve(L):
#create an array to store the new values of the solution
U_T = U_data
new_U = zeros((M,), float32)
Backward = range(M-2, 0,-1)
Forward = range(1, M-1)
#loop over each gridpoint, except for k= M-1 and k = 0
for i, k in zip(Forward, Backward):
der = -(Dt/Dx)*(a)*U_T[i]*(L[k+1] - L[k-1])
new_U[k] = L[k] + der
#applying transversality boundary conditions
new_U[0] = new_U[1]
new_U[M-1] = new_U[M-2]
return new_U
Evolve_in_One_Timestep (U)
def Iterate_Adjoint(N):
U_Values = Initial_Values()
for z in range(N-1,-1, -1):
U_Values = Evolve(U_Values)
t = z*Dt
return U_Values
def f():
vec = U_data*Iterate_Adjoint(N)
return vec
def trapezoidal(a, b, M):
grad = f()
grad[0] /= 2.0
grad[-1] /= 2.0
h = (b-a)/float(M)
grad *= h
I = sum(grad)
return I
def Update_Design_Parameter(U):
DU0 = -0.25*trapezoidal(f()[0], f()[-1], M)
U0_new = DU0 + Initialize(U)
return U0_new
if __name__ == "__main__":
for i in range(5):
print "Gradient is:", trapezoidal(f()[0], f()[-1], M)
if abs((trapezoidal(f()[0], f()[-1], M))) < 10e-7:
print "======================================"
print "abs value is:" , abs ((trapezoidal(f()[0], f()[-1], M)))
print "======================================"
print "The value of U:", U
print "======================================"
break
else:
#Update U
U = Initialize(Update_Design_Parameter(U))
print Evolve_in_One_Timestep(U)
counter = 0; U_dat = Initialize(Update_Design_Parameter(U))
for n in range(0,N):
U_dat = Evolve_in_One_Timestep(U_dat)
if(n%100==0): #Plot the solution for every 100 timesteps
t = n*Dt
plotlabel= "t = " + str(n*Dt) #time labelling
plot(linspace(-0.37,0.37,75),U_d,'g0.3',linspace(-0.37,0.37,75),U_dat,'r0.3',label=plotlabel,legend='t=%4.2f' % t, hardcopy = 'Optimizer%04d.png' % counter )
counter += 1
movie('Optimizer*.png')
legend()
show()
raw_input('Press Enter:')
```