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)

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

return vec

def trapezoidal(a, b, M):

h = (b-a)/float(M)

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

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 …

## All 3 Replies

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.

commented: agree +15

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!

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.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts learning and sharing knowledge.