### 'Product' objects attributes, project function speed

225
views
0
8 months ago by
Hi,

I am having issues with FEniCS structures and speed. I have a variable  $\psi=e^mu$ψ=emu , that I need to use/recreate frequently. I have no issues if I use psi = dp.project(dp.exp(m)*u, Vu ) to create  $\psi$ψ, but it runs slowly and I get numerical errors if I do not use dp.project for every instance of psi (it is for an optimization problem that requires multiple creations of psi multiple times). If I try to use psi = dp.exp(m)*u I get the error,
--->65 psi_bad.vector()
AttributeError: 'Product' object has no attribute 'vector'​

I have not been able to find any explicit documentation on what the attributes of 'Product' objects are; if someone could direct me to it I would appreciate it. The code below highlights the issues I am having.

import dolfin as dp
from mshr import *
import math
import numpy as np
dp.set_log_active(False)
#######################################################################################
# 1. Create the mesh and define the finite element spaces
x0 = 0.0; x1 = 1.0; y0 = 0.0; y1 = 1.0
nx = 10; ny = 10
mesh = dp.RectangleMesh(x0 , y0 , x1 , y1 , nx , ny)
Vm = dp.FunctionSpace( mesh, 'Lagrange' , 2 )
Vu = dp.FunctionSpace( mesh, 'Lagrange' , 2 )
#######################################################################################
# 2. Define the true parameter
m_true = dp.interpolate( dp.Expression( 'log( 2 + 6 * ( pow( x[0] - 0.5 , 2 ) + pow( x[1] - 0.5 , 2 ) < pow( 0.2 , 2 ) ) )' ) , Vm )
#######################################################################################
# 6. Generate the synthetic observations by solving the forward problem.

#Really u_true is the solution of a pde, but the code is not needed to illustrate
#the issues so I am making the field constant
u_true = dp.interpolate( dp.Expression( '2' ) , Vu )

psi_true = dp.project(u_true *dp.exp(m_true), Vu)

psi_array = np.exp(m_true.vector().array()[:]) * u_true.vector().array()[:]
psi = dp.Function( Vu )
psi.vector()[:] = np.copy( psi_array )

noise_level = 0.0

def generate_data(u, V, noise_level):
u_obs = dp.Function( V )
u_obs.assign( u )
# Perturb u and create synthetic observations
MAX = u_obs.vector().norm( "linf" )
noise = dp.Vector(u.vector())
noise.set_local( noise_level * MAX * np.random.normal(0 , 1 , len( u_obs.vector().array() ) ) )
u_obs.vector().axpy( 1. , noise )

return u_obs , MAX

[psi_obs, MAX] = generate_data(psi_true, Vu, noise_level)

proj_eq = ( dp.project(dp.exp(m_true)*u_true,Vu) - psi_obs) * dp.dx #<--- Works, but slow
proj = dp.assemble( proj_eq )
print "The norm with projection is ", np.linalg.norm(proj) #e-18

NO_proj_eq = ( dp.exp(m_true)*u_true - psi_obs) * dp.dx
NO_proj = dp.assemble( NO_proj_eq )
print "The norm without projection is ", np.linalg.norm(NO_proj) # e-05

NO_vect_eq = ( dp.exp(m_true)*u_true - psi_bad) * dp.dx
NO_vect = dp.assemble( NO_vect_eq )
print "The norm without vector attribute ", np.linalg.norm(NO_vect) # 0.0

array_eq = ( dp.exp(m_true)*u_true - psi) * dp.dx
array = dp.assemble( array_eq )
print "The norm of array is ", np.linalg.norm(array) #0.0441989439107

Thanks

Community: FEniCS Project

0
8 months ago by
OK, where to begin.

First, you don't need to define two function spaces, Vu and Vm are identical.

You do not need to import math; use fenics math functions with floats etc and numpy with numpy stuff.

Next, just make u_true = Constant(2)

The problem with taking psi_bad.vector() is that it has no method "vector".

If you want to add noise to the observations, take u_noisey.vector().set_local( u_true.vector().array() + eps ) where eps is a numpy array of noise of the same size as the degrees of freedom of Vm.  This will eliminate the need to project entirely.

Edit:

From the downvote I recieved, I assume this answer is insufficient.

Just project exp(m)*u_true once.  With this in hand, you can just add or subtract your psi terms using the local array values I mentioned with .set_local().

The Product is used by dolfin to discretize integrals.

When you assemble a form without including a test space, a scalar is returned.  So I'm not certain what you are trying to do with the np:norm.

I have already said that psi_bad had no attribute 'vector'; I want to know what attributes/methods it does have.
Also, how can you use psi_bad.vector().set_local( u_true.vector().array() + eps ) to add noise when as both you and I have pointed out psi_bad has no attribute 'vector'.
written 8 months ago by Katrina
Sorry, I am answering your question on my phone and the code is sort of difficult to make out.  See above edit.
written 8 months ago by pf4d
Also, to the two that downvoted, I hope to see two answers in response.
written 8 months ago by pf4d