Resizing guard celss

[mredenti]

In the boutdata.restart.resize method one can pass in the number of guard cells as argument

boutdata/boutdata/restart.py

Lines 107 to 108 in 1cb6d08

	mxg, myg : int, optional
	Number of ghost points in x, y (default: 2)

but nothing is done with them in the remaining of the code. It seems that you are resizing guard cells but guard cells should simply be pointing to adjacent grids edges and this not done in the code. What happens when we restart the simulations?

To me, it seems the safest route using this method is to resize a single restart file. In other words, one has to first create a unique restart file from the simulation (BOUT.restart.0.nc) and resize that. Multiple restart files can be created after the resizing method if needed. You should add a warning to users.

[johnomotani]

Yes, this does seem wrong. The 'new coordinates' ought to be rescaled based on the interval than includes the grid points (not guard points) and 1/2 a grid spacing on either side.

Also, this part is not needed any more as nz can have any value (although products of small primes are more efficient for FFTs).

boutdata/boutdata/restart.py

Lines 148 to 154 in 1cb6d08

	def is_pow2(x):
	"""Returns true if x is a power of 2"""
	return (x > 0) and ((x & (x - 1)) == 0)
	if not is_pow2(newNz):
	print("ERROR: New Z size {} must be a power of 2".format(newNz))
	return False

[mredenti]

Yes, this does seem wrong. The 'new coordinates' ought to be rescaled based on the interval than includes the grid points (not guard points) and 1/2 a grid spacing on either side.

I agree with the statement " The 'new coordinates' ought to be rescaled based on the interval than includes the grid points (not guard points)" but I do not understand "and 1/2 a grid spacing on either side"

[johnomotani]

Notionally the grid points are in the centre of a cell, which touches the grid on a neighbouring processor (or a physical boundary) at a 'cell face' half way between the adjacent grid points.

[mredenti]

Notionally the grid points are in the centre of a cell, which touches the grid on a neighbouring processor (or a physical boundary) at a 'cell face' half way between the adjacent grid points.

I see, but this is not an issue if I have a single BOUT.restart.0.nc file (1 processor for the all simulation) right?

[johnomotani]

I see, but this is not an issue if I have a single BOUT.restart.0.nc file (1 processor for the all simulation) right?

Even with one file it's the same thing. Boundary conditions are applied at 'cell faces' so to keep the physical size of your simulation the same, you need to account for the half-grid-spacing distance between the first/last grid point and the 'boundary'.

[mredenti]

Even with one file it's the same thing. Boundary conditions are applied at 'cell faces' so to keep the physical size of your simulation the same, you need to account for the half-grid-spacing distance between the first/last grid point and the 'boundary'.

I think though this is not an issue with the "nearest" interpolation algorithm as it does not use the co-ordinates? The "linear" method does not work, it yield NaN values.

I will try fix this. But what about if I have two guard cells in the x-direction. I assume the second guard cell is now at the next cell face, i.e. a distance dx away from the location of the first guard cell or 2dx?

[johnomotani]

Let me try to ASCII-art this...
grid points are x
guard cell points are o
cell faces are |
the physical domain boundary is I

o | o I x | x | x | ... | x | x | x I o | o

The distance between adjacent grid points x is dx. The distance between guard points o is also dx, and between adjacent grid points x and guard points o is likewise dx.

Zooming in near the physical boundary, just to show that the distance between each grid/guard point and the nearest cell face is 0.5*dx

o <-0.5*dx-> | <-0.5*dx-> o <-0.5*dx-> I <-0.5*dx-> x <-0.5*dx-> | <-0.5*dx-> x

The physical x-coordinate runs between the 'physical domain boundaries' I. So for example if you interpolate to higher resolution, the domain boundaries I remain in the same place and the left-most grid point x on the higher resolution grid has to be a bit to the left of the left-most grid point x on the lower resolution grid. Similarly the right-most grid point x on the higher resolution grid has to be a bit to the right of the right-most grid point x on the lower resolution grid.

Edited to add: maybe it would be helpful to add, if we took some normalised coordinate X covering the whole domain, we could say X=0 at the left physical boundary and X=1 at the right physical boundary. Then the left-most grid point would be at X=0.5/nx, the next point at X=1.5/nx, etc. The guard points would correspondingly be at X<0 or X>1.

[mredenti]

Just to clarify, is the physical boundary domain I an actual grid point or not?

Also suppose I specify Nz=4 and dz=1. BOUT++ outputs the z-coordinates as 0, 1, 2, 3 which I guess are equivalent to outputting 0.5, 1.5, 2.5, 3.5 since it is just a translation. However, if I double the resolution Nz=8 and dz=0.5 then the z-coordinates are outputted as 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5 whereas they are actually located at 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75 This is a problem because the relative co-ordinate location between the two resolutions is wrong. The shift/translation is by 0.5 for Nz=4 but by 0.25 for Nz=8
I am not sure the resizing as it is implemented takes into account the correct relative co-ordinate location when doing the interpolation

[johnomotani]

Just to clarify, is the physical boundary domain I an actual grid point or not?

No I is not a grid point.

The z-coordinate is special case as it is always (at least for the moment...) periodic. What did you mean by 'BOUT++ outputs the z-coordinate? In the BOUT.dmp.*.ncfiles there are no x/y/z coordinate values. xBOUT does construct some coordinate values... However, I don't think there is a consistent definition of whether 'the z-coordinate' starts at 0 on the grid point at index zero or at the cell face dz/2 away from that, as we've never (as far as I know) needed one. My vote would be for the coordinate to start on the cell face, for consistency, but that is a personal opinion. The resizing as implemented inboutdata.restart` probably does not take into account any of what I have said about coordinates. It was written 6 years ago and doesn't seem to have been touched since...

[mredenti]

The z-coordinate is special case as it is always (at least for the moment...) periodic. What did you mean by 'BOUT++ outputs the z-coordinate? In the BOUT.dmp.*.nc`files there are no x/y/z coordinate values. xBOUT does construct some coordinate values...

Yeah, sorry. I meant xBOUT.

However, I don't think there is a consistent definition of whether 'the z-coordinate' starts at 0 on the grid point at index zero or at the cell face dz/2 away from that, as we've never (as far as I know) needed one. My vote would be for the coordinate to start on the cell face, for consistency, but that is a personal opinion.

I would also think it does not matter when looking at a single resolution simulation. But if I am for example doubling the resolution we should preserve the relative position of the grid points between the coarse and the fine resolution.

What about the x-coordinate instead? The co-ordinate is at the cell centre correct?

[johnomotani]

I would also think it does not matter when looking at a single resolution simulation. But if I am for example doubling the resolution we should preserve the relative position of the grid points between the coarse and the fine resolution.

If you're implementing it can do what you want it to do, just needs documenting what the function does. If someone wants a different implementation in future then they can add an option...

What about the x-coordinate instead? The co-ordinate is at the cell centre correct?

What do you mean? A coordinate is not a thing at a point, it's a continuous thing (not sure what word to use instead of 'thing' - 'property of the space' maybe?).

[mredenti]

What about the x-coordinate instead? The co-ordinate is at the cell centre correct?

What do you mean? A coordinate is not a thing at a point, it's a continuous thing (not sure what word to use instead of 'thing' - 'property of the space' maybe?).

Sorry I meant the x grid points. I would like to understand if I double the resolution where are the new grid points placed relative to the old grid points. I would think it would be something like this:

Suppose I specify Nx=4 + 2MXG and dx=1. The grid points are located at 0.5, 1.5, 2.5, 3.5. If I double the resolution Nx=8 + 2MXG and dx=0.5 then the grid points are located at 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75?

[johnomotani]

Right, yes the grid points are at cell centres (variables can also be defined on staggered grids in BOUT++ but it's fine to ignore that - xBOUT still doesn't have any support for staggered grids).

Your example of x-coordinate values of grid points when changing Nx looks correct to me.

[mredenti]

The z-coordinate is special case as it is always (at least for the moment...) periodic. What did you mean by 'BOUT++ outputs the z-coordinate? In the BOUT.dmp.*.ncfiles there are no x/y/z coordinate values. xBOUT does construct some coordinate values... However, I don't think there is a consistent definition of whether 'the z-coordinate' starts at 0 on the grid point at index zero or at the cell face dz/2 away from that, as we've never (as far as I know) needed one. My vote would be for the coordinate to start on the cell face, for consistency, but that is a personal opinion. The resizing as implemented inboutdata.restart` probably does not take into account any of what I have said about coordinates. It was written 6 years ago and doesn't seem to have been touched since...

But BOUT must create a set of z-values in order to evaluate an initial condition that we specify in BOUT.inp. For example, if I want to initialise a variable as mixmode(z) I know z \in [0, 2pi] but which points are picked in this interval? If I say nz=4, I suppose it takes nz values equally spaced between 0 and 2pi including 0 but not 2pi (since it's periodic) or else?

[johnomotani]

But BOUT must create a set of z-values in order to evaluate an initial condition that we specify in BOUT.inp.

Ah, the expressions in BOUT.inp are a slightly different thing again. The x, y, and z as evaluated in the expressions are not defined according to the actual simulation coordinates (as defined by dx, dy, dz). Instead x always goes between 0 and 1 and z always goes between 0 and 2pi. y is defined in different ways depending on the topology of the magnetic equilibrium.

To answer your question anyway: yes, as you suspected the z values resulting from BOUT.inp expressions will take nz equally spaced values between 0 and 2pi, including 0 and excluding 2pi (2pi would be at nz+1).

[mredenti]

But BOUT must create a set of z-values in order to evaluate an initial condition that we specify in BOUT.inp.

Ah, the expressions in BOUT.inp are a slightly different thing again. The x, y, and z as evaluated in the expressions are not defined according to the actual simulation coordinates (as defined by dx, dy, dz). Instead x always goes between 0 and 1 and z always goes between 0 and 2pi. y is defined in different ways depending on the topology of the magnetic equilibrium.

To answer your question anyway: yes, as you suspected the z values resulting from BOUT.inp expressions will take nz equally spaced values between 0 and 2pi, including 0 and excluding 2pi (2pi would be at nz+1).

Right, whereas for x not-periodic the function is evaluated at both x=0 and x=1? Or slightly to the right and left respectively if the x-grid points are cell-centred? (I could not find where is this in the code)

[johnomotani]

Right, whereas for x not-periodic the function is evaluated at both x=0 and x=1? Or slightly to the right and left respectively if the x-grid points are cell-centred? (I could not find where is this in the code)

Unfortunately I think the documentation is actually wrong at the moment (going to make a PR to fix this now!)... The x evaluated in input file expressions will be (by default, for CELL_CENTRE variables) 0.5/(nx-2*MXG) on the first (non-boundary) grid point, and (nx-2*MXG-0.5)/(nx-2*MXG) on the last (non-boundary) grid point.

[johnomotani]

The x evaluated in input file expressions will be (by default, for CELL_CENTRE variables) 0.5/(nx-2MXG) on the first (non-boundary) grid point, and (nx-2MXG-0.5)/(nx-2*MXG) on the last (non-boundary) grid point.

...hang on, I was checking in order to update the documentation and I'm not sure this is actually true...

[johnomotani]

..hang on, I was checking in order to update the documentation and I'm not sure this is actually true...

What I wrote above was correct, panic over! I confused myself looking at the output of a test simulation because the boundary points were set by a Neumann boundary condition, so did not contain the negative or >1 values that I'd expected...

[mredenti]

Right, whereas for x not-periodic the function is evaluated at both x=0 and x=1? Or slightly to the right and left respectively if the x-grid points are cell-centred? (I could not find where is this in the code)

Unfortunately I think the documentation is actually wrong at the moment (going to make a PR to fix this now!)... The x evaluated in input file expressions will be (by default, for CELL_CENTRE variables) 0.5/(nx-2*MXG) on the first (non-boundary) grid point, and (nx-2*MXG-0.5)/(nx-2*MXG) on the last (non-boundary) grid point.

Is there any chance you could point me to where it does this in the code? I would like to convince myself of this too and in particular how BOUT treats non-periodic x and z differently

[mredenti]

..hang on, I was checking in order to update the documentation and I'm not sure this is actually true...

What I wrote above was correct, panic over! I confused myself looking at the output of a test simulation because the boundary points were set by a Neumann boundary condition, so did not contain the negative or >1 values that I'd expected...

If we have 0 Dirichlet boundary conditions and two guard cells in the x-direction, the 1st guard cell left of the boundary would take the negative of the value at the first grid point right of the boundary whereas the second guard cell would take value of zero. I am wondering whether when importing the data with XBOUT I should perhaps drop the guard cells and append zero values to left and right of the grid that would account for the zero boundary. I am thinking this is more correct when doing interpolation (xarray datasets have a .interp method that calls the scipy interpolation)

[johnomotani]

With 0-value Dirichlet boundary conditions, the first boundary point will be as you say the negative of the first grid point. The second boundary point though should be an extrapolated value (using cubic polynomial interpolation from the first boundary cell and first 2 grid cells). Keeping those boundary cells should be the most consistent/accurate thing to do. Setting the boundary cells to 0 in the input would result in the first one or two grid points having values larger (assuming the input values are positive on the grid points) than they should be in the output - the output would be (more or less) as if you imposed a Dirichlet boundary condition at a location a bit to the left of the actual left boundary (i.e. at the location of the first boundary cell in the input grid).

[mredenti]

With 0-value Dirichlet boundary conditions, the first boundary point will be as you say the negative of the first grid point. The second boundary point though should be an extrapolated value (using cubic polynomial interpolation from the first boundary cell and first 2 grid cells). Keeping those boundary cells should be the most consistent/accurate thing to do. Setting the boundary cells to 0 in the input would result in the first one or two grid points having values larger (assuming the input values are positive on the grid points) than they should be in the output - the output would be (more or less) as if you imposed a Dirichlet boundary condition at a location a bit to the left of the actual left boundary (i.e. at the location of the first boundary cell in the input grid).

I see so I should be fine by keeping the boundary/guard cells and every grid/guard point is a distance dx away from one another as you described earlier. Since there no actual boundary points (physical boundary domain) there is no need to specify 'dx/2' at any of the points

[johnomotani]

I see so I should be fine by keeping the boundary/guard cells and every grid/guard point is a distance dx away from one another as you described earlier. Since there no actual boundary points (physical boundary domain) there is no need to specify 'dx/2' at any of the points

No, just need to have a consistent definition of the coordinate values at the grid points (as in your comment above #78 (comment)).