Generate synthetic data!

1

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
any help please i cant generate the synthetic data !!!

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
from dolfin_adjoint import  *

# Parameters
n = 20

mesh = UnitSquareMesh(n, n)  # "crossed"

# 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))

bcu_d1 = DirichletBC(W.sub(0), Constant((0, 0)), boundary_markers, 2)
bcu_d2 = DirichletBC(W.sub(0), Constant((0, 0)), boundary_markers, 3)
bc_problem = [bcu_d1, bcu_d2]

bc = bc_problem

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 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

u_true, p_true, goal_A, goal_B = Stokes_Solver(R_true)

xi = 1 - 2 * random()

ud = Function(W)
ud.assign(u_true)
# perturb state solution and create synthetic measurements ud
noise = Vector()
#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()))))
2

I would strongly suggest that you simplify your problem to a minimal code, following the examples in:

A general issue I see with your current approach is that you use compute_vertex_values()

As your function is not a CG 1 function, but a CG2 function, you should use ud.vector().get_local() to get the array of all degrees of freedom.

In general, I would suggest using a function assigner (Usage of FunctionAssigner - #2 by dokken) to assign the noise to the mixed solution.

As a last advice, note that you are not using a P2-P1 pair, but a P2-P2 pair:

3

Hi Dokken
this script work for me I don’t know is it true or not!

noise_level = 1
noise = Vector()
goal_A.init_vector(noise, 1)
noise_rand = noise_level * (1 - 2 * np.random.rand(len(ud.vector().get_local())))
print(noise_rand)
noise.set_local(noise_rand)
for i in range(len(ud.vector())):
    ud.vector()[i] = ud.vector()[i] + noise[i]*ud.vector()[i]