I’m a new Fenics user and still getting used to formulating PDEs in weak form. While I understand all simple examples, there are terms that I don’t know what to do with.
For instance, let’s take a vector function that evolves like this
\dot u=\nabla_j \Sigma_{ij}
where \Sigma can be a stress in a NS equation. Now, if \Sigma has derivatives of other field, for instance a function like
where \phi is another field. What would be the best way of dealing with something like this in Fenics?
Since \Sigma has a second derivative, would I have to create another tensorial function space in which I solve for \Sigma? There I could solve the following weak form for a test function K_{ij}:
Is this the smartest way to solve this? Could it be done without having to use a new tensorial function space, and keeping a single weak form for u (and of course another for \phi)?
Since this is very related to nematic order parameters, are there any good Fenics examples out there for solving for nematic fields?
Since Σ has a second derivative, would I have to create another tensorial function space in which I solve for Σ ?
The standard approach to a higher-order PDE would indeed be a mixed method with an additional field. You might look at the Cahn–Hilliard demo for an example of this:
Could it be done without having to use a new tensorial function space, and keeping a single weak form for u (and of course another for ϕ )?
Just brainstorming now, you might consider a mixed method involving only a scalar-valued auxiliary field representing \nabla^2\phi (since that’s the only problematic term), instead of a tensor-valued field representing \Sigma.
You can also directly discretize higher-order PDEs without additional fields by using smoother basis functions. There has been some work on solving Cahn–Hilliard and other phase field problems using isogeometric analysis (IGA) (e.g., here), but that is not directly supported in FEniCS. I did write a library, tIGAr, to do IGA with FEniCS, but the usage differs somewhat from standard FEniCS.
Thank you for your answer, and thank you for the idea of using \nabla^2\phi as the auxiliary function, I’ll probably use that since it seems like the simplest solution.
I have a question though about your last remark of using smoother basis functions, does that mean I can always just use a higher degree function space to have a single equation? For instance, if in the Cahn-Hilliard example, I define my space as
V = FunctionSpace(mesh, "Lagrange", 2)
Then can I just solve the whole thing with a single equation having two derivatives in some of the fields and test functions? What would the difference be between using a mixed formulation with two spaces of degree 1 as in the tutorial, and the single space I wrote above?
No, what you’d need is more continuous derivatives between elements, not higher polynomial degree within each element. The Lagrange basis (of any polynomial degree) is still only C^0, i.e., derivatives up to 0^{\text{th}}-order are continuous. (You might take a look at my answer here for more discussion of inter-element continuity and variational formulations, in the simpler context of the Poisson equation.) The spline spaces from IGA allow you to control the number of continuous derivatives between elements. You would need at least C^1 continuity to discretize your PDE using a standard Bubnov–Galerkin method.
P.S. Another possible direction for higher-order PDEs is to use a discontinuous Galerkin method, such as the one used in the biharmonic tutorial, where there is a single C^0 field, but extra terms are added to the variational form, integrated on interior mesh facets. However, this is probably the most difficult to formulate.