-
Notifications
You must be signed in to change notification settings - Fork 90
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
PeriodicX option in cyclic Laplace solver #2557
Conversation
Upside-down "w" in executable name
If `periodicX = true` then the (kx,kz) = (0,0) mode is unconstrained and the Laplacian matrix inversion of vorticity to get potential becomes singular. The solution is to force solutions to have a zero average potential. This fix retains the tridiagonal form of the matrix by treating the kz=0 mode specially: - In the tridiagonal solve fix one location (x=xstart) to zero - After the solve, calculate and subtract the X average of the kz=0 modes i.e. set the (0,0) mode to zero. For the Hasegawa-wakatani model this results in solutions which are periodic in X, have zonal (kz=0) modes, and run at least as quickly as the non-periodic simulations.
Default periodic in X for easier analysis avoiding boundary effects
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Implementation looks good to me (minor comment below).
I think we should add a test. It's unfortunate that tests/integrated/test-laplace
tests against expected data saved in a binary file, so not ideal to add new cases to that test. Maybe we could test in tests/MMS/laplace
?
Also would be good to add a section to the manual mentioning periodicX
and how to use it. E.g. that you should only use laplace:type = cyclic
.
Should enable all FFT-based solvers using these tridiagonal matrix coefficients to produce sensible results for doubly-periodic (XZ) domains. Adds a zperiodic flag, because DST-based solvers don't need the special treatment.
@ZedThree I think this PR still needs a test, and |
Hi there, I just tried to compile this branch on Marconi, but it fails (strangely the master works). I then tried to merge on my local repo the branch into the master but I get several conflicts (I guess due to "next" being very different from "master"). Could you please look into it? Do you have intention to merge this anytime soon? |
So I'm having to guess, but if your error was
then the problem is the ancient version of CMake that's available on Marconi. A fix/workaround has been added to |
Hi @johnomotani, thanks a lot for your help. The quick fix suggested works (despite my poor instructions, sorry!). However, testing now the hasegawa-wakatani example and getting the following error after "gmake" Any hint? |
I suspect you may need to try using cmake to build the example if you have configured BOUT++ with cmake as well. @ZedThree may be able to help? |
Put into input_grids section for now.
Adapted from the Petsc3dAmg unit test. Doesn't yet include any actual tests.
No actual tests, but fixes compilation errors
Fails if no tests are defined, so add an empty test
Not available for 3D metrics
On MacOS it seems necessary to close and re-open the file to see the updated contents.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
clang-tidy made some suggestions
@@ -246,6 +247,23 @@ FieldPerp LaplaceCyclic::solve(const FieldPerp& rhs, const FieldPerp& x0) { | |||
cr->setCoefs(a, b, c); | |||
cr->solve(bcmplx, xcmplx); | |||
|
|||
if (localmesh->periodicX) { | |||
// Subtract X average of kz=0 mode | |||
BoutReal local[2] = { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: do not declare C-style arrays, use std::array<> instead [cppcoreguidelines-avoid-c-arrays]
BoutReal local[2] = {
^
for (int ix = xs; ix <= xe; ix++) { | ||
local[0] += xcmplx(0, ix - xs).real(); | ||
} | ||
BoutReal global[2]; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: do not declare C-style arrays, use std::array<> instead [cppcoreguidelines-avoid-c-arrays]
BoutReal global[2];
^
|
||
if (localmesh->periodicX) { | ||
// Subtract X average of kz=0 mode | ||
BoutReal local[ny + 1]; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: do not declare C VLA arrays, use std::vector<> instead [cppcoreguidelines-avoid-c-arrays]
BoutReal local[ny + 1];
^
local[ny] = static_cast<BoutReal>(xe - xs + 1); | ||
|
||
// Global reduce | ||
BoutReal global[ny + 1]; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: do not declare C VLA arrays, use std::vector<> instead [cppcoreguidelines-avoid-c-arrays]
BoutReal global[ny + 1];
^
Order of class members, some unused variables.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
clang-tidy made some suggestions
const Field2D *a, const Field2D *c1coef, const Field2D *c2coef, | ||
const Field2D *d, | ||
bool includeguards) { | ||
void Laplacian::tridagMatrix(dcomplex* avec, dcomplex* bvec, dcomplex* cvec, dcomplex* bk, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: function 'tridagMatrix' has cognitive complexity of 228 (threshold 25) [readability-function-cognitive-complexity]
void Laplacian::tridagMatrix(dcomplex* avec, dcomplex* bvec, dcomplex* cvec, dcomplex* bk,
^
src/invert/laplace/invert_laplace.cxx:434: +1, including nesting penalty of 0, nesting level increased to 1
ASSERT1(a->getLocation() == location);
^
include/bout/assert.hxx:38: expanded from macro 'ASSERT1'
if(!(condition)) { \
^
src/invert/laplace/invert_laplace.cxx:435: +1, including nesting penalty of 0, nesting level increased to 1
ASSERT1(c1coef->getLocation() == location);
^
include/bout/assert.hxx:38: expanded from macro 'ASSERT1'
if(!(condition)) { \
^
src/invert/laplace/invert_laplace.cxx:436: +1, including nesting penalty of 0, nesting level increased to 1
ASSERT1(c2coef->getLocation() == location);
^
include/bout/assert.hxx:38: expanded from macro 'ASSERT1'
if(!(condition)) { \
^
src/invert/laplace/invert_laplace.cxx:437: +1, including nesting penalty of 0, nesting level increased to 1
ASSERT1(d->getLocation() == location);
^
include/bout/assert.hxx:38: expanded from macro 'ASSERT1'
if(!(condition)) { \
^
src/invert/laplace/invert_laplace.cxx:448: +1, including nesting penalty of 0, nesting level increased to 1
if(!includeguards) {
^
src/invert/laplace/invert_laplace.cxx:449: +2, including nesting penalty of 1, nesting level increased to 2
if(!localmesh->firstX() || localmesh->periodicX)
^
src/invert/laplace/invert_laplace.cxx:449: +1
if(!localmesh->firstX() || localmesh->periodicX)
^
src/invert/laplace/invert_laplace.cxx:451: +2, including nesting penalty of 1, nesting level increased to 2
if(!localmesh->lastX() || localmesh->periodicX)
^
src/invert/laplace/invert_laplace.cxx:451: +1
if(!localmesh->lastX() || localmesh->periodicX)
^
src/invert/laplace/invert_laplace.cxx:462: +1, including nesting penalty of 0, nesting level increased to 1
if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) {
^
src/invert/laplace/invert_laplace.cxx:462: +1
if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) {
^
src/invert/laplace/invert_laplace.cxx:465: +1, including nesting penalty of 0, nesting level increased to 1
if(inner_boundary_flags & INVERT_BNDRY_ONE)
^
src/invert/laplace/invert_laplace.cxx:467: +1, including nesting penalty of 0, nesting level increased to 1
if(outer_boundary_flags & INVERT_BNDRY_ONE)
^
src/invert/laplace/invert_laplace.cxx:472: +1, including nesting penalty of 0, nesting level increased to 1
for(int ix=0;ix<=ncx;ix++) {
^
src/invert/laplace/invert_laplace.cxx:475: +2, including nesting penalty of 1, nesting level increased to 2
if (a != nullptr)
^
src/invert/laplace/invert_laplace.cxx:481: +1, including nesting penalty of 0, nesting level increased to 1
if(!localmesh->periodicX) {
^
src/invert/laplace/invert_laplace.cxx:482: +2, including nesting penalty of 1, nesting level increased to 2
if(localmesh->firstX()) {
^
src/invert/laplace/invert_laplace.cxx:487: +3, including nesting penalty of 2, nesting level increased to 3
if(!(inner_boundary_flags & (INVERT_RHS | INVERT_SET))) {
^
src/invert/laplace/invert_laplace.cxx:488: +4, including nesting penalty of 3, nesting level increased to 4
for(int ix=0;ix<inbndry;ix++)
^
src/invert/laplace/invert_laplace.cxx:493: +3, including nesting penalty of 2, nesting level increased to 3
if(kz == 0) {
^
src/invert/laplace/invert_laplace.cxx:495: +4, including nesting penalty of 3, nesting level increased to 4
if(inner_boundary_flags & INVERT_DC_GRAD && (inner_boundary_flags & INVERT_SET || inner_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:495: +1
if(inner_boundary_flags & INVERT_DC_GRAD && (inner_boundary_flags & INVERT_SET || inner_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:495: +1
if(inner_boundary_flags & INVERT_DC_GRAD && (inner_boundary_flags & INVERT_SET || inner_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:497: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:503: +1, nesting level increased to 4
else if(inner_boundary_flags & INVERT_DC_GRAD) {
^
src/invert/laplace/invert_laplace.cxx:505: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:511: +1, nesting level increased to 4
else if(inner_boundary_flags & INVERT_DC_GRADPAR) {
^
src/invert/laplace/invert_laplace.cxx:512: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++) {
^
src/invert/laplace/invert_laplace.cxx:518: +1, nesting level increased to 4
else if(inner_boundary_flags & INVERT_DC_GRADPARINV) {
^
src/invert/laplace/invert_laplace.cxx:519: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++) {
^
src/invert/laplace/invert_laplace.cxx:525: +1, nesting level increased to 4
else if (inner_boundary_flags & INVERT_DC_LAP) {
^
src/invert/laplace/invert_laplace.cxx:528: +5, including nesting penalty of 4, nesting level increased to 5
if (a != nullptr) {
^
src/invert/laplace/invert_laplace.cxx:530: +6, including nesting penalty of 5, nesting level increased to 6
if(ksq < 0.0)
^
src/invert/laplace/invert_laplace.cxx:534: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:540: +1, nesting level increased to 4
else if (inner_boundary_flags & INVERT_IN_CYLINDER){
^
src/invert/laplace/invert_laplace.cxx:576: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:583: +1, nesting level increased to 4
else {
^
src/invert/laplace/invert_laplace.cxx:585: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:593: +1, nesting level increased to 3
else {
^
src/invert/laplace/invert_laplace.cxx:595: +4, including nesting penalty of 3, nesting level increased to 4
if(inner_boundary_flags & INVERT_AC_GRAD && (inner_boundary_flags & INVERT_SET || inner_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:595: +1
if(inner_boundary_flags & INVERT_AC_GRAD && (inner_boundary_flags & INVERT_SET || inner_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:595: +1
if(inner_boundary_flags & INVERT_AC_GRAD && (inner_boundary_flags & INVERT_SET || inner_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:597: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:603: +1, nesting level increased to 4
else if(inner_boundary_flags & INVERT_AC_GRAD) {
^
src/invert/laplace/invert_laplace.cxx:605: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:611: +1, nesting level increased to 4
else if(inner_boundary_flags & INVERT_AC_LAP) {
^
src/invert/laplace/invert_laplace.cxx:613: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++) {
^
src/invert/laplace/invert_laplace.cxx:619: +1, nesting level increased to 4
else if (inner_boundary_flags & INVERT_IN_CYLINDER) {
^
src/invert/laplace/invert_laplace.cxx:622: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:626: +6, including nesting penalty of 5, nesting level increased to 6
if ((kz % 2) == 0){
^
src/invert/laplace/invert_laplace.cxx:629: +1, nesting level increased to 6
else {
^
src/invert/laplace/invert_laplace.cxx:634: +1, nesting level increased to 4
else {
^
src/invert/laplace/invert_laplace.cxx:637: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:645: +2, including nesting penalty of 1, nesting level increased to 2
if(localmesh->lastX()) {
^
src/invert/laplace/invert_laplace.cxx:650: +3, including nesting penalty of 2, nesting level increased to 3
if(!(outer_boundary_flags & (INVERT_RHS | INVERT_SET))) {
^
src/invert/laplace/invert_laplace.cxx:651: +4, including nesting penalty of 3, nesting level increased to 4
for (int ix=0;ix<outbndry;ix++) {
^
src/invert/laplace/invert_laplace.cxx:657: +3, including nesting penalty of 2, nesting level increased to 3
if(kz==0) {
^
src/invert/laplace/invert_laplace.cxx:659: +4, including nesting penalty of 3, nesting level increased to 4
if(outer_boundary_flags & INVERT_DC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:659: +1
if(outer_boundary_flags & INVERT_DC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:659: +1
if(outer_boundary_flags & INVERT_DC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:661: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<outbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:667: +1, nesting level increased to 4
else if(outer_boundary_flags & INVERT_DC_GRAD) {
^
src/invert/laplace/invert_laplace.cxx:669: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<outbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:675: +1, nesting level increased to 4
else if(inner_boundary_flags & INVERT_DC_GRADPAR) {
^
src/invert/laplace/invert_laplace.cxx:676: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++) {
^
src/invert/laplace/invert_laplace.cxx:682: +1, nesting level increased to 4
else if(inner_boundary_flags & INVERT_DC_GRADPARINV) {
^
src/invert/laplace/invert_laplace.cxx:683: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++) {
^
src/invert/laplace/invert_laplace.cxx:689: +1, nesting level increased to 4
else if (inner_boundary_flags & INVERT_DC_LAP) {
^
src/invert/laplace/invert_laplace.cxx:692: +5, including nesting penalty of 4, nesting level increased to 5
if (a != nullptr) {
^
src/invert/laplace/invert_laplace.cxx:694: +6, including nesting penalty of 5, nesting level increased to 6
if(ksq < 0.0)
^
src/invert/laplace/invert_laplace.cxx:698: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<inbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:704: +1, nesting level increased to 4
else {
^
src/invert/laplace/invert_laplace.cxx:707: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<outbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:715: +1, nesting level increased to 3
else {
^
src/invert/laplace/invert_laplace.cxx:717: +4, including nesting penalty of 3, nesting level increased to 4
if(outer_boundary_flags & INVERT_AC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:717: +1
if(outer_boundary_flags & INVERT_AC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:717: +1
if(outer_boundary_flags & INVERT_AC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) {
^
src/invert/laplace/invert_laplace.cxx:719: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<outbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:725: +1, nesting level increased to 4
else if(outer_boundary_flags & INVERT_AC_GRAD) {
^
src/invert/laplace/invert_laplace.cxx:727: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<outbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:733: +1, nesting level increased to 4
else if(outer_boundary_flags & INVERT_AC_LAP) {
^
src/invert/laplace/invert_laplace.cxx:735: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<outbndry;ix++) {
^
src/invert/laplace/invert_laplace.cxx:741: +1, nesting level increased to 4
else {
^
src/invert/laplace/invert_laplace.cxx:744: +5, including nesting penalty of 4, nesting level increased to 5
for (int ix=0;ix<outbndry;ix++){
^
src/invert/laplace/invert_laplace.cxx:752: +1, nesting level increased to 1
} else if (zperiodic and (kz == 0) and localmesh->firstX()) {
^
src/invert/laplace/invert_laplace.cxx:752: +1
} else if (zperiodic and (kz == 0) and localmesh->firstX()) {
^
Checking if this is the source of test suite failure with segfault.
Not the cause of Github test failures
Not set for all tests, resulting in warning about OMP_NUM_THREADS environment variable.
Attempt to identify origin of segfault on Github Actions
Possible linker confusion, with two ForwardOperator classes defined in separate files.
The segfault seems to be resolved.
If
periodicX = true
then the (kx,kz) = (0,0) mode is unconstrained and the Laplacian matrix inversion of vorticity to get potential becomes singular.The solution is to force solutions to have a zero average potential.
This fix retains the tridiagonal form of the matrix by treating the kz=0 mode specially:
modes i.e. set the (0,0) mode to zero.
For the Hasegawa-wakatani model this results in solutions which are periodic in X, have zonal (kz=0) modes, and run at least as quickly as the non-periodic simulations.
Fixes issue #2548