Skip to content

Commit cf71910

Browse files
malamastmalamast
committed
boundary_standard: Added a legacy implementation of 1st order Dirichlet and Neumann BCs.
- This is what we used to do in Bout++ version 3. - Higher order extrapolation can be unstable sometimes so I added back the 1st order BC. You can use thm in the input file by setting e.g. dirichlet_o1. - Implementation for staggered grids is not done.
1 parent df922ab commit cf71910

3 files changed

Lines changed: 321 additions & 2 deletions

File tree

‎include/bout/boundary_standard.hxx‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,30 @@ private:
3333
BoutReal val;
3434
};
3535

36+
/// Dirichlet (set to zero) boundary condition
37+
class BoundaryDirichlet_O1 : public BoundaryOp {
38+
public:
39+
BoundaryDirichlet_O1() : gen(nullptr) {}
40+
BoundaryDirichlet_O1(BoundaryRegion* region, std::shared_ptr<FieldGenerator> g)
41+
: BoundaryOp(region), gen(std::move(g)) {}
42+
43+
using BoundaryOp::clone;
44+
BoundaryOp* clone(BoundaryRegion* region, const std::list<std::string>& args) override;
45+
46+
using BoundaryOp::apply;
47+
void apply(Field2D& f) override;
48+
void apply(Field2D& f, BoutReal t) override;
49+
void apply(Field3D& f) override;
50+
void apply(Field3D& f, BoutReal t) override;
51+
52+
using BoundaryOp::apply_ddt;
53+
void apply_ddt(Field2D& f) override;
54+
void apply_ddt(Field3D& f) override;
55+
56+
private:
57+
std::shared_ptr<FieldGenerator> gen; // Generator
58+
};
59+
3660
/// Dirichlet (set to zero) boundary condition
3761
class BoundaryDirichlet : public BoundaryOp {
3862
public:
@@ -163,6 +187,30 @@ public:
163187
void apply(Field3D& f) override;
164188
};
165189

190+
/// Neumann (zero-gradient) boundary condition, using 1st order on boundary
191+
class BoundaryNeumann_O1 : public BoundaryOp {
192+
public:
193+
BoundaryNeumann_O1() : gen(nullptr) {}
194+
BoundaryNeumann_O1(BoundaryRegion* region, std::shared_ptr<FieldGenerator> g)
195+
: BoundaryOp(region), gen(std::move(g)) {}
196+
197+
using BoundaryOp::clone;
198+
BoundaryOp* clone(BoundaryRegion* region, const std::list<std::string>& args) override;
199+
200+
using BoundaryOp::apply;
201+
void apply(Field2D& f) override;
202+
void apply(Field2D& f, BoutReal t) override;
203+
void apply(Field3D& f) override;
204+
void apply(Field3D& f, BoutReal t) override;
205+
206+
using BoundaryOp::apply_ddt;
207+
void apply_ddt(Field2D& f) override;
208+
void apply_ddt(Field3D& f) override;
209+
210+
private:
211+
std::shared_ptr<FieldGenerator> gen;
212+
};
213+
166214
/// Neumann boundary condition set half way between guard cell and grid cell at 2nd order accuracy
167215
class BoundaryNeumann_2ndOrder : public BoundaryOp {
168216
public:

‎src/mesh/boundary_factory.cxx‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,14 +17,24 @@ using std::string;
1717
BoundaryFactory* BoundaryFactory::instance = nullptr;
1818

1919
BoundaryFactory::BoundaryFactory() {
20-
add(new BoundaryDirichlet(), "dirichlet");
20+
add(new BoundaryDirichlet_O1(), "dirichlet_o1"); // Old implementation in v3
21+
add(new BoundaryDirichlet_O1(), "dirichlet_O1");
22+
add(new BoundaryDirichlet(), "dirichlet"); // Default
2123
add(new BoundaryDirichlet(), "dirichlet_o2"); // Synonym for "dirichlet"
24+
add(new BoundaryDirichlet(), "dirichlet_O2"); // Synonym for "dirichlet"
2225
add(new BoundaryDirichlet_O3(), "dirichlet_o3");
26+
add(new BoundaryDirichlet_O3(), "dirichlet_O3");
2327
add(new BoundaryDirichlet_O4(), "dirichlet_o4");
28+
add(new BoundaryDirichlet_O4(), "dirichlet_O4");
2429
add(new BoundaryDirichlet_4thOrder(), "dirichlet_4thorder");
25-
add(new BoundaryNeumann(), "neumann");
30+
31+
add(new BoundaryNeumann_O1(), "neumann_o1"); // Old implementation in v3
32+
add(new BoundaryNeumann_O1(), "neumann_O1");
33+
add(new BoundaryNeumann(), "neumann"); // Default
34+
add(new BoundaryNeumann(), "neumann_o2"); // Synonym for "neumann"
2635
add(new BoundaryNeumann(), "neumann_O2"); // Synonym for "neumann"
2736
add(new BoundaryNeumann_4thOrder(), "neumann_4thorder");
37+
add(new BoundaryNeumann_O4(), "neumann_o4");
2838
add(new BoundaryNeumann_O4(), "neumann_O4");
2939
add(new BoundaryNeumannPar(), "neumannpar");
3040
add(new BoundaryNeumann_NonOrthogonal(), "neumann_nonorthogonal");

‎src/mesh/boundary_standard.cxx‎

Lines changed: 261 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,135 @@ void verifyNumPoints(BoundaryRegion* region, int ptsRequired) {
116116
void verifyNumPoints(BoundaryRegion*, int) {}
117117
#endif
118118

119+
120+
///////////////////////////////////////////////////////////////
121+
122+
BoundaryOp* BoundaryDirichlet_O1::clone(BoundaryRegion* region,
123+
const std::list<std::string>& args) {
124+
verifyNumPoints(region, 1);
125+
126+
std::shared_ptr<FieldGenerator> newgen;
127+
if (!args.empty()) {
128+
// First argument should be an expression
129+
newgen = FieldFactory::get()->parse(args.front());
130+
}
131+
return new BoundaryDirichlet_O1(region, newgen);
132+
}
133+
134+
void BoundaryDirichlet_O1::apply(Field2D& f) { BoundaryDirichlet_O1::apply(f, 0.); }
135+
136+
void BoundaryDirichlet_O1::apply(Field2D& f, BoutReal t) {
137+
// Set (at 1st order) the value at the grid cell to the guard cells.
138+
139+
Mesh* mesh = bndry->localmesh;
140+
ASSERT1(mesh == f.getMesh());
141+
bndry->first();
142+
143+
// Decide which generator to use
144+
std::shared_ptr<FieldGenerator> fg = gen;
145+
if (!fg) {
146+
fg = f.getBndryGenerator(bndry->location);
147+
}
148+
149+
BoutReal val = 0.0;
150+
151+
// Check for staggered grids
152+
153+
CELL_LOC loc = f.getLocation();
154+
if (mesh->StaggerGrids) {
155+
// Staggered
156+
throw BoutException("dirichlet_o1 BC is not implementated for staggered grids.");
157+
158+
} else {
159+
// Non-staggered, standard case
160+
for (; !bndry->isDone(); bndry->next1d()) {
161+
162+
if (fg) {
163+
val = fg->generate(Context(bndry, loc, t, mesh));
164+
}
165+
f(bndry->x, bndry->y) = val;
166+
167+
// Need to set second guard cell, as may be used for interpolation or upwinding derivatives
168+
// This is not very efficient. Both boundary cells can be treated in one loop.
169+
for (int i = 1; i < bndry->width; i++) {
170+
int xi = bndry->x + i * bndry->bx;
171+
int yi = bndry->y + i * bndry->by;
172+
f(xi, yi) = val;
173+
}
174+
175+
}
176+
}
177+
}
178+
179+
void BoundaryDirichlet_O1::apply(Field3D& f) { BoundaryDirichlet_O1::apply(f, 0.); }
180+
181+
void BoundaryDirichlet_O1::apply(Field3D& f, BoutReal t) {
182+
// Set (at 1st order) the value at the grid cell to the guard cells.
183+
184+
Mesh* mesh = bndry->localmesh;
185+
ASSERT1(mesh == f.getMesh());
186+
bndry->first();
187+
188+
// Decide which generator to use
189+
std::shared_ptr<FieldGenerator> fg = gen;
190+
if (!fg) {
191+
fg = f.getBndryGenerator(bndry->location);
192+
}
193+
194+
BoutReal val = 0.0;
195+
196+
// Check for staggered grids
197+
198+
CELL_LOC loc = f.getLocation();
199+
if (mesh->StaggerGrids) {
200+
// Staggered.
201+
throw BoutException("dirichlet_o1 BC is not implementated for staggered grids.");
202+
203+
} else {
204+
// Standard (non-staggered) case
205+
for (; !bndry->isDone(); bndry->next1d()) {
206+
for (int zk = 0; zk < mesh->LocalNz; zk++) {
207+
if (fg) {
208+
val = fg->generate(Context(bndry, zk, loc, t, mesh));
209+
}
210+
f(bndry->x, bndry->y, zk) = val;
211+
}
212+
213+
// This is not very efficient. Both boundary cells can be treated in one loop.
214+
for (int i = 1; i < bndry->width; i++) {
215+
// Set any other guard cells using the values on the cells
216+
int xi = bndry->x + i * bndry->bx;
217+
int yi = bndry->y + i * bndry->by;
218+
for (int zk = 0; zk < mesh->LocalNz; zk++) {
219+
if (fg) {
220+
val = fg->generate(Context(bndry, zk, loc, t, mesh));
221+
}
222+
f(xi, yi, zk) = val;
223+
}
224+
}
225+
}
226+
}
227+
}
228+
229+
void BoundaryDirichlet_O1::apply_ddt(Field2D& f) {
230+
Field2D* dt = f.timeDeriv();
231+
for (bndry->first(); !bndry->isDone(); bndry->next()) {
232+
(*dt)(bndry->x, bndry->y) = 0.; // Set time derivative to zero
233+
}
234+
}
235+
236+
void BoundaryDirichlet_O1::apply_ddt(Field3D& f) {
237+
Mesh* mesh = bndry->localmesh;
238+
ASSERT1(mesh == f.getMesh());
239+
Field3D* dt = f.timeDeriv();
240+
241+
for (bndry->first(); !bndry->isDone(); bndry->next()) {
242+
for (int z = 0; z < mesh->LocalNz; z++) {
243+
(*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero
244+
}
245+
}
246+
}
247+
119248
///////////////////////////////////////////////////////////////
120249

121250
BoundaryOp* BoundaryDirichlet::clone(BoundaryRegion* region,
@@ -1715,6 +1844,138 @@ void BoundaryNeumann_NonOrthogonal::apply(Field3D& f) {
17151844

17161845
///////////////////////////////////////////////////////////////
17171846

1847+
BoundaryOp* BoundaryNeumann_O1::clone(BoundaryRegion * region,
1848+
const std::list<std::string>& args) {
1849+
verifyNumPoints(region, 1);
1850+
std::shared_ptr<FieldGenerator> newgen = nullptr;
1851+
if (!args.empty()) {
1852+
// First argument should be an expression
1853+
newgen = FieldFactory::get()->parse(args.front());
1854+
}
1855+
return new BoundaryNeumann_O1(region, newgen);
1856+
}
1857+
1858+
void BoundaryNeumann_O1::apply(Field2D & f) { BoundaryNeumann_O1::apply(f, 0.); }
1859+
1860+
void BoundaryNeumann_O1::apply(Field2D & f, BoutReal t) {
1861+
// Set (at 1st order) the gradient/value at the grid cell to the guard cells.
1862+
1863+
1864+
#if not(BOUT_USE_METRIC_3D)
1865+
Mesh* mesh = bndry->localmesh;
1866+
ASSERT1(mesh == f.getMesh());
1867+
Coordinates* metric = f.getCoordinates();
1868+
1869+
bndry->first();
1870+
1871+
// Decide which generator to use
1872+
std::shared_ptr<FieldGenerator> fg = gen;
1873+
if (!fg) {
1874+
fg = f.getBndryGenerator(bndry->location);
1875+
}
1876+
1877+
BoutReal val = 0.0;
1878+
1879+
// Check for staggered grids
1880+
1881+
CELL_LOC loc = f.getLocation();
1882+
if (mesh->StaggerGrids) {
1883+
// Staggered.
1884+
throw BoutException("neumann_o1 BC is not implementated for staggered grids.");
1885+
1886+
} else {
1887+
// Non-staggered, standard case
1888+
1889+
for (bndry->first(); !bndry->isDone(); bndry->next1d()) {
1890+
BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y)
1891+
+ bndry->by * metric->dy(bndry->x, bndry->y);
1892+
1893+
if (fg) {
1894+
val = fg->generate(Context(bndry, loc, t, mesh));
1895+
}
1896+
1897+
f(bndry->x, bndry->y) =
1898+
f(bndry->x - bndry->bx, bndry->y - bndry->by) + delta * val;
1899+
if (bndry->width == 2) {
1900+
f(bndry->x + bndry->bx, bndry->y + bndry->by) = f(bndry->x, bndry->y) + delta * val;
1901+
}
1902+
}
1903+
}
1904+
#else
1905+
throw BoutException("Applying boundary condition 'neumann' to Field2D "
1906+
"not compatible with 3D metrics in all cases.");
1907+
#endif
1908+
}
1909+
1910+
void BoundaryNeumann_O1::apply(Field3D & f) { BoundaryNeumann_O1::apply(f, 0.); }
1911+
1912+
void BoundaryNeumann_O1::apply(Field3D & f, BoutReal t) {
1913+
Mesh* mesh = bndry->localmesh;
1914+
ASSERT1(mesh == f.getMesh());
1915+
Coordinates* metric = f.getCoordinates();
1916+
1917+
bndry->first();
1918+
1919+
// Decide which generator to use
1920+
std::shared_ptr<FieldGenerator> fg = gen;
1921+
if (!fg) {
1922+
fg = f.getBndryGenerator(bndry->location);
1923+
}
1924+
1925+
BoutReal val = 0.0;
1926+
1927+
// Check for staggered grids
1928+
1929+
CELL_LOC loc = f.getLocation();
1930+
if (mesh->StaggerGrids) {
1931+
// Staggered.
1932+
throw BoutException("neumann_o1 BC is not implementated for staggered grids.");
1933+
1934+
} else {
1935+
for (; !bndry->isDone(); bndry->next1d()) {
1936+
#if BOUT_USE_METRIC_3D
1937+
for (int zk = 0; zk < mesh->LocalNz; zk++) {
1938+
BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y, zk)
1939+
+ bndry->by * metric->dy(bndry->x, bndry->y, zk);
1940+
#else
1941+
BoutReal delta = bndry->bx * metric->dx(bndry->x, bndry->y)
1942+
+ bndry->by * metric->dy(bndry->x, bndry->y);
1943+
for (int zk = 0; zk < mesh->LocalNz; zk++) {
1944+
#endif
1945+
if (fg) {
1946+
val = fg->generate(Context(bndry, zk, loc, t, mesh));
1947+
}
1948+
f(bndry->x, bndry->y, zk) =
1949+
f(bndry->x - bndry->bx, bndry->y - bndry->by, zk) + delta * val;
1950+
if (bndry->width == 2) {
1951+
f(bndry->x + bndry->bx, bndry->y + bndry->by, zk) =
1952+
f(bndry->x, bndry->y, zk) + delta * val;
1953+
}
1954+
}
1955+
}
1956+
}
1957+
}
1958+
1959+
void BoundaryNeumann_O1::apply_ddt(Field2D & f) {
1960+
Field2D* dt = f.timeDeriv();
1961+
for (bndry->first(); !bndry->isDone(); bndry->next()) {
1962+
(*dt)(bndry->x, bndry->y) = 0.; // Set time derivative to zero
1963+
}
1964+
}
1965+
1966+
void BoundaryNeumann_O1::apply_ddt(Field3D & f) {
1967+
Mesh* mesh = bndry->localmesh;
1968+
ASSERT1(mesh == f.getMesh());
1969+
Field3D* dt = f.timeDeriv();
1970+
for (bndry->first(); !bndry->isDone(); bndry->next()) {
1971+
for (int z = 0; z < mesh->LocalNz; z++) {
1972+
(*dt)(bndry->x, bndry->y, z) = 0.; // Set time derivative to zero
1973+
}
1974+
}
1975+
}
1976+
1977+
///////////////////////////////////////////////////////////////
1978+
17181979
BoundaryOp* BoundaryNeumann::clone(BoundaryRegion * region,
17191980
const std::list<std::string>& args) {
17201981
verifyNumPoints(region, 1);

0 commit comments

Comments
 (0)