Implement indicator function in the weak formulation

340
views
0
10 months ago by
Hello everyone,

I have a term in my diffusion pde which represents the area where my function h is non-zero. I.E. ind(u) = 1 if u>=0 and ind(u) = 0 otherwise.

I have tried a couple of things and cannot figure out how to fit this into FEniCS. See the following:

def ind(u_n):
if u_n >=0:
return 1
else:
return 0

doesn't work, neither does

Can anyone help me? Full code attached at the end.

Thanks, Zhen.

from fenics import *
import time
import scipy

def q(u):
"return nonlinear coefficient"
return u**3

T = 2.0 # final time
num_steps = 20 # number of time steps
dt = T / num_steps # time step size
E = 0.1 # evaporation constant

# Create mesh and define function space
nx = ny = 30
mesh = RectangleMesh(Point(-3, -3), Point(3, 3), nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

# Define initial value
u_0 = Expression('x[0]>=-1 && x[0]<=1 && x[1] >=-1 && x[1] <= 1 ? exp(-3*a*pow(x[0], 2) - a*pow(x[1], 2)):0',
degree=2, a=5)

#u_0 = Constant(0)
u_n = interpolate(u_0, V)

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Constant(0)

#note whether there is evaporation or not
def ind(u_n):
if u_n >=0:
return 1
else:
return 0

# Create VTK file for saving solution

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
t += dt

# Compute solution
solve(F == 0, u, bc)

# Save to file and plot solution
vtkfile << (u, t)
plot(u)

# Update previous solution
u_n.assign(u)

# Hold plot
interactive()

4
10 months ago by
For your first approach consider this:

import numpy as np

def ind(u):
u_ind = Function(V)
pos = np.where(u.vector()>0.0)[0]
u_ind.vector()[pos] = 1.0

return u_ind.vector()

u_ind = Function(V)
u_ind.vector()[:] = ind(u_n)

And for the second one:

import numpy as np

u_ind = Function(V)
pos = np.where(u_n.vector()>0.0)[0]
u_ind.vector()[pos] = u_n.vector()[pos]
Thanks Hernan! May I ask how should I incorporate this into my F? Ultimately, I want to solve a pde where there is a term proportional to the area of the region where u is non zero.

best,
zhen
written 10 months ago by zhen shao
I was thinking in something like this:

u_ind = Function(V)

# Create VTK file for saving solution

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
t += dt

# Compute solution
pos = np.where(u_n.vector()>0.0)[0]
u_ind.vector()[pos] = u_n.vector()[pos]
solve(F == 0, u, bc)

.
.
.

​
written 10 months ago by Hernán Mella
That works! Thanks!

Thank you for helping me Mella, I am new to Python, do you recommend any online tutorial that I could go through?
written 10 months ago by zhen shao
Take a look here

Regards!
written 10 months ago by Hernán Mella