Skip to content

Commit

Permalink
Merge pull request #2649 from boutproject/par-bc-improve
Browse files Browse the repository at this point in the history
Improve ParallelBoundary code
  • Loading branch information
bendudson committed May 16, 2024
2 parents 85438dc + d99183e commit 22f0a03
Show file tree
Hide file tree
Showing 33 changed files with 463 additions and 280 deletions.
10 changes: 5 additions & 5 deletions examples/fci-wave-logn/boundary/BOUT.inp
Original file line number Diff line number Diff line change
Expand Up @@ -20,23 +20,23 @@ expand_divergence = false
background = 1e-06 # Background density

[all]
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2
bndry_all = neumann

[n]

zl = z / (2*pi)
function = fciwave:background + 1e-3*exp(-((x-0.7)/0.1)^2 - ((zl-0.3)/0.1)^2)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[logn]

function = log(n:function)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[v]

Expand Down
10 changes: 5 additions & 5 deletions examples/fci-wave-logn/div-integrate/BOUT.inp
Original file line number Diff line number Diff line change
Expand Up @@ -20,23 +20,23 @@ expand_divergence = false
background = 1e-06 # Background density

[all]
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2
bndry_all = neumann

[n]

zl = z / (2*pi)
function = fciwave:background + 1e-3*exp(-((x-0.7)/0.1)^2 - ((zl-0.3)/0.1)^2)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[logn]

function = log(n:function)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[v]

Expand Down
10 changes: 5 additions & 5 deletions examples/fci-wave-logn/expanded/BOUT.inp
Original file line number Diff line number Diff line change
Expand Up @@ -20,23 +20,23 @@ expand_divergence = true
background = 1e-06 # Background density

[all]
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2
bndry_all = neumann

[n]

zl = z / (2*pi)
function = fciwave:background + 1e-3*exp(-((x-0.7)/0.1)^2 - ((zl-0.3)/0.1)^2)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[logn]

function = log(n:function)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[v]

Expand Down
2 changes: 1 addition & 1 deletion examples/fci-wave-logn/fci-wave.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ class FCIwave : public PhysicsModel {

// Neumann boundaries simplifies parallel derivatives
Bxyz.applyBoundary("neumann");
Bxyz.applyParallelBoundary("parallel_neumann");
Bxyz.applyParallelBoundary("parallel_neumann_o2");
SAVE_ONCE(Bxyz);

Options::getRoot()->getSection("fciwave")->get("expand_divergence", expand_divergence,
Expand Down
10 changes: 5 additions & 5 deletions examples/fci-wave/div-integrate/BOUT.inp
Original file line number Diff line number Diff line change
Expand Up @@ -21,23 +21,23 @@ log_density = false # Evolve log(n)?
background = 1e-06 # Background density

[all]
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2
bndry_all = neumann

[n]

zl = z / (2*pi)
function = fciwave:background + 1e-3*exp(-((x-0.7)/0.1)^2 - ((zl-0.3)/0.1)^2)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[logn]

function = log(n:function)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[v]

Expand Down
10 changes: 5 additions & 5 deletions examples/fci-wave/div/BOUT.inp
Original file line number Diff line number Diff line change
Expand Up @@ -21,23 +21,23 @@ log_density = false # Evolve log(n)?
background = 1e-06 # Background density

[all]
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2
bndry_all = neumann

[n]

zl = z / (2*pi)
function = fciwave:background + 1e-3*exp(-((x-0.7)/0.1)^2 - ((zl-0.3)/0.1)^2)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[logn]

function = log(n:function)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[v]

Expand Down
2 changes: 1 addition & 1 deletion examples/fci-wave/fci-wave.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ class FCIwave : public PhysicsModel {

// Neumann boundaries simplifies parallel derivatives
Bxyz.applyBoundary("neumann");
Bxyz.applyParallelBoundary("parallel_neumann");
Bxyz.applyParallelBoundary("parallel_neumann_o2");
SAVE_ONCE(Bxyz);

SOLVE_FOR(nv);
Expand Down
10 changes: 5 additions & 5 deletions examples/fci-wave/logn/BOUT.inp
Original file line number Diff line number Diff line change
Expand Up @@ -21,23 +21,23 @@ log_density = true # Evolve log(n)?
background = 1e-06 # Background density

[all]
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2
bndry_all = neumann

[n]

zl = z / (2*pi)
function = fciwave:background + 1e-3*exp(-((x-0.7)/0.1)^2 - ((zl-0.3)/0.1)^2)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[logn]

function = log(n:function)

bndry_par_yup = parallel_neumann
bndry_par_ydown = parallel_neumann
bndry_par_yup = parallel_neumann_o2
bndry_par_ydown = parallel_neumann_o2

[nv]

Expand Down
6 changes: 3 additions & 3 deletions examples/laplace-petsc3d/data/BOUT.inp
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ mz = 128
function = mixmode(x, 1.)*mixmode(y, 2.)*mixmode(z, 3.)
bndry_xin = none
bndry_xout = none
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2

[rhs]
function = mixmode(x, 4.)*mixmode(y, 5.)*mixmode(z, 6.)
Expand All @@ -22,7 +22,7 @@ function = 1. + .1*mixmode(x, 10.)*mixmode(y, 11.)*mixmode(z, 12.)
[C2]
#function = 0.
function = .1*mixmode(x, 13.)*mixmode(y, 14.)*mixmode(z, 15.)
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2

[A]
function = 0.0
Expand All @@ -46,7 +46,7 @@ transform_from_field_aligned = false
[initial]
bndry_xin = neumann
bndry_xout = neumann
bndry_par_all = parallel_neumann
bndry_par_all = parallel_neumann_o2

[input1]
function = mixmode(x, 1.)*mixmode(z, 2.)
Expand Down
2 changes: 1 addition & 1 deletion include/bout/field_data.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ class BoundaryOpPar;
class Coordinates;
class Mesh;

class BoundaryRegion;
#include "bout/boundary_region.hxx"
class BoundaryRegionPar;
enum class BndryLoc;

Expand Down
6 changes: 3 additions & 3 deletions include/bout/mesh.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -55,16 +55,16 @@ class Mesh;
#include "bout/field_data.hxx"
#include "bout/options.hxx"

#include "fieldgroup.hxx"
#include "bout/fieldgroup.hxx"

class BoundaryRegion;
class BoundaryRegionPar;

#include "sys/range.hxx" // RangeIterator
#include "bout/sys/range.hxx" // RangeIterator

#include <bout/griddata.hxx>

#include "coordinates.hxx" // Coordinates class
#include "bout/coordinates.hxx" // Coordinates class

#include "bout/unused.hxx"

Expand Down
67 changes: 45 additions & 22 deletions include/bout/parallel_boundary_op.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ protected:
BoutReal getValue(const BoundaryRegionPar& bndry, BoutReal t);
};

template <class T>
template <class T, bool isNeumann = false>
class BoundaryOpParTemp : public BoundaryOpPar {
public:
using BoundaryOpPar::BoundaryOpPar;
Expand Down Expand Up @@ -89,51 +89,74 @@ public:
throw BoutException("Can't apply parallel boundary conditions to Field2D!");
}
void apply(Field3D& f) override { return apply(f, 0); }

void apply(Field3D& f, BoutReal t) override {
f.ynext(bndry->dir).allocate(); // Ensure unique before modifying

auto dy = f.getCoordinates()->dy;

for (bndry->first(); !bndry->isDone(); bndry->next()) {
BoutReal value = getValue(*bndry, t);
if (isNeumann) {
value *= dy[bndry->ind()];
}
static_cast<T*>(this)->apply_stencil(f, bndry, value);
}
}
};

//////////////////////////////////////////////////
// Implementations

class BoundaryOpPar_dirichlet : public BoundaryOpParTemp<BoundaryOpPar_dirichlet> {
class BoundaryOpPar_dirichlet_o1 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o1> {
public:
using BoundaryOpParTemp::BoundaryOpParTemp;

using BoundaryOpParTemp::apply;
void apply(Field3D& f, BoutReal t) override;
static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
bndry->dirichlet_o1(f, value);
}
};

class BoundaryOpPar_dirichlet_O3 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_O3> {
class BoundaryOpPar_dirichlet_o2 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o2> {
public:
using BoundaryOpParTemp::BoundaryOpParTemp;

using BoundaryOpParTemp::apply;
void apply(Field3D& f, BoutReal t) override;
static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
bndry->dirichlet_o2(f, value);
}
};

class BoundaryOpPar_dirichlet_interp
: public BoundaryOpParTemp<BoundaryOpPar_dirichlet_interp> {
class BoundaryOpPar_dirichlet_o3 : public BoundaryOpParTemp<BoundaryOpPar_dirichlet_o3> {
public:
using BoundaryOpParTemp::BoundaryOpParTemp;

using BoundaryOpParTemp::apply;
void apply(Field3D& f, BoutReal t) override;
static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
bndry->dirichlet_o3(f, value);
}
};

class BoundaryOpPar_neumann : public BoundaryOpParTemp<BoundaryOpPar_neumann> {
class BoundaryOpPar_neumann_o1
: public BoundaryOpParTemp<BoundaryOpPar_neumann_o1, true> {
public:
using BoundaryOpParTemp::BoundaryOpParTemp;

using BoundaryOpParTemp::apply;
void apply(Field3D& f, BoutReal t) override;
static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
bndry->neumann_o1(f, value);
}
};

class BoundaryOpPar_neumann_c2_simple
: public BoundaryOpParTemp<BoundaryOpPar_neumann_c2_simple> {
class BoundaryOpPar_neumann_o2
: public BoundaryOpParTemp<BoundaryOpPar_neumann_o2, true> {
public:
using BoundaryOpParTemp::BoundaryOpParTemp;
static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
bndry->neumann_o2(f, value);
}
};

using BoundaryOpParTemp::apply;
void apply(Field3D& f, BoutReal t) override;
class BoundaryOpPar_neumann_o3
: public BoundaryOpParTemp<BoundaryOpPar_neumann_o3, true> {
public:
using BoundaryOpParTemp::BoundaryOpParTemp;
static void apply_stencil(Field3D& f, const BoundaryRegionPar* bndry, BoutReal value) {
bndry->neumann_o3(f, value);
}
};

#endif // BOUT_PAR_BNDRY_OP_H

0 comments on commit 22f0a03

Please sign in to comment.