Normalized vector of trial function

variational formulation
1

Dear all,

I am dealing with a 3D nonlinear PDE which looks like, with u - scalar solution, f - vector solution:
div(f) = 0;
(a+b*norm2(f)**c)*f/norm2(f)=-grad(u)

where norm2(f) is the L2 norm. and a b c are predefined constant.

As you can see dividing f by its norm can cause dividing-by-zero, as a can be nonzero. In simple numpy scripts I can just do a if statement to make such case output a 0 vector. i.e. I tried

sele = FiniteElement(“CG”, mesh.ufl_cell(), 1)
vele = VectorElement(“CG”, mesh.ufl_cell(), 1)
mixed_space = FunctionSpace(mesh, MixedElement([sele, vele]))

U = Function(mixed_space)
u,f = split(U)
V = TestFunction(mixed_space)
v,tau = split(V)

def cfunc(f):
return conditional(le(norma(f),Constant(1e-6)),0f,(a+bnorma(f)**c)/norma(f))

However this woudn’t compile with error
raise self._exception_type(self._format_raw(*message))

UFLException: Shape mismatch between conditional branches.

How should I get around with this or does fenics have its all function to deal with it?

Thanks a lot in advance!

2

Hi, to wrap constants as vectors please try either:
as_vector((0,0)) or Constant((0,0)). However, Im not sure this will fix your problem. Could you please supply the full code?

3

The error is very likely in the line

return conditional(le(norma(f),Constant(1e-6)),0*f,(a+b* norma(f)**c)/norma(f))

and the error message you get is pretty verbose about it.
In the true path of your conditional you return a vector and in the false path you return a scalar.

Instead of 0*f just write Constant(0.)

Another simpler solution would be to add a very small number to the denominator for stabilization, e.g.

(a+b*norm2(f)**c)*f/(norm2(f)+1e-15)=-grad(u)

And please format your code properly. A multiplication symbol * will emphasize text in markdown and won’t be visible.

4

Hi,
Your equation looks like a mixed formulation of viscoplastic fluids. There are some specific methods which are well suited for this kind of problem that I can recommend if you are interested.
Best

5

Hi Klunkean, thank you very much for your answer, I wasn’t careful enough to catch the mismatch.
Also thank you for pointing out the posting format, it was my first time posting.

6

Hello bleyerj,
Thank you very much for your reply, it would be great if you can recommend it and I’d be sincerely grateful!

7

The problem comes from the a*f/norm2(f) term which is highly non-smooth and makes standard Newton solvers extremely difficult to use in practice.
One of the most popular way to deal with such an issue is through Augmented Lagrangian approaches although they may require a large number of iterations to converge. I recently proposed an interior-point method for this kind of problem which is much more efficient. These two approaches are discussed in the following paper: Advances in the simulation of viscoplastic fluid flows using interior-point methods
You can also find companion FEniCS code here: https://zenodo.org/record/1038520#.XQH3ICZS85k
Hope this will be useful

8

Hi bleyerj,

You are absolutely right! I have been wrestling with the this term for days and not having a way out. Thank you so much for the tremendous help!

9

Hi Bleyerj,

Thank you for your help! Although I haven’t implemented the method you suggested, I acquired some stable solution using Newton solver with relaxation parameter, at the cost of slow convergence.
My next step to impose a constraint over the maximum density, using lagrangian, and hopefully speed up the calculation.