Add region arguments for toFieldAligned/fromFieldAligned by johnomotani · Pull Request #1432 · boutproject/BOUT-dev

johnomotani · 2018-12-10T16:45:40Z

Add region arguments (where it is valid to give RGN_NOX or RGN_NOBNDRY) so functions like toFieldAligned/fromFieldAligned can be used both when y-guard cells have been set and when they have not. e.g. when using fromFieldAligned on the results of interpolation or derivatives, y-guard cells cannot be included, but when using toFieldAligned on the inputs to interpolation or derivatives, y-guard cells must be included.

Pass RGN_NOBNDRY to fromFieldAligned in various places where the guard cells aren't valid.

Pulled out of #1176; is independent of other stuff in that PR, which now probably won't be merged.

Add region arguments (where it is valid to give RGN_NOX or RGN_NOBNDRY) so functions like toFieldAligned/fromFieldAligned can be used both when y-guard cells have been set and when they have not. e.g. when using fromFieldAligned on the results of interpolation or derivatives, y-guard cells cannot be included, but when using toFieldAligned on the inputs to interpolation or derivatives, y-guard cells must be included.

bendudson · 2018-12-12T07:16:19Z

-const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle) {
+const Field3D ShiftedMetric::shiftZ(const Field3D &f, const Field2D &zangle, const REGION region) {
  ASSERT1(&mesh == f.getMesh());
+  ASSERT1(region == RGN_NOX || region == RGN_NOBNDRY); // Never calculate x-guard cells here


Not sure this is necessary; it may be inefficient to perform shifting in the X guard cells, but shouldn't cause a problem. There are some cases where it may be useful e.g. calculating D2DXDY

I think the problem I was originally trying to avoid was when shifting in the corner guard cells. I remember getting 'uninitialized value' warnings from Valgrind that I wanted to protect against. They must have been coming from interpolated zShift fields though, so probably the right solution is to make sure all points in zShifts get initialized to something even when checking is turned off, instead of trying to force the corners to not be used here.

I'll get rid of this check 👍

Following discussion in PR #1432: it is not necessary to restrict the regions for to/fromFieldAligned or shiftZ to RGN_NOX and RGN_NOBNDRY as was being done previously.

d7919 · 2018-12-13T11:39:57Z

-      shiftZ(f(jx,jy), phs[jx][jy], result(jx,jy));
-    }
+
+  for (const auto &i : mesh.getRegion2D(REGION_STRING(region))) {


Could this use the BOUT_FOR(_SERIAL )macro rather than the range based loop?

d7919 · 2018-12-13T11:40:07Z

+  // (Note valgrind complains about corner guard cells if we try to loop over
+  // the whole grid, because zShift is not initialized in the corner guard
+  // cells.)
+  for(const auto &i : mesh.getRegion2D(REGION_STRING(region))) {


d7919 · 2018-12-14T12:57:35Z

+  BOUT_FOR(i, zShift.getRegion("RGN_ALL")) {
+    for(int jz=0;jz<nmodes;jz++) {
+      BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+      fromAlignedPhs[i.x()][i.y()][jz] = dcomplex(cos(kwave*zShift(i.x(),i.y())) , -sin(kwave*zShift(i.x(),i.y())));


I think you can just index zShift using i (i.e. zShift[i]) here.

d7919 · 2018-12-14T12:58:17Z

-      for(int jz=0;jz<nmodes;jz++) {
-  	BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+  BOUT_FOR(i, zShift.getRegion("RGN_ALL")) {
+    BoutReal yupShift = zShift(i.x(),i.y()) - zShift(i.x(),i.y()+1);


Index with i and i.yp() possible here I think.

d7919 · 2018-12-14T12:58:51Z

-      shiftZ(&(f(jx,jy+1,0)), yupPhs[jx][jy], &(yup(jx,jy+1,0)));
-    }
+  BOUT_FOR(i, f.getRegion("RGN_NOX")) {
+    shiftZ(&(f(i.x(),i.y()+1,0)), yupPhs[i.x()][i.y()], &(yup(i.x(),i.y()+1,0)));


Again I think you can use i.yp().

@d7919: This loop should be using getRegion2D, so can't use i.yp for f. I've updated to follow the rest of your comments.

@johnomotani it's possible to index a Field3D with an Ind2D, you just need to specify the z index as well. This uses a multiply and add to get the corresponding Ind3D as opposed to three divides and a modulus to call i.x() and i.y(), as such I'd expect the Ind2D approach to be more efficient. I'll add these changes to the branch that fixes the OpenMP stuff

d7919 · 2018-12-14T12:59:03Z

-      shiftZ(&(f(jx,jy-1,0)), ydownPhs[jx][jy], &(ydown(jx,jy-1,0)));
-    }
+  BOUT_FOR(i, f.getRegion("RGN_NOX")) {
+    shiftZ(&(f(i.x(),i.y()-1,0)), ydownPhs[i.x()][i.y()], &(ydown(i.x(),i.y()-1,0)));


d7919 · 2018-12-14T13:00:51Z

-    }
+
+  BOUT_FOR(i, f.getRegion(region)) {
+    shiftZ(f(i.x(),i.y()), phs[i.x()][i.y()], result(i.x(),i.y()));


Could you replace f(i.x(), i.y()) with &f(i.x(), i.y(), 0) and the same for result. Note I think you want f.getRegion2D(region) (or f.getMesh()->getRegion2D(region)) as currently you're looping over z as well.

d7919 · 2018-12-14T13:01:09Z

+  // the whole grid, because zShift is not initialized in the corner guard
+  // cells.)
+  BOUT_FOR(i, f.getRegion(region)) {
+    shiftZ(f(i.x(), i.y()), mesh.LocalNz, zangle(i.x(),i.y()), result(i.x(), i.y()));


See previous comments.

…omFieldAligned

Regions used in BOUT_FOR loops must exclude y-guard cells, like the previous nested loops.

johnomotani · 2018-12-16T16:35:50Z

test-yupdown is failing with OpenMP. The test passes if I change the loop in this shiftZ method from BOUT_FOR to BOUT_FOR_SERIAL

BOUT-dev/src/mesh/parallel/shiftedmetric.cxx

Lines 138 to 152 in baa630c

    
           const Field3D ShiftedMetric::shiftZ(const Field3D &f, const arr3Dvec &phs, const REGION region) { 
        
             ASSERT1(&mesh == f.getMesh()); 
        
             if(mesh.LocalNz == 1) 
        
               return f; // Shifting makes no difference 
        
             Field3D result(&mesh); 
        
             result.allocate(); 
        
             BOUT_FOR(i, mesh.getRegion2D(REGION_STRING(region))) { 
        
               shiftZ(&f(i.x(), i.y(), 0), phs[i.x()][i.y()], &result(i.x(), i.y(), 0)); 
        
             } 
        
             return result; 
        
           }

I can't see why the loop should not parallelize though. The loop body for each i should be unaffected by the result for any other i... Anyone got any hints for debugging?

d7919 · 2018-12-17T07:42:12Z

I don't have the code in front of me but done of these methods use a class member working space (vector) which would lead to race competition etc. I think this could probably be removed from class and made a routine local Array without performance penalty. Can look later if any help

Those should help improve OpenMP safety.

d7919 · 2018-12-17T12:16:51Z

I've pushed a branch at https://github.com/boutproject/BOUT-dev/tree/region-arguments-tofromFieldAligned-shiftZ-openmp which might help fix this.

… possible.

…dAligned-shiftZ-openmp Use function local storage in shiftZ rather than member storage

ZedThree · 2019-02-12T15:57:49Z

This needs updating now #1345 is in

…omFieldAligned

Need to set number of y-guard cells to 2 before creating default regions. Otherwise BOUT_FOR loops go out of bounds.

johnomotani · 2019-02-12T19:03:20Z

I've merged next now.

This PR does two things:

adds and makes use of REGION arguments for the toFieldAligned fromFieldAligned methods.
Uses Tensor for the cached phases in ShiftedMetric instead of arr3Dvec

I think 1 is good and should be merged. @ZedThree do you want 2, and if so should I split it into a separate PR?

RGN_ALL is no longer the default for toFieldAligned/fromFieldAligned. Force it to be used for the unit test.

johnomotani · 2019-02-12T22:48:29Z

This is failing with OpenMP enabled again. If I switch the BOUT_FOR to a BOUT_FOR_SERIAL in this loop, it passes

BOUT-dev/src/mesh/parallel/shiftedmetric.cxx

Lines 193 to 213 in e1857ec

    
           for (auto& phase : phases) { 
        
             // In C++17 std::vector::emplace_back returns a reference, which 
        
             // would be very useful here! 
        
             results.emplace_back(&mesh); 
        
             auto& current_result = results.back(); 
        
             current_result.allocate(); 
        
             current_result.setLocation(f.getLocation()); 
        
             BOUT_FOR(i, mesh.getRegion2D("RGN_NOY")) { 
        
               // Deep copy the FFT'd field 
        
               Array<dcomplex> shifted_temp(nmodes); 
        
               for (int jz = 1; jz < nmodes; ++jz) { 
        
                 shifted_temp[jz] = f_fft(i.x(), i.y() + phase.y_offset, jz) 
        
                                    * phase.phase_shift(i.x(), i.y(), jz); 
        
               } 
        
               irfft(shifted_temp.begin(), mesh.LocalNz, 
        
                     &current_result(i.yp(phase.y_offset), 0)); 
        
             } 
        
           }

Can anyone see the race condition or other error?

d7919 · 2019-02-13T08:40:32Z

Probably unrelated but

BOUT-dev/src/mesh/parallel/shiftedmetric.cxx

Line 205 in e1857ec

for (int jz = 1; jz < nmodes; ++jz) {

doesn't set the jz=0 point which I think we probably should as it's not initialised anywhere (this could just be shifted_temp[0] = 0 before the loop)

d7919 · 2019-02-13T08:44:57Z

The only place inside the loop that I could see an issue arising is with

BOUT-dev/src/mesh/parallel/shiftedmetric.cxx

Line 211 in e1857ec

&current_result(i.yp(phase.y_offset), 0));

As we're writing to current_result and this is a shared object. However, it looks like it should be fine as i.yp(phase.y_offset) should be unique for each i which is the OpenMP loop variable.

ZedThree · 2019-02-13T09:31:24Z

Probably unrelated but

BOUT-dev/src/mesh/parallel/shiftedmetric.cxx

Line 205 in e1857ec

for (int jz = 1; jz < nmodes; ++jz) {
doesn't set the jz=0 point which I think we probably should as it's not initialised anywhere (this could just be shifted_temp[0] = 0 before the loop)

L232-233 has been clobbered, it should look like:

BOUT-dev/src/mesh/parallel/shiftedmetric.cxx

Lines 232 to 233 in ace3f77

    
           Array<dcomplex> shifted_temp(f_fft[jx][jy + phase.y_offset]); 
        
           shifted_temp.ensureUnique();

Prevents errors with OpenMP.

johnomotani · 2019-02-13T10:32:15Z

@ZedThree I'd changed f_fft to a Tensor so couldn't get an Array from it to initialize shifted_temp with. I thought it would be OK to move getting values out of f_fft down into the loop instead, using shifted_temp[jz] = f_fft(i.x(), i.y() + phase.y_offset, jz) * phase.phase_shift(i.x(), i.y(), jz); instead of shifted_temp[jz] *= phase.phase_shift(i.x(), i.y(), jz);.

I've just pushed a commit replacing Tensor<dcomplex> f_fft with Matrix<Array<dcomplex>> f_fft so the structure looks a lot more like it did previously. The unit tests now pass with OpenMP. I still don't understand what the problem was though... I tried adding shifted_temp.ensureUnique() and shifted_temp[0] = 0., giving the method below, but still got the same errors

std::vector<Field3D>
ShiftedMetric::shiftZ(const Field3D& f,
                      const std::vector<ParallelSlicePhase>& phases) const {

  const int nmodes = mesh.LocalNz / 2 + 1;

  // FFT in Z of input field at each (x, y) point
  Tensor<dcomplex> f_fft(mesh.LocalNx, mesh.LocalNy, nmodes);

  BOUT_FOR(i, mesh.getRegion2D("RGN_ALL")) {
    rfft(&f(i, 0), mesh.LocalNz, &f_fft(i.x(), i.y(), 0));
  }

  std::vector<Field3D> results{};

  for (auto& phase : phases) {
    // In C++17 std::vector::emplace_back returns a reference, which
    // would be very useful here!
    results.emplace_back(&mesh);
    auto& current_result = results.back();
    current_result.allocate();
    current_result.setLocation(f.getLocation());
  
    BOUT_FOR(i, mesh.getRegion2D("RGN_NOY")) {
      // Deep copy the FFT'd field
      Array<dcomplex> shifted_temp(nmodes);
      shifted_temp.ensureUnique();
    
      shifted_temp[0] = 0.;
      for (int jz = 1; jz < nmodes; ++jz) {
        shifted_temp[jz] = f_fft(i.x(), i.y() + phase.y_offset, jz)
                           * phase.phase_shift(i.x(), i.y(), jz);
      }
  
      irfft(shifted_temp.begin(), mesh.LocalNz,
            &current_result(i.yp(phase.y_offset), 0));
    }
  } 
     
  return results;
}

Swap some array entries around so that 'input' in ShiftedMetricTest is not constant in x or y. Being constant could mask errors.

johnomotani · 2019-02-13T22:35:25Z

The previous test failure was because I'd messed up an index in shiftZ; I fixed it in c2f9eff, but was concerned that it hadn't been caught by the unit tests. It was missed because the input for the unit tests was constant in x and y, so I swapped some pairs of entries in the input field, and corresponding ones in the expected results. The updated unit test would catch the error.

bendudson · 2019-02-21T11:24:37Z

Is this PR good to go in now, if the conflict is fixed?

bendudson · 2019-02-21T11:45:31Z

  /// Transform a field into field-aligned coordinates
-  const Field3D toFieldAligned(const Field3D &f) {
-    return getParallelTransform().toFieldAligned(f);
+  const Field3D toFieldAligned(const Field3D &f, const REGION region = RGN_NOX) {


Should the default just be RGN_ALL since (I think) that's what it was before?

ZedThree · 2019-02-21T12:05:03Z

+  BOUT_FOR(i, mesh.getRegion2D("RGN_ALL")) {
+    for(int jz=0;jz<nmodes;jz++) {
+      BoutReal kwave=jz*2.0*PI/zlength; // wave number is 1/[rad]
+      fromAlignedPhs(i.x(), i.y(), jz) =


Is it worth caching the values of i.x(), i.y() outside the z-loop? Probably won't make much difference in the constructor, but could be more important in the other methods

(also clang-format please 🙂)

johnomotani mentioned this pull request Dec 10, 2018

Shiftedmetric updates #1176

Closed

bendudson reviewed Dec 12, 2018

View reviewed changes

johnomotani added 2 commits December 12, 2018 12:03

Remove region restrictions in ShiftedMetric::shiftZ

0f6a0fe

Following discussion in PR #1432: it is not necessary to restrict the regions for to/fromFieldAligned or shiftZ to RGN_NOX and RGN_NOBNDRY as was being done previously.

Add a few more explicit region arguments to fromFieldAligned calls

468f816

d7919 reviewed Dec 13, 2018

View reviewed changes

Replace nested loops in shiftedmetric.cxx with BOUT_FOR

5c6c0a0

d7919 reviewed Dec 14, 2018

View reviewed changes

johnomotani added 3 commits December 14, 2018 16:05

Fix ShiftedMetric BOUT_FOR loops

bd23c71

Merge remote-tracking branch 'origin/next' into region-arguments-tofr…

69c3fc6

…omFieldAligned

Fix regions for loops in ShiftedMetric

baa630c

Regions used in BOUT_FOR loops must exclude y-guard cells, like the previous nested loops.

johnomotani mentioned this pull request Dec 16, 2018

Coordinate system labelling #1459

Closed

Use function local storage in shiftZ rather than member storage

2c31bed

Those should help improve OpenMP safety.

d7919 and others added 4 commits December 17, 2018 12:46

Use Tensor for 3D storage in shifted metric rather than nested vectors

d49cb3a

Use Ind2D indexing in shifted metric rather than three integers where…

bd1f5e2

… possible.

Merge pull request #1462 from boutproject/region-arguments-tofromFiel…

5f6b407

…dAligned-shiftZ-openmp Use function local storage in shiftZ rather than member storage

Merge branch 'next' into region-arguments-tofromFieldAligned

b46942d

ZedThree mentioned this pull request Jan 29, 2019

Add multiple parallel slices #1345

Merged

johnomotani added 2 commits February 12, 2019 18:51

Merge remote-tracking branch 'origin/next' into region-arguments-tofr…

dee496b

…omFieldAligned

Fix ShiftedMetric unit test

1cb55e1

Need to set number of y-guard cells to 2 before creating default regions. Otherwise BOUT_FOR loops go out of bounds.

Fix ShiftedMetric unit test: use RGN_ALL explicitly

e1857ec

RGN_ALL is no longer the default for toFieldAligned/fromFieldAligned. Force it to be used for the unit test.

Use Matrix<Array<dcomplex>> instead of Tensor<dcomplex> for f_fft

52b6b7e

Prevents errors with OpenMP.

johnomotani added 2 commits February 13, 2019 21:16

Restore missing y_offset in ShiftedMetric::shiftZ

c2f9eff

Make input for ShiftedMetric unit tests not constant in y

a308542

Swap some array entries around so that 'input' in ShiftedMetricTest is not constant in x or y. Being constant could mask errors.

boutproject deleted a comment from codecov-io Feb 13, 2019

Merge branch 'next' into region-arguments-tofromFieldAligned

21667d2

bendudson reviewed Feb 21, 2019

View reviewed changes

ZedThree reviewed Feb 21, 2019

View reviewed changes

Comment thread tests/unit/mesh/parallel/test_shiftedmetric.cxx

ZedThree reviewed Feb 21, 2019

View reviewed changes

johnomotani added 4 commits February 21, 2019 12:43

Set the default REGION of toFieldAligned and fromFieldAligned to RGN_ALL

ebdbc77

Add comment on shuffling input for ShiftedMetric unit test

be9cd55

Cache i.x() and i.y() outside z-loops in ShiftedMetric

adf7c66

Clang formatting

33801ce

bendudson approved these changes Feb 21, 2019

View reviewed changes

ZedThree merged commit d113036 into next Feb 21, 2019

ZedThree deleted the region-arguments-tofromFieldAligned branch February 21, 2019 16:17

Uh oh!

Conversation

johnomotani commented Dec 10, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

johnomotani commented Dec 16, 2018

Uh oh!

d7919 commented Dec 17, 2018

Uh oh!

d7919 commented Dec 17, 2018

Uh oh!

ZedThree commented Feb 12, 2019

Uh oh!

johnomotani commented Feb 12, 2019

Uh oh!

johnomotani commented Feb 12, 2019

Uh oh!

d7919 commented Feb 13, 2019

Uh oh!

d7919 commented Feb 13, 2019

Uh oh!

ZedThree commented Feb 13, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

johnomotani commented Feb 13, 2019

Uh oh!

johnomotani commented Feb 13, 2019

Uh oh!

bendudson commented Feb 21, 2019

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

johnomotani commented Dec 10, 2018 •

edited

Loading

ZedThree commented Feb 13, 2019 •

edited

Loading