Higher derivatives of trial and test functions.

[martenlienen]

Hey, I am trying to solve equations such as the Korteweg-de Vriews equation $u_t + uu_x + u_{xxx} = 0$ and the Kuramoto-Shivashinsky equation $u_t + u_{xx} + u_{xxxx} + uu_x = 0$ in 1D. With ElementLinePp I can create a trial space with enough smoothness for 3rd and 4th derivatives and I have seen that DiscreteField also has fields for 3rd through 6th derivatives. However, they are not used or filled anywhere except for ElementGlobal.

Is it possible to solve these equations with scikit-fem as of today? If not, what would it take?

[kinnala]

I think it would be possible to create Element classes for such purposes. What kind of variational formulation are you envisioning? What are the boundary conditions and the time discretization techniques like?

[kinnala]

It looks to me that ElementLineHermite could be enough to solve either of those. However, these are nonlinear equations and there could be some special considerations, e.g., irregular points that are not well suited to such basis functions.

[martenlienen]

I would like to solve in time with adaptive time stepping. On the boundary, I would apply Dirichlet conditions. For the variational formulation, I don't have any special considerations, just the simplest possible <v, rhs(u)> for a start. I will check out ElementLineHermite, thanks.

[kinnala]

The variational formulation decides what kind of boundary conditions are natural and essential. It is perhaps misleading to talk about 'just' Dirichlet boundary conditions as whenever the order of the boundary value problem increases, so does the number of boundary conditions required to render the solution unique. E.g., for Euler-Bernoulli beam equation $u_{xxxx} = 1$ one can set various orders of the derivatives to given values on the boundary and it is not clear what is "Dirichlet condition" as one example could be $u=u_{xx}=0$ on the left and $u_{x}=u_{xxx}=0$ on the right. Technically both are somewhat similar to Dirichlet conditions as we are prescribing DOFs.

Concerning time discretization, I was more into knowing whether you are planning on implicit or explicit time stepping? We don't have anything for time stepping here but you can discretize the space part with scikit-fem and deal with the time part however you like. Nevertheless, in case of explicit time stepping you might be able to cheat a bit in the linearization of the problem. But I guess it comes back to the properties of the PDE's and how strict they are about details such as energy preservation.

If you have an example setup, e.g., a reference problem, I could also try it out.

[martenlienen]

Thanks for the explanation! I will have to take a deeper look to understand, which of these pieces are relevant for me. I don't have a specific problem to solve and I am rather just poking around while thinking about a new research problem.

[kinnala]

Here is an example with periodic boundary conditions:

https://www.chebfun.org/examples/pde/KSWave.html

We have quite recently added support for periodic meshes, perhaps it might work for this.

[kinnala]

Unfortunately seems that MeshLine1DG was not tested with ElementGlobal, see #902. Otherwise I got this to run, however, I don't have a good test case for it.

[kinnala]

I think the periodicity is somehow important here. With fixed boundary conditions I'm only getting solutions that go rather quickly to zero. Here is the code I used, it requires master because that automatic differentiation thing is not released:

from skfem import *
from skfem.experimental.autodiff import NonlinearForm
from skfem.experimental.autodiff.helpers import *
import numpy as np


m = MeshLine(np.linspace(0, 1, 200))
basis = Basis(m, ElementLineHermite())

x = basis.project(lambda x: np.cos(4. * np.pi * x[0]) - 1)
xinit = x.copy()

basis.plot(x, color='k-').show()

dt = 1e-4
t = 0

for itr in range(10):

    uprev = basis.interpolate(x)

    @NonlinearForm
    def step(u, v, w):
        return ((u[0] - uprev) / dt * v[0]
                + u[0] * grad(u)[0] * v[0]
                + dd(u)[0][0] * v[0]
                + ddot(dd(u), dd(v)))

    for itr in range(50):

        J, f = step.assemble(x, basis)

        xprev = x.copy()
        x += solve(*condense(J, -f, D=basis.get_dofs()))

        res = np.linalg.norm(x - xprev)
        print(res)
        if res < 1e-8:
            break

    t += dt
    print('time: {}'.format(t))

ax = basis.plot(xinit, color='k-')
basis.plot(x, color='b-', ax=ax).show()