I've spotted what I believe are two problems with the stretching algorithm in subroutine grd7.
The first has to do with how grd7 defines the spacing of grid points when the z stretching factor zsf is set to 0.
If zsf is set to zero then presumably the idea is to have no stretching and thus a uniform grid across the entire z axis.
I have not been getting this. The reason I know this is because I included a check on the grid spacing in my Fourier space plotting routine to safeguard against non-uniform sampling intervals - since DFT's only make sense when the sampling intervals are uniform (one should use a "Z-transform" if spacing is non-uniform).
Recently I started going line-by-line again through the entire code to shore up my understanding of it and to see if I can track down the causes of some of the unexpected things I've been seeing. Yesterday I I spotted the problem in grd7 but wanted to test things first before writing an issue again.
What I find is that the algorithm as written fails to safeguard against two special cases. The first case is that of complete uniformity. It turns out that even if I set zsf to 0 in inputs.dat I still get a non-negligible spacing difference between interior spaces and the spaces between the first and second grid points. I haven't tracked down the reason for this (whatever it is, it's subtle) , but as a test I added a conditional to grd7 which does the following:
if (zsf .gt. 0) then
…evaluate dz_grid and qz in the usual way
else
do k = 1, nz
qz(k) = k * delz
dz_grid(k) = delz
enddo
endif
If i do this and run a test simulation and then subsequently put the results through my Fourier color map function, I no longer get the complaints from the function about uneven spacing.
I'm trying to determine if this issue could shed light on why I've been having to include that factor of two in my Alfvén speed calculation. The jury is still out on that.
The other special case where a problem arises is when there is a mismatch between the choices for the "grad" and "height_stretch" parameters on the one hand, and the size of the domain in z (i.e. the value of nz) on the other.
Here's what I'm seeing. the quantities "grad" and "height_stretch" which are input via inputs.dat together determine the value of the integer nrgrd which is used to delimit the "up from center" and "down from center" do loops. In each case there are two loops for initializing dz_grid. The first loops over the portion of dz_grid starting from the "center" element - i.e. where k = rk = nz/2 - out to where stretching is to begin as determined by the sum rk + nrgrd.
The problem is that - depending on how nz, grad, and height_stretch are chosen - you can run into a situation where the limits on the do loops don't make any sense, and it may be difficult or impossible to predict the outcome.
I'll give you an example from my own experiments.
The copy of the code I initially cloned had the values 3000 and 1.5 set for grad and height_stretch respectively. This leads to a value of 4500 for nrgrd.
Meanwhile my value of nz was set to 2001. So rk is 1000.
Hence the limits on k for the first "up from center" do loop are k = 1000, 5500
and for the second they are k = 5501, 2001
The problem with this is that dz_grid is dimensioned with the value of nz which is
So the first of these loops goes out of the bounds of dz_grid and overwrites whatever is in those memory addresses. Maybe this causes trouble, maybe it doesn't but it's risky.
If I remember my fortran correctly the second loop won't execute because the first index is larger than the second so that should not be a problem.
Meanwhile in the "down from center" portion, the first loop has limits k = -3500,999
and as far as I know this should execute so once again there is some writing of data beyond the proper bounds of dz_grid.
The second loop has limits k = 1, -3501 so once again this loop should not execute.
Maybe these lines cause trouble or maybe though don't but they give me the willies and a simple check ensuring that nrgrd doesn't go beyond some maximum possible value given the size of the grid could prevent the loops from any overwriting. So I thought I'd better point this out.
I've spotted what I believe are two problems with the stretching algorithm in subroutine grd7.
The first has to do with how grd7 defines the spacing of grid points when the z stretching factor zsf is set to 0.
If zsf is set to zero then presumably the idea is to have no stretching and thus a uniform grid across the entire z axis.
I have not been getting this. The reason I know this is because I included a check on the grid spacing in my Fourier space plotting routine to safeguard against non-uniform sampling intervals - since DFT's only make sense when the sampling intervals are uniform (one should use a "Z-transform" if spacing is non-uniform).
Recently I started going line-by-line again through the entire code to shore up my understanding of it and to see if I can track down the causes of some of the unexpected things I've been seeing. Yesterday I I spotted the problem in grd7 but wanted to test things first before writing an issue again.
What I find is that the algorithm as written fails to safeguard against two special cases. The first case is that of complete uniformity. It turns out that even if I set zsf to 0 in inputs.dat I still get a non-negligible spacing difference between interior spaces and the spaces between the first and second grid points. I haven't tracked down the reason for this (whatever it is, it's subtle) , but as a test I added a conditional to grd7 which does the following:
if (zsf .gt. 0) then
else
enddo
endif
If i do this and run a test simulation and then subsequently put the results through my Fourier color map function, I no longer get the complaints from the function about uneven spacing.
I'm trying to determine if this issue could shed light on why I've been having to include that factor of two in my Alfvén speed calculation. The jury is still out on that.
The other special case where a problem arises is when there is a mismatch between the choices for the "grad" and "height_stretch" parameters on the one hand, and the size of the domain in z (i.e. the value of nz) on the other.
Here's what I'm seeing. the quantities "grad" and "height_stretch" which are input via inputs.dat together determine the value of the integer nrgrd which is used to delimit the "up from center" and "down from center" do loops. In each case there are two loops for initializing dz_grid. The first loops over the portion of dz_grid starting from the "center" element - i.e. where k = rk = nz/2 - out to where stretching is to begin as determined by the sum rk + nrgrd.
The problem is that - depending on how nz, grad, and height_stretch are chosen - you can run into a situation where the limits on the do loops don't make any sense, and it may be difficult or impossible to predict the outcome.
I'll give you an example from my own experiments.
The copy of the code I initially cloned had the values 3000 and 1.5 set for grad and height_stretch respectively. This leads to a value of 4500 for nrgrd.
Meanwhile my value of nz was set to 2001. So rk is 1000.
Hence the limits on k for the first "up from center" do loop are k = 1000, 5500 and for the second they are k = 5501, 2001
The problem with this is that dz_grid is dimensioned with the value of nz which is
If I remember my fortran correctly the second loop won't execute because the first index is larger than the second so that should not be a problem.
Meanwhile in the "down from center" portion, the first loop has limits k = -3500,999 and as far as I know this should execute so once again there is some writing of data beyond the proper bounds of dz_grid.
The second loop has limits k = 1, -3501 so once again this loop should not execute.
Maybe these lines cause trouble or maybe though don't but they give me the willies and a simple check ensuring that nrgrd doesn't go beyond some maximum possible value given the size of the grid could prevent the loops from any overwriting. So I thought I'd better point this out.