Multiple boundary conditions and third derivative

[ThitiACHR]

Since my work relates to the calculation of gradient and curvature, the first and second derivative, respectively, I have to set the boundary conditions as follows:

1. p'(y=0) = n'(y=0) = 0 --> bndry_ydow = robin(0,1,0) for both variable.
2. p(y=a) = 0.01, n(y=a) = 0.1 --> bndry_yup = dirichlet(0.01) and dirichlet(0.1), respectively.
3. p'''(y=0=a) = n'''(y=0=a) = 0, which is the third derivative. This one is used to control the profile curvature behavior.

My questions are

1. Can we specify multiple boundary conditions at y = 0 and y = a? For example, at y = a, p(y=a) = 0.01 and p'''(y=a) = 0. If yes, how can I do that?
2. According to the third derivative (3rd bullet) of p and n, when I want to implement the boundary conditions, how do I do that? Should I create my own boundary condition in this case?

I am new to BOUT++. If anyone knows how to do it, please help. All answers are welcome.

[dschwoerer]

I do not think we have methods for setting the third derivative. I guess your best bet is to use costum boundary conditions.

See e.g. STORM or hermes-3 for examples of custom boundary conditions.

[bendudson]

Hi @ThitiACHR, thanks @dschwoerer !

BOUT++ sets boundary conditions by setting the value of the boundary cells: The boundary is at cell edges, half-way between cell centers.

The default methods are 2nd-order accurate, so only use one boundary cell. In that case only one boundary condition can be applied. To apply two boundary conditions on the same boundary (e.g on p and p''') there must be two boundary cells, which are only used for 4th-order finite difference methods.

e.g The final cell in Y inside the domain is at y = mesh->yend. We have cells like this:

... | yend-1 | yend |X| yend+1 | yend+2 |

where | marks cell edges, and |X| is the boundary location.
The value of p at yend-1 and yend is provided by the time integrator. The boundary condition needs to specify the values of p at yend+1 and yend+2.

See e.g.
https://github.com/boutproject/hermes/blob/master/hermes-1.cxx#L1758

This iterates over the boundary, ensuring that only the processors at the boundary apply the boundary condition:

for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) {
        for (int jz = 0; jz < mesh->ngz; jz++) {

Then sets each of the boundary conditions, starting with yend+1

for (int jy = mesh->yend + 1; jy < mesh->ngy; jy++) {

and sets the boundary condition using a mid-point method:

Ne(r.ind, jy, jz) = 2. * nesheath - Ne(r.ind, mesh->yend, jz);

That boundary condition is only applying one (Dirichlet) boundary condition, setting Ne = nesheath at the boundary. For your case you would have to calculate what the values of the two boundary cells should be, to impose your two boundary conditions.

[ThitiACHR]

Hi @dschwoerer and @bendudson

Thank you for your very informative reply. I would like to discuss some more about the model. Could you send me the Slack invitation link to my email: [email protected]? @bendudson

[ThitiACHR]

Hi @bendudson
Concerning your suggestion, that means I have to set the value of p and p'' in the boundary cells. For simplicity, I started with the lower-y boundary where boundary conditions are p'=0 and p'''=0.
This is my setting in the input file.

[mesh]
MYG = 2 # 2 guard cells on each boundary

[mesh:ddy] # using 4th order central finite different for DDY and D2DY2
first = C4
second = C4
...

[p]
bndry_ydown = none # prevent overwritten
bndry_yup = dirichlet_o4(0.01)
...

To impose the boundary conditions, I set the values in the boundary cells as shown below, to make
p and p'' at y=mesh->ystart-1 = y=mesh->ystart+1 and at y=mesh->ystart-2 = y=mesh->ystart+2,
so that p' and p''' will be 0 at y=mesh->ystart, respectively.
Note that my PDEs only evolve in the y-direction, so no need to loop in the x and z directions.

for(RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++){
       for(int jy = mesh->ystart - 1; jy >= 0; jy--){
        if(jy==mesh->ystart-1){
          p(0,jy,0) = p(0,mesh->ystart+1,0);
          D2DY2(p)(0,jy,0) = D2DY2(p)(0,mesh->ystart+1,0);
        }
        if(jy==mesh->ystart-2){
          p(0,jy,0) = p(0,mesh->ystart+2,0);
          D2DY2(p)(0,jy,0) = D2DY2(p)(0,mesh->ystart+2,0);
        }       
      }
    }

The result is shown below.
Note that the figure shows the zoom-in at the lower-y boundary of p, p', p'', and p''' as a function of the grid point index, where y=mesh->ystart locates at index = 2, indices 1 and 0 represent first and second boundary cells, respectively.

The simulation gives me such that p' = 0 at y=mesh->ystart as expected, but p''' is not 0. I think this is because the value of D2DY2(p) in the boundary cells is not followed as I set above.

This is what you meant. Do I understand correctly? If you have time, please point this out.

update 1
The result without boundary cells looks like this

[dschwoerer]

You need to manually do the taylor expansion around the point of the boundary, and then set the appropriate terms for the inversion, to then extrapolate into the boundary.

I tried to write it up here:
https://github.com/boutproject/BOUT-dev/pull/1179/files#diff-6de22069f1b200b5983d8abbb9347305e6f63724fecc39c98c6e5643e6aba746

That you could use as a starting point, to set all the desired quantities. Note that for setting third order derivatives, you need a relative high order scheme.

calling DDX/D2DX2 computes the derivate and gives you a field that contains the derivative. Of course you can change that field, but there is no code in BOUT++ to do the inverse and get back to the original field. So updating the original field via the operators will not work.

[ThitiACHR]

I want to make sure I understand correctly. There will be two things I have to be aware of.

1. When I do the Taylor expansion around the boundary point. The coefficient of the Taylor expansion and the number of neighbor points depend on the accuracy of the derivative scheme. For example, if I use the 4th-order accuracy central finite difference to calculate the third derivative at the boundary, I should expand 3 points on each side of ystart and yend. On one side, the information is provided by the time integrator, and on the other side (i.e., guard cells), I should use the interior points in the simulation domain to extrapolate the value into the guard cells so that the boundary conditions are met.

2. I should create the D3DY3 operator and use the corresponding order of accuracy (i.e., 4th-order in this case) to calculate and verify the third derivative at each boundary.

In this context what does the inversion mean?

Thank you in advance.

[dschwoerer]

I am not sure I get what you are trying to say.

Consider the following 1d grid. The numbers are the cells, the | are the boundaries:
0 1 | 2 3 4 5 6 7 | 8 9
^

Assume you want to write a BC for where the ^ is.
You take a taylor expansion around ^ of a given order n. Then you know your function at ^ for the first and 3rd derivative. Then you also take n-2 data points in your domain, where you know the value.

This gives you a linear system of equation, for your coefficients.
Assuming n is small enough you can solve this analytically (what I referred to as inverting)
Once you have your coefficients, you can put them into the taylor expansion to get the values for 0 and 1.

I have linked a script that does the inversion for you. It does not right now allow for setting 3rd derivatives, but should be easily extendable to do so.

I hope that helps.

Verifying the derivative is of course helpful :-)

[ThitiACHR]

I apologize for making you confused. This is what I meant

1. I want to set the BCs as follows:
   At LowerY p' = 0 and p'''= 0
   At UpperY p = 0.01 and p'''= 0

2. As you mentioned when setting the third derivative, I need a relatively high-order scheme. Suppose I use the 4th-order accuracy scheme. So the expression will look like this

$$\dfrac{dp^3}{dy^3} = \dfrac {1\cdot p(y-3h) - 8\cdot p(y-2h) + 13\cdot p(y-h) +0\cdot p(y) - 13\cdot p(y+h) + 8\cdot p(y+2h) - 1\cdot p(y+3h)} {8\cdot h^{3}}$$

where the values of p at y, y+h, y+2h, and y+3h are provided by time integrator.
As you can see we need 3 guard cells to calculate the third derivative, so the stencil should look like this (let's focus on the LowerY boundary)

| ystart-3 | ystart-2 | ystart-1 |x| ystart | ystart+1| ystart+2 | ystart+3 | ...
^

3. The goal is to find the value of p at ystart-3, ystart-2, and ystart-1 so the the BCs are met.
4. It seems like this is done by doing the Taylor expansion with order n (I choose n=3) around |x|, so that

$$p(y-|x|) = p(|x|) + \dfrac{p'(|x|)}{1!}(y-|x|) + \dfrac{p''(|x|)}{2!}(y-|x|)^{2} + \dfrac{p'''(|x|)}{3!}(y-|x|)^{3}$$

where we know that at |x| we have p' = 0 and p''' = 0, so the expression become

$$p(y-|x|) = p(|x|) + 0 + \dfrac{p''(|x|)}{2!}(y-|x|)^{2} + 0$$

I have seen in your stencils.md, you call p(|x|), p'(|x|), and so on as $f_{i}$ (coefficients).
This is what I understand right now. I still wonder how can we use this approach to obtain the value of the function (i.e. p in the guard cells). I want to make sure I understand your concept.
To build a linear system of equations, does it mean I have to do the Taylor expansion around points ystart-3, ystart-2, and ystart-1?