How to calculate linf norm of an expression both globally and locally and update it during calculation?

Question

Hi All,

I try to solve an stabilized advection dominated transient transport equation (my previous post about SUPG. In addition to SUPG, I am currently experimenting artificial viscosity as indicated in this tutorial (dealii step 31).

$R_{\alpha} = (\frac{\partial T}{\partial t} + \bf{u} \cdot \nabla T - \nabla \cdot k \nabla T) T^{\alpha - 1}$ with additional viscosity(diffusion) term as $\nu_{\alpha} = \beta ||u||_{l\inf(K)} \min({h_K, h^{\alpha}_K \frac{||R_{\alpha}||_{linf(K)}}{c_R||u||_{l\inf(\Omega)} var(T) diam(\Omega)^{\alpha - 2}}})$.

As you can see, the additional viscosity term dependents on local residual on cell $K$, and is large where residual is large and small where residual is small. This involves calculation of $||u||_{l\inf}$ both in cell $K$, and in entire domain $\Omega$, and $||R_{\alpha}||_{linf(K)}$ for cell $K$.

Here are my problems, and I am sorry if it is too basic.

1. How do I calculate velocity norm using expression ? Here is my attempt, but it does not work. I have also tried many other methods (project instead of interpolate, add .array() in the end..... etc), it does not work. I wonder what is the appropriate way to calculate norm given a expression ?

   mesh = RectangleMesh(Point(-1,-1), Point(1, 1), 100, 100, "crossed")
   W = VectorFunctionSpace(mesh, "CG", 1)
   velocity = Expression(("x[1]" , "-x[0]"), domain = mesh, degree=3) # convecting velocity
   velocity_vector = interpolate(velocity, W)
   velocity_norm = norm(velocity_vector,'linf')
   

2. My guess above calculation is for entire domain. I wonder is there any way to calculate locally? Further, because residual changes respect to time and iterations. Thus, do I need to implement it during the time stepping loop ?

3. I have the similar question regarding to the output of cellSize(mesh). Is it a local mesh size, or averaged global mesh size ?

Any help will be greatly appreciated! I have attached my complete code as follow:

from dolfin import *

#read mesh info
mesh = RectangleMesh(Point(-1,-1), Point(1, 1), 100, 100, "crossed")

V = FunctionSpace(mesh, "CG", 1)
W = VectorFunctionSpace(mesh, "CG", 1)

tol = 1e-14

#  initial condition
ic= Expression("((pow(x[0]+0.3,2)+pow(x[1]+0.3,2))<0.2*0.2)?(1.0):(0.0)", domain=mesh, degree = 3)

mesh = RectangleMesh(Point(-1,-1), Point(1, 1), 100, 100, "crossed")
velocity = Expression(("x[1]" , "-x[0]"), domain = mesh, degree=3) # convecting velocity
velocity_vector = interpolate(velocity, W)
velocity_norm = norm(velocity_vector,'linf')


# Define trial and test function and solution at previous time-step
u = Function(V) # set as nonlinear term
v = TestFunction(V)
u0 = Function(V)

#initial condition
ic= Expression("((pow(x[0]-0.2,2)+pow(x[1]-0.2,2))<0.1*0.1)?(1.0):(0.0)", domain=mesh, degree = 3)
u0.interpolate(ic)

# Crank Nicolson
u_mid = 0.5*(u0 + u)

# Equation coefficients
c = Constant(0.000000001)
h = CellSize(mesh)
t_end = 0.5
dt = 0.001

# Calculation of artifical viscosity
alpha = 1.0 # stablization parameters
beta = 0.017 * 2.0 # stablization parameters
R = abs((u - u0 + dt*(dot(velocity, grad(u_mid)) - c*div(grad(u_mid))))*pow(u_mid, alpha - 1.0)) # residual
c_R = pow(2.0, (4.0 - 2.0*alpha)/2.0) # constant
var_T = u.vector().array().max() - u.vector().array().min()# Note: u max - u min over entire domain
h_alpha = pow(h, alpha)*norm(R, 'linf')/(c_R*norm(velocity_vector,'linf')* var_T * pow(h, alpha - 2)) # Note: Here Linf of velocity is over entire domian but norm (R) should be local to each cell
nu = beta * norm(velocity_vector,'linf') * min(h, h_alpha)

# Galerkin variational problem
F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx + (c + nu)*dot(grad(v), grad(u_mid))*dx)


# Time-stepping
t = 0.0
k = 0
while t < t_end:

    # Solve the problem
    solve(F == 0, u)

     # Move to next time step
    u0.assign(u)
    t += dt
    k += 1

Answer 1

1. The 'linf' norm applies only to vectors (and just computes the max of the abs), so you need velocity_norm = norm(velocity_vector.vector(),'linf')
2. If you want to calculate this locally, find the indices of the dofs you need, then index the array with them local_values = velocity_vector.vector().array()[dof_indices], then you can compute local_velocity_norm = np.max(np.abs(local_values)))
3. CellSize() returns a mesh cell Function, i.e. one whose value varies with the cell. I think it uses the circumradius, but check the documentation.

Comments

Hi mdbenito, thanks for taking the time answering the question. I have a few more general question regarding the "data type" used in fenics.

1. As you can see, I project my velocity expression onto the VectorFunctionSpace as "velocity_vector". But as indicated by your answer, it is not a "vector". And in order to get the vector, I need "velocity_vecotor.vector()". So my question is what is the data type of "velocity_vector" ?

2. To setup the formulation, I used following line code:

   F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx + (c + nu)*dot(grad(v), grad(u_mid))*dx)
   

   As you notice, I directly use the velocity expression instead of "vector". So I guess my question is what is the "data type" of the function expression "F" ? Or I should say what is the acceptable data type for "Solve" function. The dolfin documentation seemdoes not explicitly indicate this part.
   3

---

Sorry, due to some browser error. I split third part into following:

I see you point about get to the array and loop through the dof_index. Since this part changes with my local solution, time. Do you think I could create a mesh cell function to let it update automatically ? I think my question is that I see how this value can be calculated now. But How exactly can I update this local value back to "nu" at every time/iteration ?

Again, thanks !

---

1. velocity is a Function over the finite element space.
2. F is an UFL Form.

3. A CellFunction is attached to a Mesh, not any Function. However you might be able to define one which marks cells where your solution attains some specific value of interest, then compute any forms over those cells.

   class InterestingPart(SubDomain):
       def __init__(self, solution, threshold):
           self.u = solution
           self.threshold = threshold
           super().__init__()
   
       def inside(self, x, on_boundary):
           return u(x) > self.threshold
   
   # Assuming u is your (scalar) solution over FunctionSpace V...
   
   interesting_domains = CellFunction("uint", V.mesh())
   marker = InterestingPart(u, 0)
   marker.mark(interesting_domains, 1)
   
   integral = u*dx(subdomain_data=interesting_domains, subdomain_id=1)
   assemble(integral)
   

   I guess this will be quite slow, though.

I strongly suggest that you read the FEniCS tutorial in order to clarify all those questions about data types, etc.

Good luck!