Hypnotoad: double precision and y-boundary guard cells by johnomotani · Pull Request #1647 · boutproject/BOUT-dev

johnomotani · 2019-03-20T22:56:18Z

Adds the option to write y-boundary guard cells from Hypnotoad. This includes the 'upper' target in double null configurations; the number of y-points in the gridfile may therefore be ny+2*n_yguards in single-null or ny+4*n_yguards in double null configurations. Requires the updates to GridFile in #1646 to read these grid files correctly if n_yguards != 0, so the number of y-guards defaults to zero.

Upgrades Hypnotoad to use double-precision arithmetic everywhere. The motivation for this is that the change to number of y-boundary guard cells changes some of the indexing in Hypnotoad, which effectively chages some of the rounding errors. I generated one grid file with a version of Hypnotoad with double-precision but without the y-boundary changes here and another with this version: the relative errors between the two were mostly below 1.e-6, but bxcvx and ShiftTorsion were ~3e-6. That was for an 'orthogonal' grid - field-aligned grids were significantly worse because of the integrated shear appearing in metric components. Comparing the single-precision version in next to the double precision (with or without the y-boundary changes) the relative errors are up to 1.8 in the metric components!

Also a few other fixes:

turn on refinement of the starting contour for non-orthogonal grids in all cases
pass through the simple setting if create_grid is re-run
rename gen_surface to gen_surface_hypnotoad so it doesn't get masked by the version from idllib (only matters because gen_surface_hypnotoad needed changing in this PR)
H was calculated as the y-derivative of qinty, but qinty is a y-integral, so can just set H from the integrand
calculate curvature with I instead of sinty so that when writing out x-z orthogonal (non-field-aligned) grids the curvature does not include the integrated shear (the option to write x-z orthogonal output sets I=0)
fix a place where nnpol was re-calculated wrongly. This did not cause a problem because the value would always be too big, and only had the effect of making a few output arrays too long (but they were just filled with harmless zeros).
don't smooth beta - not sure this is a fix as such, but smoothing propagates differences from the boundary points further into the grid and made it harder to compare grids with/without y-boundary guard cells - hopefully with double precision beta won't be too noisy.

Previously, for nonorthogonal grids the starting positions (points along the separatrix) were refined only when the grid was being created in the core. This commit moves the refinement to a position (previously commented-out) in grid_region_nonorth which ensures it is done in every case. This makes non-orthogonal grid generation more robust, especially around x-points.

When create_grid fails, sometimes it changes some options and tries calling create_grid again. Previously the 'simple' setting was not being passed through to the new calls.

Previously Hypnotoad would use the version from idllib if that was in the $IDL_PATH.

The quantity 'H' is the y-derivative of qinty, but qinty is an integral in y. Instead of calculating numerically the integral and derivative, just set 'H' equal to the integrand. This removes errors at the y-boundaries.

When calculating 'bxcvz' the integrated shear 'sinty' was used. This commit updates to use 'I' instead, so that when 'I' is set to zero for orthogonal coordinates the curvature is calculated consistently.

Calculate and save values at user-chosen number of guard cells beyond the end of the grid at the divertor targets. Default value is 0 to ensure compatibility with older versions of BOUT++ that don't check the (newly added) 'y_boundary_guards' value in grid files.

Change my_int_y to make its result zero at the first grid point, i.e. subtract the contribution of the lower y-boundary guard cells (which then have negative values in the result). Make qinty, sinty consistent on open surfaces, by setting them to zero at the first y-grid point in the domain for open flux-surface regions. This makes the results the same when the number of y-boundary guard cells is changed.

1.1.4

Previously Hypnotoad used single precision floats in most places. Using double precision improves reproducibility of grids by reducing rouding errors. In particular rounding errors due to limited precision (e.g. in interpolations) are amplified by numerical derivatives; there are some quantities, e.g. zShift (aka qinty) that require interpolation, then differentiation, then integration so that errors can accumulate quite badly. To ensure double precision: - pass /DOUBLE argument to all calls to INTERPOLATE - use double precision for constants, by entering them with a 'D', e.g. 1.5D*f... - use DBLARR instead of FLTARR to create arrays

1.2.0

Makes results more consistent with/without y-boundary guard cells.

bendudson · 2019-03-24T22:54:47Z

Thanks @johnomotani this looks good. My one worry is what happens if a new grid with n_yguards != 0 is read into an older version of the code. Is there a way to make it fail rather than reading incorrectly? Is there a way to catch this with version number checking? (I guess not for older versions of BOUT++).

johnomotani · 2019-03-25T10:10:57Z

@bendudson I think the n_yguards != 0 case should be mostly caught even by old versions of BOUT++ (hooray for input checking!):

BOUT-dev/src/mesh/data/gridfromfile.cxx

Lines 200 to 210 in 1c8d769

    
           ///Check if field dimensions are correct. y-direction 
        
           if (field_dimensions[1] == m->GlobalNy) { ///including ghostpoints 
        
             ny_to_read = m->LocalNy; 
        
             yd = 0; 
        
           } else if ( field_dimensions[1] == m->GlobalNy - 2*myg ) {///including ghostpoints 
        
             ny_to_read = m->LocalNy - 2*myg; 
        
             yd = myg; 
        
           } else { 
        
             throw BoutException("Could not read '%s' from file: y-dimension = %i do neither match ny = %i" 
        
                         "nor ny-2*myg = %i ", name.c_str(), field_dimensions[1], m->GlobalNy, m->GlobalNy-2*myg); 
        
           }

Single null cases with 0 or MYG y-boundary guards will read correctly, double null cases with MYG y-boundary guards will throw an exception (because field_dimensions[1] > GlobalNy). The only plausible case that would be read incorrectly is a double-null with 1 y-boundary guard cell when the simulation has MYG=2, which would slip into the first case.

I can't think of a way to check more strictly? I meant to make people think before using grids with n_yguards > 0 with the comment about backward compatibility here

BOUT-dev/tools/tokamak_grids/gridgen/hypnotoad.pro

Lines 1072 to 1075 in 3049f4a

    
           l = WIDGET_LABEL(tab1, value = '(default 0 for backward compatibility,' + STRING(10B) $ 
        
                                          + 'recommended to set to number of' + STRING(10B)      $ 
        
                                          + 'y-guards in your simulation, e.g. 2)',              $ 
        
                                          /ALIGN_LEFT)

bendudson · 2019-03-25T10:13:34Z

tools/tokamak_grids/gridgen/adjust_jpar.pro

    SURFACE, data.jpar0, tit="Input Jpar0", chars=2
    SURFACE, jpar, tit="New Jpar0", chars=2
-    PLOT, data.jpar0[0,*], tit="jpar at x=0. Solid=input", yr=[MIN([data.jpar0[0,*],jpar[0,*]]), $
+    PLOT, data.jpar0[0,*], tit="jpar at x=0.D Solid=input", yr=[MIN([data.jpar0[0,*],jpar[0,*]]), $


This 0.D probably not needed, though doesn't hurt

johnomotani · 2019-03-25T10:16:49Z

Actually, I just noticed that BoutMesh tries to read MXG and MYG from the grid file before reading from options. So if Hypnotoad wrote MYG = n_yguards into the gridfile, then the check at l.200-210 should ensure that an exception is thrown if the grid would not be read correctly.

Don't need double precision specification in a plot label.

bendudson · 2019-03-25T10:22:42Z

tools/tokamak_grids/gridgen/create_nonorthogonal.pro

-      start_ri[i] = ri1
-      start_zi[i] = zi1
-    ENDFOR
+    ;FOR i=0, np-1 DO BEGIN


Is this refinement not needed, or does it lead to incorrect results?

I thought it was not needed, because I switched on refinement in grid_region_nonorth. Looking again, they're not quite the same thing, because the refinement I switched on is of the line after interpolating onto the poloidal grid, whereas this was for the line that is input to the interpolation. I guess more refinement should never hurt and might help, so I'll put back this stage of refinement (possibly at the beginning of grid_region_nonorth if that reduces code duplication.

bendudson · 2019-03-25T10:23:47Z

tools/tokamak_grids/gridgen/create_nonorthogonal.pro

@@ -1450,20 +1516,20 @@ FUNCTION create_nonorthogonal, F, R, Z, in_settings, critical=critical, $
      ;;   REPEAT BEGIN


I guess all this commented out code can be removed?

bendudson · 2019-03-25T10:29:30Z

tools/tokamak_grids/gridgen/leg_separatrix.pro


  ; Use finite-differencing
-  ;di = 2.
+  ;di = 2.DD


DD? Anyway this commented-out code can go I think.

bendudson

Thanks @johnomotani that's a lot of changes!

Previously was done before calling grid_region_nonorth in the closed-field line only case (though this was recently commented out), but not in other cases. Now done for all cases.

If y_boundary_guards>0, then write 'MYG=y_boundary_guards' into the grid file. This ensures that checks present in BOUT++ versions earlier than 4.3 will only allow compatible grids to be read; i.e. single null with y_boundary_guards=0 or MYG, or double null with y_boundary_guards=0.

johnomotani · 2019-03-25T22:10:00Z

With the last commit I pushed, Hypnotoad writes MYG to grid files. I've checked that this makes BOUT-4.2.2 throw an exception where with grid files generated before that commit (no MYG) they could be read incorrectly (if a strange value of MYG was chosen). So Hypnotoad grid files should now be backwards-safe, although only backwards compatible for single null or if y_boundary_guards=0.

johnomotani · 2019-03-25T22:13:53Z

@bendudson: In ca936c6 I turned back on the grid refinement that I commented out before, but moved into grid_region_nonorth so it is used in every case. The grids before and after look similar and sensible (tested on a MAST disconnected-double-null case), but the grid points move by up to ~2cm. Large changes in metrics, etc. are mostly just around the X-point where it's not surprising that small differences in position would cause large changes in value. So I think the update is OK, presumably the new grids are slightly 'better' (I wish we had a way to measure what 'better' means...).

bendudson

Thanks @johnomotani !

Agreed a way to test these grids would be nice, but is a problem for another PR. This one is definitely an improvement.

johnomotani added 12 commits March 20, 2019 19:20

Hypnotoad: pass 'simple' setting to create_grid re-tries

98089fd

When create_grid fails, sometimes it changes some options and tries calling create_grid again. Previously the 'simple' setting was not being passed through to the new calls.

Hypnotoad: rename gen_surface->gen_surface_hypnotoad to avoid name clash

1a1eb26

Previously Hypnotoad would use the version from idllib if that was in the $IDL_PATH.

Hypnotoad: 'H' is derivative of integral

e642121

The quantity 'H' is the y-derivative of qinty, but qinty is an integral in y. Instead of calculating numerically the integral and derivative, just set 'H' equal to the integrand. This removes errors at the y-boundaries.

Hypnotoad: use I instead of sinty to calculate curvature

e801fbe

When calculating 'bxcvz' the integrated shear 'sinty' was used. This commit updates to use 'I' instead, so that when 'I' is set to zero for orthogonal coordinates the curvature is calculated consistently.

Hypnotoad: fix 'nnpol' update

d1902c5

Bump Hypnotoad version number

c77eb7c

1.1.4

Bump Hypnotoad version number

30bb0cb

1.2.0

Hypnotoad: don't smooth 'beta'

3049f4a

Makes results more consistent with/without y-boundary guard cells.

johnomotani mentioned this pull request Mar 20, 2019

Fix y-boundaries of Coordinates fields #1646

Merged

johnomotani added feature bugfix labels Mar 20, 2019

johnomotani mentioned this pull request Mar 20, 2019

Fix Coordinates initialization at boundaries #1607

Closed

johnomotani added the auxiliary-tools For items relating to the tools that come along with BOUT++ rather than the core library label Mar 21, 2019

bendudson reviewed Mar 25, 2019

View reviewed changes

Undo change to plot label 0 -> 0.D

53d18c8

Don't need double precision specification in a plot label.

bendudson reviewed Mar 25, 2019

View reviewed changes

bendudson previously approved these changes Mar 25, 2019

View reviewed changes

johnomotani added 2 commits March 25, 2019 18:44

Hypnotoad: refine position of starting line in grid_region_nonorth

ca936c6

Previously was done before calling grid_region_nonorth in the closed-field line only case (though this was recently commented out), but not in other cases. Now done for all cases.

Hypnotoad: remove some commented-out code

aba34cf

johnomotani dismissed bendudson’s stale review via aba34cf March 25, 2019 20:18

bendudson approved these changes Mar 28, 2019

View reviewed changes

bendudson merged commit 3b131a8 into next Mar 28, 2019

bendudson deleted the hypnotoad-y-boundary-guards branch March 28, 2019 21:15

johnomotani mentioned this pull request May 8, 2019

Hypnotoad non-orthogonal 'H' #1714

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Hypnotoad: double precision and y-boundary guard cells#1647

Hypnotoad: double precision and y-boundary guard cells#1647
bendudson merged 16 commits intonextfrom
hypnotoad-y-boundary-guards

johnomotani commented Mar 20, 2019

Uh oh!

bendudson commented Mar 24, 2019

Uh oh!

johnomotani commented Mar 25, 2019 •

edited

Loading

Uh oh!

bendudson Mar 25, 2019

Uh oh!

johnomotani commented Mar 25, 2019

Uh oh!

bendudson Mar 25, 2019

Uh oh!

johnomotani Mar 25, 2019

Uh oh!

bendudson Mar 25, 2019

Uh oh!

bendudson Mar 25, 2019

Uh oh!

bendudson left a comment

Uh oh!

johnomotani commented Mar 25, 2019

Uh oh!

johnomotani commented Mar 25, 2019

Uh oh!

bendudson left a comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

		@@ -1450,20 +1516,20 @@ FUNCTION create_nonorthogonal, F, R, Z, in_settings, critical=critical, $
		;; REPEAT BEGIN

Conversation

johnomotani commented Mar 20, 2019

Uh oh!

bendudson commented Mar 24, 2019

Uh oh!

johnomotani commented Mar 25, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

bendudson Mar 25, 2019

Choose a reason for hiding this comment

Uh oh!

johnomotani commented Mar 25, 2019

Uh oh!

bendudson Mar 25, 2019

Choose a reason for hiding this comment

Uh oh!

johnomotani Mar 25, 2019

Choose a reason for hiding this comment

Uh oh!

bendudson Mar 25, 2019

Choose a reason for hiding this comment

Uh oh!

bendudson Mar 25, 2019

Choose a reason for hiding this comment

Uh oh!

bendudson left a comment

Choose a reason for hiding this comment

Uh oh!

johnomotani commented Mar 25, 2019

Uh oh!

johnomotani commented Mar 25, 2019

Uh oh!

bendudson left a comment

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

johnomotani commented Mar 25, 2019 •

edited

Loading