Fast normalization of a vector field

1

Hi!

I need to normalize a vector field (coming from the solution of another problem) and have come up with several ways of doing it in the following small example. If the norm is very small, the vector is set (arbitrarily) to (10,0,0).

from dolfin import *
import numpy as np
import time

np.set_printoptions(formatter={'float': '{: 0.6f}'.format})

N=10
mesh = BoxMesh(Point((0,0,0)),Point((1,1,1)),N,N,N)

for order in range(1,3):
    print('Using Lagrange element order',order)
    Vv = VectorFunctionSpace(mesh, "Lagrange", order)
    du=interpolate(Expression(("cos(3*x[0]+x[1]+x[2])+0.1", "sin(3*x[1]-x[0]-x[2])-0.5", "x[0]+x[1]+x[2]+0.2"), degree=order), Vv)

    print('Original du:           ',du.vector()[0:6])
    print()

    # Normalize using project
    from ufl import le,conditional
    t0=time.time()
    normalized_du=project( conditional( le(dot(du,du),1e-8),Constant((1.0e1,0.0,0.0)),du/sqrt(dot(du,du)) ), Vv, solver_type="gmres",preconditioner_type="hypre_amg")
    dt=dt=time.time()-t0
    print('project normalized dt: ',dt)
    print('project normalized du: ',normalized_du.vector()[0:6])
    print()

    # Normalize using numpy
    t0=time.time()
    norm=np.sqrt(du.vector()[0::3]**2+du.vector()[1::3]**2+du.vector()[2::3]**2)
    m=np.zeros_like(du.vector()[:])
    m[0::3]=np.where(norm[:]<1e-8, 10 , du.vector()[0::3]/norm[:] )
    m[1::3]=np.where(norm[:]<1e-8,  0 , du.vector()[1::3]/norm[:] )
    m[2::3]=np.where(norm[:]<1e-8,  0 , du.vector()[2::3]/norm[:] )
    normalized_du=Function(Vv)
    normalized_du.vector()[:]=m[:]
    normalized_du.vector().apply('')
    dt=dt=time.time()-t0
    print('numpy normalized dt:   ',dt)
    print('numpy normalized du:   ',normalized_du.vector()[0:6])
    print()

    t0=time.time()
    n=Expression(("(mx*mx+my*my+mz*mz)<1e-8 ? 10 : mx/sqrt(mx*mx+my*my+mz*mz) ","(mx*mx+my*my+mz*mz)<1e-8 ? 0 : my/sqrt(mx*mx+my*my+mz*mz) ","(mx*mx+my*my+mz*mz)<1e-8 ? 0 : mz/sqrt(mx*mx+my*my+mz*mz) "),mx=du.sub(0),my=du.sub(1),mz=du.sub(2), degree=order)
    normalized_du = interpolate(n, Vv)
    dt=dt=time.time()-t0
    print('interpol normalized dt:',dt)
    print('interpol normalized du:',normalized_du.vector()[0:6])
    print()

    t0=time.time()
    n=Expression(("(mx*mx+my*my+mz*mz)<1e-8 ? 10 : mx/sqrt(mx*mx+my*my+mz*mz) ","(mx*mx+my*my+mz*mz)<1e-8 ? 0 : my/sqrt(mx*mx+my*my+mz*mz) ","(mx*mx+my*my+mz*mz)<1e-8 ? 0 : mz/sqrt(mx*mx+my*my+mz*mz) "),mx=du.sub(0),my=du.sub(1),mz=du.sub(2), degree=order)
    normalized_du=Function(Vv)
    LagrangeInterpolator.interpolate(normalized_du,n)
    dt=dt=time.time()-t0
    print('Lagrange normalized dt:',dt)
    print('Lagrange normalized du:',normalized_du.vector()[0:6])
    print()

    if order==1:
        # Normalize using assemble
        # Only works with 1st order elements
        t0=time.time()
        normalized_du=Function(Vv)
        w=TestFunction(Vv)
        assemble( inner(   conditional( le(dot(du,du),1e-8), Constant((1.0e1,0.0,0.0)) , du/sqrt(dot(du,du)) )   , w )*dP,normalized_du.vector())
        dt=dt=time.time()-t0
        print('assemble normalized dt:',dt)
        print('assemble normalized du:',normalized_du.vector()[0:6])
        print()

I get the following output:

Using Lagrange element order 1
Original du:            [ 0.640302 -1.341471  1.200000  0.721610 -1.283327  1.100000]

project normalized dt:  0.08706331253051758
project normalized du:  [ 0.336907 -0.705458  0.630546  0.395111 -0.703114  0.602121]

numpy normalized dt:    0.0008580684661865234
numpy normalized du:    [ 0.335171 -0.702203  0.628149  0.392641 -0.698281  0.598529]

interpol normalized dt: 0.6311995983123779
interpol normalized du: [ 0.335171 -0.702203  0.628149  0.392641 -0.698281  0.598529]

Lagrange normalized dt: 0.03386521339416504
Lagrange normalized du: [ 0.335171 -0.702203  0.628149  0.392641 -0.698281  0.598529]

assemble normalized dt: 0.004649162292480469
assemble normalized du: [ 0.335171 -0.702203  0.628149  0.392641 -0.698281  0.598529]

Using Lagrange element order 2
Original du:            [ 0.508487 -1.367423  1.250000  0.721610 -1.283327  1.100000]

project normalized dt:  0.8406329154968262
project normalized du:  [ 0.264656 -0.711693  0.650627  0.393022 -0.698052  0.597951]

numpy normalized dt:    0.0018808841705322266
numpy normalized du:    [ 0.264676 -0.711765  0.650644  0.392641 -0.698281  0.598529]

interpol normalized dt: 2.368515729904175
interpol normalized du: [ 0.264676 -0.711765  0.650644  0.392641 -0.698281  0.598529]

Lagrange normalized dt: 0.19696259498596191
Lagrange normalized du: [ 0.264676 -0.711765  0.650644  0.392641 -0.698281  0.598529]

As you can see, the low-level numpy normalization is clearly the fastest (also in 100 mio.++ DOF parallel computations), but I’m here assuming the numpy array to have the DOFs ordered in a specific way, so I’m not really very fond of this way of normalizing the vector field as I don’t know if my assumptions regarding the inner workings of fenics, dolfin and/or petsc are correct.
(I could not get petsc’s VecNormalize to work…)

Is there a high-level way of normalizing the vector field, which would approach the same speed as the numpy normalization?

Best Regards,
Søren

1 Like