Skip to content
Snippets Groups Projects
Commit 24f34708 authored by MALOU THIBAULT's avatar MALOU THIBAULT
Browse files

Merge branch 'main' of forgemia.inra.fr:pherosensor/pherosensor-toolbox

parents 3e8b60d1 37ad67ef
No related branches found
No related tags found
No related merge requests found
Pipeline #224727 passed
source diff could not be displayed: it is too large. Options to address this: view the blob.
%% Cell type:code id:63d68a42 tags:
``` python
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
```
%% Cell type:code id:7ff10728 tags:
``` python
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)
```
%% Cell type:code id:a80eaeed tags:
``` python
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)
```
%% Cell type:markdown id:76fea742 tags:
# Test of the adjoint operators
This notebook aims at checking numerically that the implementation of the adjoint operators are indeed the adjoint of implementation of the direct operators.
This notebook aims to check numerically that the discretized adjoint operators are indeed the adjoint of the discretized direct operators.
%% Cell type:code id:864e983f-17b8-4780-a550-435e21cf6894 tags:
``` python
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)
```
%% Cell type:code id:00c8d8af tags:
``` python
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)))
```
%% Output
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
%% Cell type:code id:cd692889-4695-483e-bada-20c4b26ab801 tags:
``` python
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)))
```
%% Output
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
%% Cell type:markdown id:6b1fd7e3 tags:
## 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\}$.
%% Cell type:code id:ef95c7ff tags:
``` python
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)
```
%% Cell type:code id:cf57659f tags:
``` python
A_flux_term = Advection(U, msh)
A_divU_term = Reaction(U.div, msh)
```
%% Cell type:code id:321294b1 tags:
``` python
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)))
```
%% Output
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
%% Cell type:markdown id:1f3c694f tags:
## 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$).
%% Cell type:code id:c1e688de tags:
``` python
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
```
%% Cell type:code id:6aaabc1f tags:
``` python
D = Diffusion(K, msh)
```
%% Cell type:code id:3249df5b tags:
``` python
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)))
```
%% Output
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
%% Cell type:markdown id:1fafa200 tags:
## 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$
%% Cell type:code id:670c9549 tags:
``` python
tau_loss = np.random.normal(0, 1, size=(msh.y.size,msh.x.size))
```
%% Cell type:code id:0cecaaeb tags:
``` python
R = Reaction(tau_loss, msh)
```
%% Cell type:code id:b026834d tags:
``` python
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)))
```
%% Output
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
%% Cell type:markdown id:0696c32e tags:
## Test of the adjoint operators of the population dynamique model operators
## Test of the adjoint operators of the population dynamic model operators
%% Cell type:markdown id:41f35d06-b48a-4d36-935d-6da5a56e8ce9 tags:
Recall that the population
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.
%% Cell type:code id:3a58e1f7 tags:
``` python
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)
```
%% Cell type:code id:0f07c242 tags:
``` python
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)))
```
%% Output
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
%% Cell type:code id:2d4d495c tags:
``` python
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)))
```
%% Output
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
......
source diff could not be displayed: it is too large. Options to address this: view the blob.
source diff could not be displayed: it is too large. Options to address this: view the blob.
source diff could not be displayed: it is too large. Options to address this: view the blob.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment