Skip to content

Commit b99de51

Browse files
einolaeinola
committed
Add open (Neumann zero) boundary conditions
Adapts the dirichletFromEdge function to be neumannFromEdge and adds two neumannZero functions - one for DGVectors and one for CGVectors. Both are inspired by the existing dirichletZero function. The neumannZero function must be called after advection - which is a bit worse, but moving it into the DGTransport class is compicated.
1 parent 327de33 commit b99de51

File tree

7 files changed

+113
-20
lines changed

7 files changed

+113
-20
lines changed

‎dynamics/src/CGDynamicsKernel.cpp‎

Lines changed: 52 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/*!
22
* @file CGDynamicsKernel.cpp
33
*
4-
* @date 14 Jan 2025
4+
* @date 21 Feb 2025
55
* @author Tim Spain <[email protected]>
66
*/
77

@@ -418,10 +418,61 @@ void CGDynamicsKernel<DGadvection>::dirichletZero(CGVector<CGdegree>& v) const
418418
}
419419
}
420420

421+
template <int DGadvection>
422+
void CGDynamicsKernel<DGadvection>::neumannZero(CGVector<CGdegree>& v) const
423+
{
424+
// the four segments bottom, right, top, left, are each processed in parallel
425+
for (size_t seg = 0; seg < 4; ++seg) {
426+
#pragma omp parallel for
427+
for (size_t i = 0; i < smesh->neumann[seg].size(); ++i) {
428+
429+
const size_t eid = smesh->neumann[seg][i];
430+
const size_t ix = eid % smesh->nx; // compute coordinates of element
431+
const size_t iy = eid / smesh->nx;
432+
433+
switch (seg) {
434+
case 0: // bottom <= top
435+
for (size_t j = 0; j < CGdegree + 1; ++j)
436+
v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0) = v(
437+
(iy + 1) * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0);
438+
break;
439+
case 1: // right <= left
440+
for (size_t j = 0; j < CGdegree + 1; ++j)
441+
v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + CGdegree
442+
+ (CGdegree * smesh->nx + 1) * j,
443+
0)
444+
= v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix
445+
+ (CGdegree * smesh->nx + 1) * j,
446+
0);
447+
break;
448+
case 2: // top <= bottom
449+
for (size_t j = 0; j < CGdegree + 1; ++j)
450+
v((iy + 1) * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0)
451+
= v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + j, 0);
452+
break;
453+
case 3: // left <= right
454+
for (size_t j = 0; j < CGdegree + 1; ++j)
455+
v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix
456+
+ (CGdegree * smesh->nx + 1) * j,
457+
0)
458+
= v(iy * CGdegree * (CGdegree * smesh->nx + 1) + CGdegree * ix + CGdegree
459+
+ (CGdegree * smesh->nx + 1) * j,
460+
0);
461+
break;
462+
default:
463+
std::cerr << "That should not have happened!" << std::endl;
464+
abort();
465+
}
466+
}
467+
}
468+
}
469+
421470
template <int DGadvection> void CGDynamicsKernel<DGadvection>::applyBoundaries()
422471
{
423472
dirichletZero(u);
424473
dirichletZero(v);
474+
neumannZero(u);
475+
neumannZero(v);
425476
// TODO Periodic boundary conditions.
426477
}
427478

‎dynamics/src/ParametricMesh.cpp‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,14 @@
11
/*!
22
* @file ParametricMesh.hpp
3-
* @date 14 Jan 2025
3+
* @date 21 Feb 2025
44
* @author Thomas Richter <[email protected]>
55
*/
66

77
#include "ParametricMesh.hpp"
88

99
#include "ParametricTools.hpp"
1010

11+
#include <algorithm>
1112
#include <fstream>
1213
#include <iostream>
1314

@@ -259,11 +260,11 @@ void ParametricMesh::dirichletFromMask()
259260
}
260261

261262
/*!
262-
* Add to the dirichlet arrays due to the domain edges according to an edge index.
263+
* Add to the neumann arrays due to the domain edges according to an edge index.
263264
*
264-
* @param edge index of the edge to add closed boundary conditions to.
265+
* @param edge index of the edge to add open boundary conditions to.
265266
*/
266-
void ParametricMesh::dirichletFromEdge(Edge edge)
267+
void ParametricMesh::neumannFromEdge(Edge edge)
267268
{
268269
// BOTTOM, RIGHT, TOP, LEFT
269270
const std::array<size_t, N_EDGE> start = { 0, nx - 1, nelements - nx, 0 };
@@ -272,10 +273,11 @@ void ParametricMesh::dirichletFromEdge(Edge edge)
272273

273274
for (size_t idx = start[edge]; idx < stop[edge]; idx += stride[edge]) {
274275
if (landmask[idx]) {
275-
dirichlet[edge].push_back(idx);
276+
neumann[edge].push_back(idx);
276277
}
277278
}
278-
sortDirichlet(edge);
279+
for (Edge edge : edges)
280+
std::sort(neumann[edge].begin(), neumann[edge].end());
279281
}
280282

281283
/*!

‎dynamics/src/include/BrittleCGDynamicsKernel.hpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/*!
22
* @file BrittleCGDynamicsKernel.hpp
33
*
4-
* @date 06 Dec 2024
4+
* @date 21 Feb 2025
55
* @author Tim Spain <[email protected]>
66
* @author Einar Ólason <[email protected]>
77
*/
@@ -12,7 +12,6 @@
1212
#include "CGDynamicsKernel.hpp"
1313

1414
#include "BBMParameters.hpp"
15-
#include "ParametricMap.hpp"
1615
#include "StressUpdateStep.hpp"
1716
#include "include/constants.hpp"
1817
#include <cmath>
@@ -37,6 +36,7 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
3736
using DynamicsKernel<DGadvection, DGstressComp>::applyBoundaries;
3837
using DynamicsKernel<DGadvection, DGstressComp>::advectionAndLimits;
3938
using DynamicsKernel<DGadvection, DGstressComp>::dgtransport;
39+
using DynamicsKernel<DGadvection, DGstressComp>::neumannZero;
4040

4141
using CGDynamicsKernel<DGadvection>::u;
4242
using CGDynamicsKernel<DGadvection>::v;
@@ -109,6 +109,8 @@ template <int DGadvection> class BrittleCGDynamicsKernel : public CGDynamicsKern
109109
Nextsim::LimitMax(damage, 1.0);
110110
Nextsim::LimitMin(damage, 1e-12);
111111

112+
neumannZero(damage);
113+
112114
prepareIteration({ { hiceName, hice }, { ciceName, cice } });
113115

114116
// The timestep for the brittle solver is the solver subtimestep

‎dynamics/src/include/CGDynamicsKernel.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/*!
22
* @file CGDynamicsKernel.hpp
33
*
4-
* @date 06 Dec 2024
4+
* @date 21 Feb 2025
55
* @author Tim Spain <[email protected]>
66
*/
77

@@ -69,6 +69,7 @@ class CGDynamicsKernel : public DynamicsKernel<DGadvection, DGstressComp> {
6969
protected:
7070
void addStressTensorCell(const size_t eid, const size_t cx, const size_t cy);
7171
void dirichletZero(CGVector<CGdegree>&) const;
72+
void neumannZero(CGVector<CGdegree>&) const;
7273
// CG ice velocity
7374
CGVector<CGdegree> u;
7475
CGVector<CGdegree> v;

‎dynamics/src/include/DynamicsKernel.hpp‎

Lines changed: 40 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/*!
22
* @file DynamicsKernel.hpp
33
*
4-
* @date Jan 5, 2024
4+
* @date 21 Feb 2025
55
* @author Tim Spain <[email protected]>
66
*/
77

@@ -53,9 +53,9 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
5353
smesh->RotatePoleToGreenland();
5454
smesh->landmaskFromModelArray(mask);
5555
smesh->dirichletFromMask();
56-
// TODO: handle periodic and open edges
57-
for (ParametricMesh::Edge edge : ParametricMesh::edges) {
58-
smesh->dirichletFromEdge(edge);
56+
// TODO: handle periodic edges
57+
for (const ParametricMesh::Edge edge : ParametricMesh::edges) {
58+
smesh->neumannFromEdge(edge);
5959
}
6060

6161
//! Initialize transport
@@ -180,6 +180,9 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
180180
Nextsim::LimitMax(cice, 1.0);
181181
Nextsim::LimitMin(cice, 0.0);
182182
Nextsim::LimitMin(hice, 0.0);
183+
184+
neumannZero(cice);
185+
neumannZero(hice);
183186
}
184187

185188
protected:
@@ -232,6 +235,39 @@ template <int DGadvection, int DGstress> class DynamicsKernel {
232235
*/
233236
virtual void prepareAdvection() = 0;
234237

238+
template <int DG> void neumannZero(DGVector<DG>& v) const
239+
{
240+
// the four segments bottom, right, top, left, are each processed in parallel
241+
for (size_t seg = 0; seg < 4; ++seg) {
242+
#pragma omp parallel for
243+
for (size_t i = 0; i < smesh->neumann[seg].size(); ++i) {
244+
245+
const size_t eid = smesh->neumann[seg][i];
246+
const size_t ix = eid % smesh->nx; // compute coordinates of element
247+
const size_t iy = eid / smesh->nx;
248+
249+
switch (seg) {
250+
case 0: // bottom <= top
251+
for (size_t j = 0; j < DG; ++j)
252+
v(iy * smesh->nx + ix, j) = v((iy + 1) * smesh->nx + ix, j);
253+
break;
254+
case 1: // right <= left
255+
for (size_t j = 0; j < DG; ++j)
256+
v(iy * smesh->nx + ix, j) = v(iy * smesh->nx + ix - 1, j);
257+
break;
258+
case 2: // top <= bottom
259+
for (size_t j = 0; j < DG; ++j)
260+
v(iy * smesh->nx + ix, j) = v((iy - 1) * smesh->nx + ix, j);
261+
break;
262+
case 3: // left <= right
263+
for (size_t j = 0; j < DG; ++j)
264+
v(iy * smesh->nx + ix, j) = v(iy * smesh->nx + ix + 1, j);
265+
break;
266+
}
267+
}
268+
}
269+
}
270+
235271
private:
236272
std::unordered_map<std::string, DGVector<DGadvection>> advectedFields;
237273

‎dynamics/src/include/ParametricMesh.hpp‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/*!
22
* @file ParametricMesh.hpp
3-
* @date 09 Jan 2025
3+
* @date 21 Feb 2025
44
* @author Thomas Richter <[email protected]>
55
*/
66

@@ -56,9 +56,9 @@ class ParametricMesh {
5656
Eigen::Matrix<Nextsim::FloatType, Eigen::Dynamic, 2> vertices; // stores the
5757

5858
/*!
59-
* Stores the Dirichlet boundary information
59+
* Stores the Dirichlet and Neumann boundary information
6060
*
61-
* Dirichlet-Information is stored in 4 vectors
61+
* Dirichlet- and Neumann-Information is stored in 4 vectors, e.g.
6262
* dirichlet[0]: List of elements that have a Dirichlet-Boundary on the bottom
6363
* dirichlet[1]: List of elements that have a Dirichlet-Boundary on the right
6464
* dirichlet[2]: List of elements that have a Dirichlet-Boundary on the top
@@ -67,6 +67,7 @@ class ParametricMesh {
6767
* Each of the vectors is processed in parallel
6868
*/
6969
std::array<std::vector<size_t>, 4> dirichlet;
70+
std::array<std::vector<size_t>, 4> neumann;
7071

7172
/*!
7273
* Periodic:
@@ -403,7 +404,7 @@ class ParametricMesh {
403404
*
404405
* @param edge index of the edge to add closed boundary conditions to.
405406
*/
406-
void dirichletFromEdge(Edge edge);
407+
void neumannFromEdge(Edge edge);
407408

408409
/*!
409410
* Sort all the dirichlet arrays, so the element indices are ordered.

‎dynamics/test/ParametricMesh_test.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
*
44
* @brief Test the ParametricMesh class, especially processing from ModelArray files.
55
*
6-
* @date Dec 15, 2023
6+
* @date 21 Feb 2025
77
* @author Tim Spain <[email protected]>
88
*/
99

@@ -122,7 +122,7 @@ TEST_CASE("Compare readmesh and landmask reading")
122122
// Dirichlet conditions
123123
fromArrays.dirichletFromMask();
124124
for (ParametricMesh::Edge edge : ParametricMesh::edges) {
125-
fromArrays.dirichletFromEdge(edge);
125+
fromArrays.neumannFromEdge(edge);
126126
}
127127
fromArrays.sortDirichlet();
128128
for (ParametricMesh::Edge edge : ParametricMesh::edges) {

0 commit comments

Comments
 (0)