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).