2D NLBVP

[Brian de Keijzer]

Hi there,

Firstly, thank you very much for having this software open-source.

I’m working on a two dimensional (k,t) NLBVP with NCC’s and periodic boundary conditions in k. See the equations below.

Is this possible to implement in Dedalus? Currently, I’m running into an error

 TypeError: __init__() missing 2 required positional arguments: 'layout' and 'func'

When adding the second equation to the problem. 

Best, 

Brian

[Brian de Keijzer]

The problem seems to arise due to the 'imag' and 'conj' general operators that I created. These do work for an IVP, though not for a NLBVP? The general operator code is found below. The error arises as soon as I add 'imag' or 'conj' in any equation, even for some test equation e.g. problem.add_equation("dy(fe) = conj(dy(E))")

def conj(field):

    # Call np conj function on the field's data

    return np.conj(field.data)

def conj_operator(field):

    # Return GeneralFunction instance that applies erf_func in grid space

    return de.operators.GeneralFunction(

        field.domain,

        layout = 'g',

        func = conj,

        args = (field,)

    )

de.operators.parseables['conj'] = conj_operator

def imag(field):

    # Call np imag function on the field's data

    return np.imag(field.data)

def imag_operator(field):

    # Return GeneralFunction instance that applies imag_func in grid space

    return de.operators.GeneralFunction(

        field.domain,

        layout = 'g',

        func = imag,

        args = (field,)

    )

de.operators.parseables['imag'] = imag_operator

Perhaps this error should move to a GitHub issue? Although I do still wonder if Dedalus is able to solve this problem. Any advice on that is greatly appreciated. 

Op donderdag 17 juni 2021 om 12:45:43 UTC+2 schreef Brian de Keijzer:

[Daniel Michael Lecoanet]

Hi,

Can you send the full error message?

Thanks,

Daniel

-- 
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/dcc41801-fc98-4a07-ba91-f389b9a50cd6n%40googlegroups.com.

[Brian de Keijzer]

Hi Daniel, 

Thank you for the quick response. Here is a minimal example including the error. 



---------------------------------- code ----------------------------------

import numpy as np

from dedalus import public as de

a=1

twidth = 25

# numerical params

kb, kt = (-np.pi/a, np.pi/a) # one Brillouin zone

tstart, tstop = (-twidth, twidth)

nk, nt = (64, 64)

# Create bases and domain

k_basis = de.Chebyshev('x', nk, interval=(kb, kt)) # x = k

t_basis = de.Chebyshev('y',nt, interval=(tstart, tstop)) # y = t , must stick to this notation otherwise the code wont work

domain = de.Domain([k_basis, t_basis], grid_dtype=np.complex128)

x,y = domain.grids(scales=1)

problem = de.NLBVP(domain, variables=['pk', 'fe', 'fh'], ncc_cutoff=1e-6)

def conj(field):

    # Call np conj function on the field's data

    return np.conj(field.data)

def conj_operator(field):

    # Return GeneralFunction instance that applies erf_func in grid space

    return de.operators.GeneralFunction(

        field.domain,

        layout = 'g',

        func = conj,

        args = (field,)

    )

de.operators.parseables['conj'] = conj_operator

def imag(field):

    # Call np imag function on the field's data

    return np.imag(field.data)

def imag_operator(field):

    # Return GeneralFunction instance that applies imag_func in grid space

    return de.operators.GeneralFunction(

        field.domain,

        layout = 'g',

        func = imag,

        args = (field,)

    )

de.operators.parseables['imag'] = imag_operator

problem.add_equation("dy(pk) = pk") # N.B. i = 1j;

problem.add_equation("dy(fe)= fe")

problem.add_equation("dy(fh)= conj(fh)") # may also replace 'conj' with 'imag'



---------------------------------- error (as a screenshot for readability) ----------------------------------


Switching from 

problem = de.NLBVP(domain, variables=['pk', 'fe', 'fh'], ncc_cutoff=1e-6)

to 

problem = de.IVP(domain, variables=['pk', 'fe', 'fh'], ncc_cutoff=1e-6)

gets rids of the error, but changes the problem of course.  


Best, 

Brian

Op donderdag 17 juni 2021 om 14:39:54 UTC+2 schreef Daniel Michael Lecoanet:

[Daniel Michael Lecoanet]

Because you’re trying to solve a nonlinear boundary value problem, I think Dedalus needs to know how to take symbolic derivatives of the conj and imag operators. Maybe Keaton or Ben or someone else with experience on NLBVPs can help you figure out how to do that.

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/e06c8973-38ae-47a7-b238-766b50729094n%40googlegroups.com.
<Screenshot_3.png>

[Geoffrey Vasil]

Alternatively, you might break the problem into real and imaginary parts as a workaround. -Geoff

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/D22FADFC-CFD6-4D5F-874D-C9A6E57A8B03%40northwestern.edu.

[Brian de Keijzer]

Hi Daniel and Geoffrey. 

Thank you both for your quick answers. Breaking the problem into real and imaginary parts brings me to the following equations.

 

Preferably the k-basis should be periodic, because the solution in k-space is periodic.  E.g. Fourier? However, then I cannot have NCC’s in k-space. Thus I need to use e.g. a Chebyshev basis, am I right? Using a Chebyshev basis for both the k and t mesh leads to a numerical unstable problem. The solutions explode rather quickly after a few perturbations. Moreover, I would like to impose periodic boundary conditions in k-space.

problem.add_bc("left(pR) - right(pR) = 0")

problem.add_bc("left(pI) - right(pI) = 0")

problem.add_bc("left(fe) - right(fe) = 0")

problem.add_bc("left(fh) - right(fh) = 0")

This causes a nan matrix as a solution. What am I missing? Relevant parts of the code are found below. Note that I substituted some variables k->x, t->y and d_k -> Dk. I'd attach a PDF but for some reason, I cannot find a way to do that.

It would be really great if this problem could be solved with Dedalus.

Best,

Brian

As for the commented out bc’s, left(variable)=right(variable) should be the actual boundary condition. Having either side (preferably the same) real number should also converge.

Initial conditions are all zero.

Op vrijdag 18 juni 2021 om 00:51:12 UTC+2 schreef Geoffrey Vasil:

[Keaton Burns]

Hi Brian,

It looks like you’re trying to do a 2D Chebyshev problem, which unfortunately Dedalus does not yet support for NLBVPs (it looks like you’re hitting some syntax errors regarding general functions first, but that is a less fundamental issue).  Dedalus also does not currently support 2D NLBVPs using Fourier / NCCs in Fourier directions, although we plan to add this capability soon. 

Best,

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/470f210d-a640-48dd-b6a7-942070b156f2n%40googlegroups.com.

[Geoffrey Vasil]

Hi Brian, 

My apologies. I would have given you a similar answer to Keaton if I had read the thread closer. 

There is one very annoying messy hack you could try. 

You could write the projection of your equations on to a Fourier basis by hand. And then you could programmaticly add the system by string manipulation in loops. 

The Chebyshev direction would still be Dedalus job. 

We actually did something like this to solve a problem on a network graph with Chebyshev problems on each edge. Because the problem is unstructured, there is no alternative in this case. We describe how to do this in the main methods paper.  There’s also an example scrip on GitHub. 

Doing this would be quite an effort. But it would be very rewarding if you got it working. 

What Keaton is saying is that we plan to have Dedalus do this automatically. But it may take a bit before we get some other stuff out of the way. 

Cheers, 

Geoff 

Sent from my iPhone

On 22 Jun 2021, at 5:03 am, Keaton Burns <keaton...@gmail.com> wrote:



Hi Brian,

It looks like you’re trying to do a 2D Chebyshev problem, which unfortunately Dedalus does not yet support for NLBVPs (it looks like you’re hitting some syntax errors regarding general functions first, but that is a less fundamental issue).  Dedalus also does not currently support 2D NLBVPs using Fourier / NCCs in Fourier directions, although we plan to add this capability soon. 

Best,

-Keaton

On June 18, 2021 at 3:23:49 PM, Brian de Keijzer (briande...@gmail.com) wrote:

Hi Daniel and Geoffrey. 

Thank you both for your quick answers. Breaking the problem into real and imaginary parts brings me to the following equations.

 

<1.png>

Preferably the k-basis should be periodic, because the solution in k-space is periodic.  E.g. Fourier? However, then I cannot have NCC’s in k-space. Thus I need to use e.g. a Chebyshev basis, am I right? Using a Chebyshev basis for both the k and t mesh leads to a numerical unstable problem. The solutions explode rather quickly after a few perturbations. Moreover, I would like to impose periodic boundary conditions in k-space.

problem.add_bc("left(pR) - right(pR) = 0")

problem.add_bc("left(pI) - right(pI) = 0")

problem.add_bc("left(fe) - right(fe) = 0")

problem.add_bc("left(fh) - right(fh) = 0")

This causes a nan matrix as a solution. What am I missing? Relevant parts of the code are found below. Note that I substituted some variables k->x, t->y and d_k -> Dk. I'd attach a PDF but for some reason, I cannot find a way to do that.

It would be really great if this problem could be solved with Dedalus.

Best,

Brian

<2.png>

<3.png>

<4.png>

<5.png>

As for the commented out bc’s, left(variable)=right(variable) should be the actual boundary condition. Having either side (preferably the same) real number should also converge.

<6.png>

Initial conditions are all zero.

<7.png>

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBgQwF3csWEfntpqKoiWu%3DjMP-tq3J9O9eaU8wiJ7isO2A%40mail.gmail.com.

[Brian de Keijzer]

Hi all, 

Thank you for the answers. Just a follow-up. A workaround that seems to work for this particular problem;

- Solve as an IVP

- Use a Fourier basis

- Use the equations of the PDE with conditions="nx!=0"

- Add boundary conditions, variables=0 with conditions="nx==0". Tiny bit unphysical, but it's only one mode at a very specific point, right? 

- Solve IVP with small enough timesteps. Each solution is a solution at timestep dt of the time dimension.

Nevertheless, I cannot wait for 2D Chebyshev support.

Best, 

Brian

Op maandag 21 juni 2021 om 23:59:40 UTC+2 schreef Geoffrey Vasil:

[MANISH KUMAR]

Hello Keaton,

  I want to solve a 2D periodic NLBVP problem in Dedalus. Is it possible to implement it in Dedalus? I plan to use the Fourier basis in the periodic direction and Chebyshev in the other direction in the Dedalus. However,  based on the online discussion I realized that Dedalus only supports 1D problems. Please, let me know if you modified the Dedalus to handle 2D NLBVPs or if there is a way to implement it.

Thanks,

Manish kumar

[Ben Brown]

Manish,

     Yes, the most recent version of dedalus (d3), now supports multidimensional NLBVPs.  These were implemented in PR 219 (https://github.com/DedalusProject/dedalus/pull/219).  I was actually trying one just as your e-mail came in.

Here's an example of how the d3 code looks:

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/e72ef34d-1a43-4136-82b9-0360057e782cn%40googlegroups.com.

[Ben Brown]

Sorry, sent to quickly.

The d3 code for the 2-D NLBVPs I'm playing with looks like this (only showing a sketch of relevant parts of it here):

coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=float)
xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(x0, x1), dealias=dealias)
zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(z0, z1), dealias=dealias)
x, z = dist.local_grids(xbasis, zbasis)

<...>

u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
p = dist.Field(name='p', bases=(xbasis, zbasis))

grad = lambda A: d3.Gradient(A, coords)
transpose = lambda A: d3.TransposeComponents(A)
e_ij = grad(u) + transpose(grad(u))

<...>

ω = (-div(skew(u)))
# Problem
problem = d3.NLBVP([u, p, τ1x, τ2x, τp], namespace=locals())
problem.add_equation('-1/Re*lap(u) + 1/η*u*ibm + grad(p) + lift(τ1x, -1) + lift(τ2x, -2) = -grad_p + skew(u)*ω')
problem.add_equation('div(u) + 1/Re*lift1(τ2x,-1)@ez + τp = 0')
problem.add_equation('integ(p) = 0')
if stress_free:
    problem.add_equation("ez@(ex@e_ij(z=z0)) = 0")
    problem.add_equation("ez@u(z=z0) = 0")
    problem.add_equation("ez@(ex@e_ij(z=z1)) = 0")
    problem.add_equation("ez@u(z=z1) = 0")
else:
    problem.add_equation('u(z=z0) = 0')
    problem.add_equation('u(z=z1) = 0')

<...>

solver = problem.build_solver(ncc_cutoff=5e-3)
pert_norm = np.inf
tol = 1e-4
while pert_norm > tol:
    solver.newton_iteration()

Here I'm trying to solve for the steady-state Poiseuille flow around an object in the flow (represented by the 2-d NCC field "ibm", which here is a circle using volume penalisation approaches from Hester et al 2021 JCP; ibm is the Γ term in their eqn 9).  

The nonlinearity is in the u.grad(u) term, which here is "skew(u)*ω" (this is cross(u, ω), with ω the vorticity).

I've done LBVPs of this, and they seem like they are working correctly.  The NLBVP is a fair bit slower than the LBVP (~10x slower per iteration on this particular run) and so is still converging, but it looks like it is also latching onto a valid solution.  I'll send an example image of the LBVP and NLBVP if it finishes converging in the next bit.

--Ben

[MANISH KUMAR]

Hello Ben,

 Thanks for your reply. I would love to see the result and the complete code for this example.

Best regards,

Manish Kumar 

[Ben Brown]

Dear Manish,

     I can confirm that 2-D LBVPs and NLBVPs appear to be working correctly in d3.

An example 2-D LBVP looks like this:

The problem here is pressure driven flow past a cylinder, with the cylinder (grey region at center) represented by an immersed boundary method.  The lines show the flow and the colors show the vorticity of the flow.  Since it's a linear BVP here, all changes are local to the cylinder.

The corresponding NLBVP looks like this (at Reynolds number Re=5) :

As expected, the bulk flow now stretches the vorticity and forms a wake, which extends beyond the object.  This particular problem is in a horizontally periodic domain, so the cylinder encounters its own wake again (wrapping from +4  to -4), so it's more analogous to an infinite array of such cylinders.  All of this lines up with my physical intuition for this problem, and when I've changed parameters, I've seen things changing in expected ways.  IVPs of the same problem produce similar solutions.

In your reply, you asked for a complete code example.  I won't be providing that here; this particular work is still in progress, and it's not appropriate to publicly release it at this time.

Generally, the core dev team works hard to release our tools as we publish scientific papers to provide concrete examples to help the community, in addition to the examples we develop and release with the code itself.  It is a bigger (and potentially inappropriate) ask to request complete code examples of in-development science.  There are many reasons for this.

You've got a sketch of the core NLBVP code in the e-mail chain above, and you should be able to work out a fully functional example from that, paired with our existing examples.

Best,  

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/a746bdf7-744e-4f66-9c59-d08caedc0524n%40googlegroups.com.

[MANISH KUMAR]

Hello Ben, 

  LBVP solver assumes linear LHS and known RHS. You mentioned that you also used LBVP to solve for steady-state Poiseuille flow around an object. Is it possible for you to mention how you linearized u.grad(u) term in the NS equations while using LBVP? I am getting a memory issue with NLBVP.  Therefore, I would like to try LBVP for my problem.

Thanks,

Manish Kumar

On Thursday, September 22, 2022 at 6:06:33 PM UTC-5 Ben Brown wrote:

[Ben Brown]

Manish,

     In the images I sent in this thread, I linearized u.grad(u) by setting u.grad(u) = 0; I was testing functionality rather than trying to use it to solve the problem.

If you know a base u_0 to work from (maybe you take poiseulle flow everywhere, or uniform flow, or something that matches the boundary conditions in the absence of the object-in-flow), then you could have a u.grad(u_0) + u_0.grad(u) term on the LHS, where u_0 is a field object and an NCC that you specify.  Again, I didn't do that here, because I was mostly playing with implementation and confirming that the NLBVP gave a physical and different answer than the LBVP.

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/86d4b7a4-2e77-4cb8-a684-d042a8d78161n%40googlegroups.com.