Hi every one
please help me to dad a random noise to the stokes solution u_true=(u1,u2)
u=u_true + xi u_true
with xi is a random number between -1 and 1:
first, I want to solve the problem to get the true solution then I want to perturb the data u_true by a random noise to use it as synthetic data to solve the inverse problem
from random import random
import matplotlib.pyplot as plt
import numpy as np
from dolfin import *
from dolfin import Function
from mpl_toolkits import mplot3d
from mshr.cpp import *
from sympy import true
import nb
from dolfin_adjoint import *
"""
Unit Circle Mesh = mesh0
Annular Mesh = mesh1
Rectangle Mesh = mesh 2
Unit square Mesh = mesh 3
Circle Mesh = mesh 4
Rectangle mesh with hole = mesh5
"""
# noise level
noise_level = 100
fichier = open("data.txt", "w")
# Parameters
n = 2
x1, x2 = 1, 4
y1, y2 = 1, 2
r1 = 1
r2 = 2
x11, y11 = 0, 0
x22, y22 = 0, 0
# Domains
domain0 = Rectangle(Point(x1, y1), Point(x2, y2))
domain1 = Circle(Point(2, 1.5), 0.2)
domain2 = Rectangle(Point(2.5, 1.25), Point(3, 1.75))
domain3 = domain0 - domain1
domain4 = Circle(Point(0, 0), 1)
domain5 = Circle(Point(x11, y11), r1)
domain6 = Circle(Point(x22, y22), r2)
domain7 = domain6 - domain5
# Mesh
mesh0 = generate_mesh(domain4, n)
mesh1 = generate_mesh(domain7, n)
mesh2 = RectangleMesh(Point(x1, y1), Point(x2, y2), abs(x1 - x2) * n, abs(y1 - y2) * n) # "crossed"
mesh3 = UnitSquareMesh(n, n) # "crossed"
mesh4 = generate_mesh(domain5, n)
mesh5 = generate_mesh(domain3, n)
# mesh = Mesh_Choice()
mesh = mesh2
# plot(mesh)
# Define function spaces
P2 = VectorElement('P', triangle, 2)
P1 = FiniteElement('P', triangle, 2)
TH = MixedElement([P2, P1])
Q = FunctionSpace(mesh, P2)
Qp = FunctionSpace(mesh, P1)
W = FunctionSpace(mesh, TH)
n1 = FacetNormal(mesh)
# Define the RHS
force_null = Constant((0, 0))
R_true = 10
# Define boundary condition
def boundary(x, on_boundary):
return on_boundary
# Define boundary subdomains
tol = 1e-14
class Gamma_N(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], x1, tol)
class Gamma_R(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], x2, tol)
class Gamma_L1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], y1, tol)
class Gamma_L2(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], y2, tol)
class Gamma_C(SubDomain):
def inside(self, x, on_boundary):
ry = sqrt(pow(x[0], 2) + pow(x[1], 2))
return on_boundary and near(ry, 0.5)
# Mark boundaries
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
G_N = Gamma_N()
G_R = Gamma_R()
G_L1 = Gamma_L1()
G_L2 = Gamma_L2()
G_C = Gamma_C()
G_N.mark(boundary_markers, 0)
G_R.mark(boundary_markers, 1)
G_L1.mark(boundary_markers, 2)
G_L2.mark(boundary_markers, 3)
G_C.mark(boundary_markers, 4)
# Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
# Define boundary conditions
u_ex = Expression(("4*pow(x[1],3) - pow( x[0] ,2)", "4 * pow( x[0] ,3 )+ 2 * x[0]*x[1]-1"),
element=VectorElement("Lagrange", triangle, 5))
g_n = Expression(("24*x[1]* x[0]", "2*x[0]"),
element=VectorElement("Lagrange", triangle, 5))
p_ex = Expression("24*x[0]*x[1]- 2 * x[0]",
element=FiniteElement("Lagrange", triangle, 5))
"""
p_eex = interpolate(p_ex,Qp)
u_eex = interpolate(u_ex,Q)
du_ex_dn = assemble(dot(grad(u_eex),FacetNormal(mesh))*ds(0))
"""
bc_total = DirichletBC(W.sub(0), Constant((0, 0)), boundary)
bcu_n_r = DirichletBC(W.sub(0), u_ex, boundary_markers, 1 and 0)
# for verifying the errors results we choose bc_dex and bc_pex
bcu_dex = DirichletBC(W.sub(0), u_ex, boundary)
bc_pex = DirichletBC(W.sub(1), p_ex, boundary)
bc_test = [bcu_dex, bc_pex]
# for the main code we choose our BC bcu_d and bc_p or bc_pex
bcu_d1 = DirichletBC(W.sub(0), Constant((0, 0)), boundary_markers, 2)
bcu_d2 = DirichletBC(W.sub(0), Constant((0, 0)), boundary_markers, 3)
bcu_d3 = DirichletBC(W.sub(0), Constant((0, 0)), boundary_markers, 4)
bc_p1 = DirichletBC(W.sub(1), Constant(0), boundary_markers, 2)
bc_p2 = DirichletBC(W.sub(1), Constant(0), boundary_markers, 3)
bc_problem = [bcu_d1, bcu_d2]
bc_ch = [bcu_d1, bcu_d2, bcu_d3]
bc = bc_problem
# bc = BC_Choice()
def Stokes_Solver_TEST():
# Define trial and test functions
(us, ps) = TrialFunctions(W)
(v, q) = TestFunctions(W)
a1 = inner(grad(us), grad(v)) * dx - div(v) * ps * dx + q * div(us) * dx
L1 = inner(force_null, v) * dx
# Sample parameters for pressure stabilization
diam = CellDiameter(mesh)
beta = 0.2
delta = beta * diam ** 2
# The additional pressure stabilization terms
a1 += delta * inner(grad(ps), grad(q)) * dx
L1 += delta * inner(force_null, grad(q)) * dx
# Solve system
Us = Function(W)
# solve(a1 == L1, U, bc)
solve(a1 == L1, Us, bc)
# Get sub-functions
us, ps = Us.split(True)
us_error_L2 = errornorm(u_ex, us, 'L2')
print("L2_Err for U = || u_h - u_ex ||_L2 = ", us_error_L2)
us_error_H1 = errornorm(u_ex, us, 'H1')
print("H1_Err for U = || u_h - u_ex ||_H1 = ", us_error_H1)
ps_error_L2 = errornorm(p_ex, ps, 'L2')
print("L2_Err for p = || p_h - p_ex ||_L2 = ", ps_error_L2)
# Compute L2 error of the divergence
Ms_div = div(us) * div(us) * dx
div_err_us = sqrt(assemble(Ms_div))
print("L2_Err for div(u_h) = || div(u_h) ||_L2 = ", div_err_us)
plt.subplot(3, 1, 1)
plot(us)
# ncols stays as 1
plt.subplot(3, 1, 2)
plot(ps)
plt.subplot(3, 1, 3)
plot(mesh)
plt.show()
return us, ps
def Stokes_Solver(ROBIN):
# Define trial and test functions
(us, ps) = TrialFunctions(W)
(v, q) = TestFunctions(W)
a1 = inner(grad(us), grad(v)) * dx - div(v) * ps * dx + q * div(us) * dx + ROBIN * inner(us, v) * ds(1)
L1 = inner(force_null, v) * dx + inner(g_n, v) * ds(0)
# Sample parameters for pressure stabilization
diam = CellDiameter(mesh)
beta = 0.002
delta = beta* diam ** 2
# The additional pressure stabilization terms
a1 += delta * inner(grad(ps), grad(q)) * dx
L1 += delta * inner(force_null, grad(q)) * dx
# solve the forward/state problem to generate synthetic observations
goal_Au, goal_bu = assemble_system(a1, L1, bc)
# Solve system
U = Function(W)
# solve(a1 == L1, U, bc)
solve(a1 == L1, U, bc)
# print('Solving Done')
# Get sub-functions
us, ps = U.split(True)
# Compute L2 error of the divergence
M_div = div(us)*div(us)*dx
div_err = sqrt(assemble(M_div))
print("L2_Err for div(u_h) = || div(u_h) ||_L2 = ", div_err)
return us, ps, goal_Au, goal_bu
def Stokes_prime_Solver(ROBIN):
d = 1
# Define trial and test functions
(ups, p1s) = TrialFunctions(W)
(vp, q1) = TestFunctions(W)
a2 = inner(grad(ups), grad(vp)) * dx - div(vp) * p1s * dx + q1 * div(ups) * dx + ROBIN * inner(ups, vp) * ds(1)
L2 = inner(force_null, vp) * dx - d * inner(u_true, vp) * ds(1)
# Sample parameters for pressure stabilization
diam = CellDiameter(mesh)
beta = 0.2
delta = beta * diam ** 2
# The additional pressure stabilization terms
a2 += delta * inner(grad(p1s), grad(q1)) * dx
L2 += delta * inner(force_null, grad(q1)) * dx
# Solve system
Up = Function(W)
# solve(a1 == L1, U, bc)
solve(a2 == L2, Up, bc)
# print('Solving Done')
# Get sub-functions
ups, p1s = Up.split(True)
# Compute L2 error of the divergence
M_div_prime = abs(div(ups)) * dx
div_err_prime = assemble(M_div_prime)
print("Div_err_prime=", div_err_prime)
return ups, p1s
def Stokes_adj_Solver(ROBIN):
h = 1
# Define trial and test functions
(uas, pds) = TrialFunctions(W)
(va, qa) = TestFunctions(W)
a3 = inner(grad(uas), grad(va)) * dx - div(va) * pds * dx + qa * div(uas) * dx + ROBIN * inner(uas, va) * ds(1)
L3 = inner(force_null, va) * dx - h * inner(u_true, va) * ds(0)
# Sample parameters for pressure stabilization
diam = CellDiameter(mesh)
beta = 0.2
delta = beta * diam ** 2
# The additional pressure stabilization terms
a3 += delta * inner(grad(pds), grad(qa)) * dx
L3 += delta * inner(force_null, grad(qa)) * dx
# Solve system
Ua = Function(W)
# solve(a1 == L1, U, bc)
solve(a3 == L3, Ua, bc)
# Get sub-functions
uas, pds = Ua.split(True)
# Compute L2 error of the divergence
M_div_adj = div(uas) * div(uas) * dx
div_err_adj = sqrt(assemble(M_div_adj))
print("Div_err_Adj=", div_err_adj)
return uas, pds
if bc == bc_test:
u, p = Stokes_Solver_TEST()
exit()
else:
u_true, p_true, goal_A, goal_B = Stokes_Solver(R_true)
up, p1 = Stokes_prime_Solver(R_true)
ua, pa = Stokes_adj_Solver(R_true)
xi = 1 - 2 * random()
ud = Function(W)
ud.assign(u_true)
# perturb state solution and create synthetic measurements ud
MAX = ud.vector().norm("linf")
noise = Vector()
goal_A.init_vector(noise, 1)
#noise.set_local(noise_level * MAX * np.random.normal(0, 1, len(ud.compute_vertex_values())))
noise.set_local(noise_level * (1 - 2 * np.random.rand(len(ud.compute_vertex_values()))))
# bc_adj.apply(noise)
ud.vector().axpy(1., noise)
plt.subplot(2, 1, 1)
plot(u_true, "State solution with atrue")
# ncols stays as 1
plt.subplot(2, 1, 2)
plot(ud, "Synthetic observations")
plt.show()
plt.subplot(4, 1, 1)
plot(u_true)
# ncols stays as 1
plt.subplot(4, 1, 2)
plot(up)
plt.subplot(4, 1, 3)
plot(ua)
plt.subplot(4, 1, 4)
plt.show()
I got this error
/usr/bin/python3.8 /home/zayeni211/Desktop/PycharmProjects/Article-1/Robin_stokes/python_Testes/matrix.py
Invalid MIT-MAGIC-COOKIE-1 keySolving linear variational problem.
L2_Err for div(u_h) = || div(u_h) ||_L2 = 0.31250464177883164
Solving linear variational problem.
Div_err_prime= 0.03061035317689818
Solving linear variational problem.
Div_err_Adj= 0.05712158395892892
Traceback (most recent call last):
File "/home/zayeni211/Desktop/PycharmProjects/Article-1/Robin_stokes/python_Testes/matrix.py", line 370, in <module>
ud.vector().axpy(1., noise)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to perform axpy operation with PETSc vector.
*** Reason: Vectors are not of the same size.
*** Where: This error was encountered inside PETScVector.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
Process finished with exit code 1
thank you