[WIP] Implementation of absorbing boundary layer

[cjvogl]

Implemented absorbing boundary layer and modified acoustics_2d_radial example to add an ABL. Makefile_abl, setrun_abl.py, and setplot_abl.py are provided to make use of the new ABL routines.

This PR depends on https://github.com/clawpack/clawutil/pull/129 and should not be merged until that PR is merged.

[rjleveque]

Thanks @cjvogl for working on this. I'll test it out a bit more and it would be great if we can get this incorporated soon.

I'm not sure why the travis tests failed, it seems to have been a server error on the apt-get update, so I'll try restarting the tests. They pass on my laptop.

[rjleveque]

OK, the tests pass!

[rjleveque]

@cjvogl: looking at this some more, it is not clear to me why you need num_dim extra aux arrays. Isn't there just a single scale factor in each cell, at least along the edges? Is it only because of corners where two absorbing layers overlap where you potentially need different scale factors in each direction?

It would be nice to avoid extra aux arrays as much as possible, particularly for the 3D case and with AMR. Especially since these these arrays just have the value 1 everywhere except in the layers I guess?

I also notice that in flux2f.90 you multiply s, amdp, apdq by aux2(num_aux-2+ixy,:) everywhere even though this value is mostly 1. Can we avoid this somehow?

[cjvogl]

@rjleveque: once in the mapped coordinate frame, the system to be solved is Q_t + s(x)AQ_x + s(y)BQ_y + s(z)*CQ_z = 0. Thus, we either need to have the values of s(x_i), s(y_j), s(z_k) or maybe x_i, y_j, and z_k available to wherever we want to implement the layer (currently in flux2.f90). You are correct that for the majority of cells in the layer, all but one of these values is 1, and that it is only in the corners where all 3 will have values less than 1.

As for the extra memory footprint, I agree that this approach sacrifices memory to minimize how invasive the ABL implementation is. @mandli and I briefly discussed some alternative approaches, where it seemed like reserving memory for these values only in the absorbing layer being the most efficient. I think we also agreed that @mjberger would need to be brought on-board to discuss the best way to do this.

As for the multiply-by-one inefficiency, I can certain replace the related do loops with a where loop. I don't have enough experience to to speak to the efficiency of one approach (nested loops without logic) over the other (a "where" statement that involves logic) once the compiler optimizations happen though. If the "where loop" is the way to go, I'll get that changed.

[mjberger]

i don’t know what we’re talking about really, but can it be computed rather than stored? Math is cheap, memory bandwidth is slow.

On Aug 20, 2018, at 2:26 PM, Chris Vogl notifications@github.com wrote:

@rjleveque https://github.com/rjleveque: once in the mapped coordinate frame, the system to be solved is (in 2d) Q_t + s(x)AQ_x + s(y)BQ_y + s(z)*CQ_z = 0. Thus, we either need to have the values of s(x_i), s(y_j), s(z_k) or maybe x_i, y_j, and z_k available to wherever we want to implement the layer (currently in flux2.f90). You are correct that for the majority of cells in the layer, all but one of these values is 1, and that it is only in the corners where all 3 will have values less than 1.

As for the extra memory footprint, I agree that this approach sacrifices memory to minimize how invasive the ABL implementation is. @mandli https://github.com/mandli and I briefly discussed some alternative approaches, where it seemed like reserving memory for these values only in the absorbing layer being the most efficient. I think we also agreed that @mjberger https://github.com/mjberger would need to be brought on-board to discuss the best way to do this.

As for the multiply-by-one inefficiency, I can certain replace the related do loops with a where loop. I don't have enough experience to to speak to the efficiency of one approach (nest loops without logic) over the other (a "where" statement that involves logic) once the compiler optimizations happen though. If the "where loop" is the way to go, I'll get that changed.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/clawpack/classic/pull/83#issuecomment-414468840, or mute the thread https://github.com/notifications/unsubscribe-auth/AA1DC0oylLCLKP_7nmCFLs27UK2Ia0Qmks5uSymMgaJpZM4VdSCS.

[mandli]

@mjberger I can give you a brief synopsis of the issue in the near future but I will try summarize. The basic problem is that we want to save the speeds coming out of the Riemann solvers so that the BCs can be set to slow down these speeds appropriately. @cjvogl has been using aux arrays for this which of course even in classic is using more memory than needed but for AMR this becomes a very large overhead penalty given we only need these at the true boundaries of the domain. The problem then becomes detecting if we are at a true domain boundary and having appropriately sized storage at the edges or redoing the calculations at the boundaries when necessary and having all the correct information required.

[cjvogl]

@mjberger: yes, I apologize for the lack of details... a write-up of this is in progress. To expand on what @mandli said, we need to be able to compute an scaling factor that is applied to the wave info that comes out of the Riemann solver. In order to compute that scaling factor, we need to know where the grid cell edge is relative to the absorbing layer (that surrounds the domain). I dont recall there being an easy way to obtain grid cell location from flux2.f90... hopefully I’m mistaken.

On Aug 21, 2018, at 9:01 AM, Kyle Mandli notifications@github.com wrote:

@mjberger I can give you a brief synopsis of the issue in the near future but I will try summarize. The basic problem is that we want to save the speeds coming out of the Riemann solvers so that the BCs can be set to slow down these speeds appropriately. @cjvogl has been using aux arrays for this which of course even in classic is using more memory than needed but for AMR this becomes a very large overhead penalty given we only need these at the true boundaries of the domain. The problem then becomes detecting if we are at a true domain boundary and having appropriately sized storage at the edges or redoing the calculations at the boundaries when necessary and having all the correct information required.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[rjleveque]

It would be nice if we knew the x-y coordinates of cells in flux2 but I don't think we do. Note that in GeoClaw we have to store info based on the latitude in aux also for this reason. Maybe something to think about more generally at some point.

Since we need to get v5.5.0 out, I suggest we not merge this just yet and have some more joint discussion of the best way to do this before proceeding.

[cjvogl]

Not sure why the Travis checks failed, but I tried a refactor of the ABL implementation based on some discussions with @rjleveque.

I tried for a good balance between multiplying-by-one in a large section of the for loops that needed to be modified and creating too many additional data structures... let me know what you all think.

[mandli]

The following was in the Travis log:

Output from Test Advection2DAnnulusTest
Started 2019/03/06-03:53.11
================================================================================
Paths:  /tmp/tmpea2ppqn3  /home/travis/build/clawpack/classic/tests/advection_2d_annulustouch /home/travis/build/clawpack/classic/clawpack/classic/src/2d/abl_module.mod; gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/abl_module.f90 -J/home/travis/build/clawpack/classic/clawpack/classic/src/2d   -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/abl_module.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/driver.f90                       -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/driver.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/claw2ez.f                    -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/claw2ez.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/claw2.f                      -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/claw2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/bc2.f                    -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/bc2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/b4step2.f90                      -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/b4step2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/step2.f90                    -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/step2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/step2ds.f90                      -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/step2ds.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/dimsp2.f                     -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/dimsp2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/flux2.f90                    -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/flux2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/copyq2.f                     -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/copyq2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/inlinelimiter.f90                    -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/inlinelimiter.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/src2.f90                     -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/src2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/out2.f                       -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/out2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/restart2.f                       -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/restart2.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/classic/src/2d/opendatafile.f                       -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/opendatafile.o
gfortran -c -cpp mapc2p.f                     -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o mapc2p.o
gfortran -c -cpp qinit.f                      -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o qinit.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/riemann/src/rpn2_vc_advection.f90                       -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/riemann/src/rpn2_vc_advection.o
gfortran -c -cpp /home/travis/build/clawpack/classic/clawpack/riemann/src/rpt2_vc_advection.f90                       -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o /home/travis/build/clawpack/classic/clawpack/riemann/src/rpt2_vc_advection.o
gfortran -c -cpp setaux.f                     -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o setaux.o
gfortran -c -cpp setprob.f                    -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o setprob.o
gfortran -c -cpp stream.f                     -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/    -o stream.o
gfortran     /home/travis/build/clawpack/classic/clawpack/classic/src/2d/abl_module.o     /home/travis/build/clawpack/classic/clawpack/classic/src/2d/driver.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/claw2ez.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/claw2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/bc2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/b4step2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/step2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/step2ds.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/dimsp2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/flux2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/copyq2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/inlinelimiter.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/src2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/out2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/restart2.o /home/travis/build/clawpack/classic/clawpack/classic/src/2d/opendatafile.o mapc2p.o qinit.o /home/travis/build/clawpack/classic/clawpack/riemann/src/rpn2_vc_advection.o /home/travis/build/clawpack/classic/clawpack/riemann/src/rpt2_vc_advection.o setaux.o setprob.o stream.o   -I/home/travis/build/clawpack/classic/clawpack/classic/src/2d/ -L/home/travis/build/clawpack/classic/tests/advection_2d_annulus/   -o xclaw                         
Reading data file: claw.data   
         first 5 lines are comments and will be skipped
 +++ mx, my =           40         120
 running...

Reading data file: setprob.data             
         first 5 lines are comments and will be skipped
*** in opendatafile, file not found:abl.data
==> runclaw: Will take data from  /tmp/tmpea2ppqn3
==> runclaw: Will write output to  /tmp/tmpea2ppqn3
==> runclaw: Removing all old fort/gauge files in  /tmp/tmpea2ppqn3
==> Running with command:
    /tmp/tmpea2ppqn3/xclaw
==> runclaw: Finished executing
==> runclaw: Done executing /tmp/tmpea2ppqn3/xclaw via clawutil.runclaw.py
==> runclaw: Output is in  /tmp/tmpea2ppqn3
Test path and class info:
  class: <class 'tests.advection_2d_annulus.regression_tests.Advection2DAnnulusTest'>
  class file: /home/travis/build/clawpack/classic/tests/advection_2d_annulus/regression_tests.py
  test path: /home/travis/build/clawpack/classic/tests/advection_2d_annulus
  temp path: /tmp/tmpea2ppqn3
Errors from Test Advection2DAnnulusTest
Started 2019/03/06-03:53.11
================================================================================
Done. Your build exited with 1.

It looks like abl.data was not found?

[cjvogl]

Ah, that's right: abl.data will not be produced unless data.py from clawutil #129 is used. Thanks for grabbing the log info.