Hi everyone,
I am currently running a coupled diffusion-advection problem. The code solves fine but there are negative values (numerical errors) that I would like to minimise if possible. Is there a way to set a condition that does not allow for negative solutions? I’ve attached an image to illustrate the problem. Any advice would be greatly appreciated.
Maybe take a look at this demo for enforcing bound constraints:
I have had similar problem when my mesh is under-resolved and the diffusivity is low. If the problem disappears when you increase diffusivity you might need higher mesh resolution in the problem area.
Just to check, you are using stabilization, e.g., SUPG, right? (Residual-based stabilization converges optimally even if the element Peclet number is low, so there’s not really any reason not to use it.)
Thanks for getting back to me. I’m quite new to the finite element method so I’m not too sure how to implement stabilization for my problem. Currently, the problem I’m trying to solve is as follows:
Taking C_f to be u1 and B to be u2 I implement the following :
#Define function space for system of concentrations
P1 = FiniteElement('P', mesh.ufl_cell(), 1)
element = MixedElement([P1,P1])
V = FunctionSpace(mesh, element)
#Define test functions
v1,v2 = TestFunction(V)
#Define functions for concentration
u = Function(V)
u_n = Function(V)
u1,u2 = split(u)
u_n1,u_n2 = split(u_n)
#Define terms
Da = 50000
D_T = 2.041*pow(10,-4)
T = 0.5
B_M = 0.356*pow(10,-6)
K_d = 0.0326*pow(10,-6)
ka = (D_T*Da)/(B_M*(T*T))
kd = ka*K_d
k1 = 0.009221
A1 = 0.24
Z_w = 966.21*pow(10,3)
tol = 1E-9
sim_time = 60.0
time_steps = 600
dt = sim_time/time_steps
#Define Initial Conditions
u_n1 = interpolate(Constant(0),V.sub(0).collapse())
u_n2 = interpolate(Constant(0),V.sub(1).collapse())
#Define Boundary Conditions
Flux = Expression("t<30+1e-9 ? ((k1*A1)/Z_w)*exp(-1*k1*t) : 0", degree=2, k1=k1, A1=A1, Z_w=Z_w, t=0)
u1_D = Constant(0)
boundaryconditions = {
1 : {"Neumann" : Flux},
2 : {"Dirichlet" : u1_D}
}
ds = Measure('ds', domain=mesh, subdomain_data=mf)
bcs = []
for i in boundaryconditions:
if "Dirichlet" in boundaryconditions[i]:
# bc = DirichletBC(V, boundaryconditions[i]["Dirichlet"], boundaries, i)
bc = DirichletBC(V.sub(0), boundaryconditions[i]["Dirichlet"], boundaries, i) # Dirichlet is applied only to 1 subspace according to the paper
bcs.append(bc)
integrals_N = []
for i in boundaryconditions:
if "Neumann" in boundaryconditions[i]:
if boundaryconditions[i]["Neumann"] != 0: #if not zero flux
g = boundaryconditions[i]["Neumann"]
integrals_N.append(v1*g*ds(i))
#Define Variational Problem
F = (((u1-u_n1 + u2-u_n2)*v1)/dt)*dx + D_T*dot(grad(u1),grad(v1))*dx + (((u2-u_n2)*v2)/dt)*dx - sum(integrals_N) \
- v2*ka*u1*B_M*dx + ka*u1*v2*u2*dx + kd*v2*u2*dx
#Solve problem
t = -1*dt
J = derivative(F, u)
for n in range(time_steps+1):
t += dt
Flux.t = t
# solve(F == 0, u, bcs, J=J)
solve(F == 0, u, bcs, J=J,
solver_parameters={"newton_solver":{"absolute_tolerance":1e-30,
"linear_solver": "gmres"}})
u1,u2 = u.split()
vtkfile_u1 << (u1,t)
vtkfile_u2 << (u2,t)
u_n.assign(u)
The example I’m familiar with with regards to SUPG as seen here deals with velocity. I’m not sure how to implement stabilization to my problem. I suppose its not accurate to refer to this as a advection problem due to the lack of a velocity field; however, I do intend on introducing one down the line.
Are there any examples I can study to stabilize my current problem? I’m not too familiar with stabilization methods and how to employ them. Any advice would be greatly appreciated and thanks for taking the time.
Also, I’m not too sure if I’ve implemented the flux (time dependent) correctly. The magnitude of the concentration in the mesh drops significantly after 30s when the flux changes to zero. Which should not be the case.
Hi thanks for taking the time. Yes the problem does resolve with a much larger diffusivity but I believe, for my problem, the refinement required to resolve this might not be possible.
Hey,
I would suggest that you take a look at this: https://home.iitm.ac.in/kanupindi/course_projects/me697f_report.pdf
Here, Galerkin Least Squares (GLS) stabilization is introduced in an abstract setting, and examples are given for flow problems. You can interpret the popular stabilization methods like SUPG, PSPG or LSIC (for Navier-Stokes equations) as certain terms of the GLS stabilization. So that should be a good starting point for you.
Hi, thanks for the reference I’ll certainly look into it. However, from what I’ve gleaned these stabilization methods are only implemented for convection problems where there is a velocity field. Here the problem is coupled and non-linear diffusion (changed the titled) so these stabilization methods may not work. I was wondering if there are methods, implemented in FEniCS, for stabilizing a pure diffusion problem? Also, sorry for mislabeling the problem as an advection problem.
You could consider doing a coordinate transform to solve for lnC, and find c=exp(lnC). It b comes impossible to get negative c values but it does change the properties of the finite element approximation somewhat.
Hi Mike,
Thanks, I will give this a try. But it turns out the error for my case was with the way I implemented the initial conditions. I specified it as a function and went on to interpolate a constant value. This prevented it from being updated with the assign function. So either I specify it as a vector and then interpolate or I specify it as a function and not assign a value since it initializes it to be zero by default. The numerical oscillation is very small and does disappear with a more refined mesh as suggested by @fkag2. I think the stabilization functions are more useful with diffusion-advection problems.
Best Regards,
Kart
If you want to assign a function with a constant value you can always do
c = 1.0
u = Function(V)
u.vector()[:] = c
atleast for Lagrangian finite elements.
Thanks. Since I am using subspaces:
u_n = Function(V)
u_n1, u_n2 = split(u_n)
For now I can just get away with this since I need to initialize to zero anyways, which is what it defaults to.
For other values (non-zero) I can do so with an expression but as a vector
u_0 = Expression (('E1','E2'), degree=1) # E1 and E2 are expressions
u_n = interpolate(u_0,V)
u_n1, u_n2 = split(u_n)
But it might be worth testing if I can do:
c = 1.0
u_n = Function(V)
u_n1, u_n2 = split(u_n)
u_n1.vector()[:] = c
u_n2.vector()[:] = c
Thanks for the tip 
If you want to have the same constant for both parts, just use
u_n.vector()[:] = c
before you use “split” in your last example.
For different values have a look here: Interpolation of Indexed Functions of Mixed Elements
oh damn, I thought I had invented that approach myself…! For every clever thought, there is a published clever thought that precedes it lol.
I got it from the Console multiphysics manual… 
Actually I found it doesn’t always behave well. My colleague explained that the finite element representation on lnC might not be as realistic as just C…
It seems to be working ok for my application, though I could well believe it might introduce trouble elsewhere under the right conditions. The reason I introduced was not to maintain positive c (that seemed to be happening anyway), but to maintain consistent convergence criteria. I was finding that with unscaled c, I had to readjust convergence tolerances each time when I changed boundary conditions (in particular BCs with large magnitudes in nonlinear system) or reference c0.
I used ln c to help construct a scaled version of the solution that had values mainly between 0 and 1, for which the same convergence tolerance works in all cases.