Adjoint optimization

[Simplicity.]

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:')

[Gribouillis]

It seems that the only place where U_data is modified are lines 48 and 51, which are only executed once when the program is loaded, so U_data can't be affected by the code at the end of the file.

You could make your code more robust: 1) don't use global variables except for some constants N = 1001; M = 75; Dx = 1.5/M ; Dt = 0.0005; a = 1.0 ('a' should be given a more explicit name), 2) don't write code at module level, (lines 11, 16, 17, 48, 50, 51 could go in the if __name__ == "__main__" section), 3) instead, write code in functions (most part of the if __name__ = "__main__" section should go in one or more functions), 4) configure your editor to indent python code with 4 space characters instead of a tab.

[Simplicity.]

Thank you Gribouillis for your contribution. I understand your 2,3, and 4 points, for example I could have moved the part under if __name__ == "__main__": from counter = 0 to counter +=1 to the place under U_data = Evolve_in_One_Timestep(U_data), and simply change U_dat to U_data, this part of the code is creating a series of images at every iteration to make a movie, but I wanted it to do that at last when some criteria under if __name__ == "__main__": are met, or if not at the last cycle of for i in range(number). I don't get your first point clearly.

But the main issue remains, the point that I also noted earlier, U_data is just affected only once after the program is loaded! more advice is welcome on how to make it affected at every cycle!

[Gribouillis]

Thank you Gribouillis for your contribution. I understand your 2,3, and 4 points, for example I could have moved the part under if __name__ == "__main__": from counter = 0 to counter +=1 to the place under U_data = Evolve_in_One_Timestep(U_data), and simply change U_dat to U_data, this part of the code is creating a series of images at every iteration to make a movie, but I wanted it to do that at last when some criteria under if __name__ == "__main__": are met, or if not at the last cycle of for i in range(number). I don't get your first point clearly.

But the main issue remains, the point that I also noted earlier, U_data is just affected only once after the program is loaded! more advice is welcome on how to make it affected at every cycle!

I don't understand: if you don't write U_data = <something> at some place in the cycle, U_data can't change. Your program's logic is not very apparent.