'Product' objects attributes, project function speed


298
views
0
11 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_bad = u_true *dp.exp(m_true)

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

1 Answer


0
11 months ago by
pf4d  
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.

Just trying to be helpful.








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 11 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 11 months ago by pf4d  
Also, to the two that downvoted, I hope to see two answers in response.
written 11 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »