import matplotlib.pyplot as plt
import numpy as np

import pherosensor

from pheromone_dispersion.geom import MeshRect2D
from pheromone_dispersion.velocity import Velocity
from pheromone_dispersion.diffusion_tensor import DiffusionTensor
from pheromone_dispersion.diffusion_operator import Diffusion
from pheromone_dispersion.advection_operator import Advection, AdvectionAdjoint
from pheromone_dispersion.reaction_operator import Reaction
from source_localization.population_dynamique import PopulationDynamicModel
from source_localization.population_dynamique import StationnaryPopulationDynamicModel

Lx = 50 # the length along the x-axis
Ly = 60 # the length along the y-axis
Delta_x = 0.4 # 0.4 the space step along the x-axis
Delta_y = 0.4 # 0.4 the space step along the y-axis
T_final = 5 # the final time of the simulation
msh = MeshRect2D(Lx, Ly, Delta_x, Delta_y, T_final)
msh.calc_dt_implicit_solver(.1)

np.random.seed(0)
v1 = np.random.normal(0, 1, size=msh.x.size*msh.y.size)
v2 = np.random.normal(0, 1, size=msh.x.size*msh.y.size)

Test of the adjoint operators

This notebook aims to check numerically that the discretized adjoint operators are indeed the adjoint of the discretized direct operators.

np.random.seed(0)
x1 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*(msh.t_array.size-1))
x2 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*msh.t_array.size)

pop_dyn_model = PopulationDynamicModel(msh)
print("test of the adjoint operator of the population dynamic model operator:")
print("<v_1,Av_2> = ",np.dot(x1,pop_dyn_model.matvec(x2)))
print("<A*v_1,v_2> = ",np.dot(pop_dyn_model.rmatvec(x1), x2))
print("|<v_1,Av_2>-<A*v_1,v_2>| = ", np.abs(np.dot(x1,pop_dyn_model.matvec(x2))-np.dot(pop_dyn_model.rmatvec(x1),x2)))

test of the adjoint operator of the population dynamic model operator:
<v_1,Av_2> =  12503.694354331254
<A*v_1,v_2> =  12503.694354331188
|<v_1,Av_2>-<A*v_1,v_2>| =  6.548361852765083e-11

stat_pop_dyn_model = StationnaryPopulationDynamicModel(msh)
print("test of the adjoint operator of the stationnary population dynamic model operator:")
print("<v_1,Av_2> = ",np.dot(x1,stat_pop_dyn_model.matvec(x2)))
print("<A*v_1,v_2> = ",np.dot(stat_pop_dyn_model.rmatvec(x1), x2))
print("|<v_1,Av_2>-<A*v_1,v_2>| = ", np.abs(np.dot(x1,stat_pop_dyn_model.matvec(x2))-np.dot(stat_pop_dyn_model.rmatvec(x1),x2)))

test of the adjoint operator of the stationnary population dynamic model operator:
<v_1,Av_2> =  12455.477432781594
<A*v_1,v_2> =  12455.47743278158
|<v_1,Av_2>-<A*v_1,v_2>| =  1.4551915228366852e-11

Test of the adjoint of the advection operator

Recall that the advection operator is: \(A:c(x,y)\mapsto \nabla\cdot(U(x,y) c(x,y))~\forall (x,y)\in\Omega\) such that \(U(x,y)c(x,y)\cdot \vec{n} = 0~\forall(x,y)\in\partial\Omega_i\) with \(\partial\Omega_i = \{ (x,y)\in\partial\Omega, U(x,y)\cdot\vec{n}<0\}\).

On the other hand, the adjoint operator of the advection operator is \(A^*:c(x,y)\mapsto -U(x,y)\cdot\nabla c(x,y)=-\nabla\cdot(U(x,y) c(x,y))+c(x,y)\nabla\cdot U(x,y)~\forall (x,y)\in\Omega\) such that \(U(x,y)c(x,y)\cdot \vec{n} = 0~\forall(x,y)\in\partial\Omega_o\) with \(\partial\Omega_o = \{ (x,y)\in\partial\Omega, U(x,y)\cdot\vec{n}>0\}\).

U_hi = np.random.normal(0, 1, size=(msh.y_horizontal_interface.size,msh.x.size,2))
U_vi = np.random.normal(0, 1, size=(msh.y.size,msh.x_vertical_interface.size,2))

U = Velocity(msh, U_vi, U_hi)

A_flux_term = Advection(U, msh)
A_divU_term = Reaction(U.div, msh)

print("test of the adjoint operator of the advection operator:")
print("<v_1,Av_2> = ",np.dot(v1,A_flux_term.matvec(v2)))
print("<A*v_1,v_2> = ",np.dot(A_flux_term.rmatvec(v1)+A_divU_term.matvec(v1),v2))
print("|<v_1,Av_2>-<A*v_1,v_2>| = ", np.abs(np.dot(v1,A_flux_term.matvec(v2))-np.dot(A_flux_term.rmatvec(v1)+A_divU_term.matvec(v1),v2)))

test of the adjoint operator of the advection operator:
<v_1,Av_2> =  -676.6346071435416
<A*v_1,v_2> =  -676.6346071435415
|<v_1,Av_2>-<A*v_1,v_2>| =  1.1368683772161603e-13

Test of the adjoint of the diffusion operator

Recall that the diffusion operator is: \(D:c\mapsto \nabla\cdot(\mathbf{K}(x,y)\nabla c(x,y))~\forall (x,y)\in\Omega\) such that \(K\nabla c\cdot \vec{n} = 0~\forall(x,y)\in\partial\Omega\).

This operator is: \(D^*:c\mapsto \nabla\cdot(\mathbf{K}^T(x,y)\nabla c(x,y))~\forall (x,y)\in\Omega\) such that \(K^T\nabla c\cdot \vec{n} = 0~\forall(x,y)\in\partial\Omega\).

Let us note that the diffusion operator is auto-adjoint (\(D^*=D\)) for a symmetric diffusion tensor (\(\mathbf{K}=\mathbf{K}^T\)).

U_hi = np.ones((msh.y_horizontal_interface.size,msh.x.size,2))
U_hi[:,:,1]=0
U_vi = np.ones((msh.y.size,msh.x_vertical_interface.size,2))
U_vi[:,:,1]=0
U = Velocity(msh, U_vi, U_hi)

K_u_t = 0.01
K_u = 0.01
K = DiffusionTensor(U, K_u, K_u_t)

# change for a

D = Diffusion(K, msh)

print("test of the adjoint operator of the diffusion operator")
print("<v_1,Dv_2> = ",np.dot(v1,D.matvec(v2)))
print("<D*v_1,v_2> = ",np.dot(D.matvec(v1),v2))
print("|<v_1,Dv_2>-<D*v_1,v_2>| = ", np.abs(np.dot(v1,D.matvec(v2))-np.dot(D.matvec(v1),v2)))

test of the adjoint operator of the diffusion operator
<v_1,Dv_2> =  13.741038678354187
<D*v_1,v_2> =  13.741038678354208
|<v_1,Dv_2>-<D*v_1,v_2>| =  2.1316282072803006e-14

Test of the adjoint of the reaction operator

Recall that the reaction operator is: \(D:c\mapsto\tau_{loss}c\).

This operator is auto-adjoint: \(R^*=R\)

tau_loss = np.random.normal(0, 1, size=(msh.y.size,msh.x.size))

R = Reaction(tau_loss, msh)

print("test of the adjoint operator of the diffusion operator")
print("<v_1,Dv_2> = ",np.dot(v1,R.matvec(v2)))
print("<D*v_1,v_2> = ",np.dot(R.matvec(v1),v2))
print("|<v_1,Dv_2>-<D*v_1,v_2>| = ", np.abs(np.dot(v1,R.matvec(v2))-np.dot(R.matvec(v1),v2)))

test of the adjoint operator of the diffusion operator
<v_1,Dv_2> =  151.0434490129216
<D*v_1,v_2> =  151.0434490129216
|<v_1,Dv_2>-<D*v_1,v_2>| =  0.0

Test of the adjoint operators of the population dynamic model operators

Recall that the population dynamic model considered in the toolbox until now is \(\partial_t s = \gamma s\), and when \(\gamma=0\), we have a stationary population dynamic model.

np.random.seed(0)
v1 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*(msh.t_array.size-1))
v2 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*msh.t_array.size)

pdm = PopulationDynamicModel(msh)

print("test of the adjoint operator of the population dynamique")
print("<v_1,Dv_2> = ",np.dot(v1,pdm.matvec(v2)))
print("<D*v_1,v_2> = ",np.dot(pdm.rmatvec(v1),v2))
print("|<v_1,Dv_2>-<D*v_1,v_2>| = ", np.abs(np.dot(v1,pdm.matvec(v2))-np.dot(pdm.rmatvec(v1),v2)))

test of the adjoint operator of the population dynamique
<v_1,Dv_2> =  12503.694354331254
<D*v_1,v_2> =  12503.694354331188
|<v_1,Dv_2>-<D*v_1,v_2>| =  6.548361852765083e-11

spdm = StationnaryPopulationDynamicModel(msh)

print("test of the adjoint operator of the stationnary population dynamique")
print("<v_1,Dv_2> = ",np.dot(v1,spdm.matvec(v2)))
print("<D*v_1,v_2> = ",np.dot(spdm.rmatvec(v1),v2))
print("|<v_1,Dv_2>-<D*v_1,v_2>| = ", np.abs(np.dot(v1,spdm.matvec(v2))-np.dot(spdm.rmatvec(v1),v2)))

test of the adjoint operator of the stationnary population dynamique
<v_1,Dv_2> =  12455.477432781594
<D*v_1,v_2> =  12455.47743278158
|<v_1,Dv_2>-<D*v_1,v_2>| =  1.4551915228366852e-11