Poisson adjoint solution is the same for both : solution & solution + noise

Legacy DOLFIN
1

Hello,

I am currently working on resolving a poisson equation for both forward and ajdoint problem. Here is my code :

from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
import numpy as np

mesh = UnitSquareMesh ( 50 , 50 )
V = FunctionSpace ( mesh , "Lagrange" , 1 )
# Define basis functions and parameters
u = TrialFunction ( V )
v = TestFunction ( V )
f = interpolate(Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2), V)
# Define variational problem
a = inner ( grad ( u ) , grad ( v ) ) * dx
L = f * v * dx
bc = DirichletBC (V , 0.0 , "on_boundary" )
# Solve variational problem
u = Function ( V )
solve ( a == L , u , bc )
plt.colorbar (plot(u))

ud = Function(V)

j = inner ( u - ud , u - ud ) * dx
J = assemble(j)

fc = Control(f)

dJdf = compute_gradient (J , fc)

# solution of the adjoint problem
tape = get_working_tape()
solve_block = tape.get_blocks()[1]
adj_u = solve_block.adj_sol.vector()[:]

adj_u_func = Function(V)
adj_u_func.assign(u)
adj_u_func.vector()[:]= adj_u
plot(adj_u_func)
plt.colorbar(plot(adj_u_func))

noise_level = 2
u_noise = Function(V)
u_noise.assign(u)
MAX = u_noise.vector().norm("linf")
noise_rand = noise_level * MAX * np.random.normal(0, 1, len(u_noise.vector().get_local()))

for i in range(len(u_noise.vector())):
    u_noise.vector()[i] = u_noise.vector()[i] + noise_rand[i]*u_noise.vector()[i]
   
plot(u_noise)
plt.colorbar(plot(u_noise))

j = inner ( u_noise - ud , u_noise - ud ) * dx
J = assemble(j)

fc = Control(f)

dJdf = compute_gradient (J , fc)

import numpy as np
# solution of the adjoint problem
tape = get_working_tape()
solve_block = tape.get_blocks()[1]
adj_u = solve_block.adj_sol.vector()[:]

adj_u_func_noise = Function(V)
adj_u_func_noise.assign(u_noise)
adj_u_func_noise.vector()[:]= adj_u
plot(adj_u_func_noise)
plt.colorbar(plot(adj_u_func_noise))

It’s weird that both the adjoint solution with and without noise are the same.

Am I misunderstaing something ?

PS : I think the tape call always get back to the adjoint solution without noise, even though we specify the noisy adjoint solution.

2

Adding some plots for better understanding.

Solution of the forward problem without noise :

image

3

Solution of the adjoint problem without noise :

image

4

Adding noise to the forward solution :

image

5

Adjoint solution for the noisy part (exactly the same as solution of the adjoint problem without noise):

image

6

You are assigning noise to the solution after solving the variational problem, thus the adjoint of the forward problem sees no noise and should be unaffected.

1 Like
7

Thank you, that was indeed the mistake !

8

Hello dokken,

I thought I’ve figured it out, but before going to the adjoint, I couldn’t see a change by adding noise before solving the variational problem.

from dolfin import *
from dolfin_adjoint import *
import matplotlib.pyplot as plt
import numpy as np

mesh = UnitSquareMesh ( 50 , 50 )
V = FunctionSpace ( mesh , "Lagrange" , 1 )
# Define basis functions and parameters
u = TrialFunction ( V )
v = TestFunction ( V )
f = interpolate(Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2), V)
# Define variational problem
a = inner ( grad ( u ) , grad ( v ) ) * dx
L = f * v * dx
bc = DirichletBC (V , 0.0 , "on_boundary" )
# Solve variational problem
u = Function ( V )
##### add noise to u ####
noise_level = 2
MAX = u.vector().norm("linf")
noise_rand = noise_level * MAX * np.random.normal(0, 1, len(u.vector().get_local()))

for i in range(len(u.vector())):
    u.vector()[i] = u.vector()[i] + noise_rand[i]*u.vector()[i]
#########################
solve ( a == L , u , bc )
plt.colorbar (plot(u))

Removing or leaving this part of the code, doesn’t make a difference in the u plot :

##### add noise to u ####
noise_level = 2
MAX = u.vector().norm("linf")
noise_rand = noise_level * MAX * np.random.normal(0, 1, len(u.vector().get_local()))

for i in range(len(u.vector())):
    u.vector()[i] = u.vector()[i] + noise_rand[i]*u.vector()[i]
#########################

Is there something I am not understanding ? (Maybe I should add noise to u when I call it as TrialFunction(V) ? I couldn’t do that since there is no way to access its values).

Thank you.