Skip to content

Dirichlet bcond #499

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.

Sign up for GitHub

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

Open
pdziekan wants to merge 10 commits into
base: master
Choose a base branch
from
Open

Dirichlet bcond #499

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
521b164
inflow sandbox test
pdziekan Dec 22, 2025
7e7ffc7
fixed boundary conditions
pdziekan Dec 22, 2025
09fea5e
Merge branch 'master' of github.com:igfuw/libmpdataxx into dirichlet
pdziekan Jan 8, 2026
dc7459d
save_edge_val sharedmem bcond
pdziekan Feb 3, 2026
666ebf6
Apply suggestions from code review
pdziekan Feb 3, 2026
c1690c6
merge
pdziekan Feb 3, 2026
bed3cdd
dirichlet fixes
pdziekan Feb 4, 2026
8ab7ce0
prs solvers laplacian calc (with err_init): move set_edges before G m…
pdziekan Feb 26, 2026
3bac12f
Merge branch 'init_laplacian_edges_fix' into dirichlet
pdziekan Feb 26, 2026
38ebf3c
fixed bcond: use open if value was not set
pdziekan Mar 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions libmpdata++/bcond/cyclic_2d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,11 @@ namespace libmpdataxx
fill_halos_sclr(a, j);
}

void save_edge_vel(const arr_t &, const rng_t &) {}
void save_edge_vel(const arr_t &, const rng_t &) override {}

void set_edge_pres(arr_t &, const rng_t &, int) {}
void save_edge_val(const arr_t &, const arr_t &, const rng_t &) override {}

void set_edge_pres(arr_t &, const rng_t &, int) override {}

void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
{
Expand Down Expand Up @@ -129,9 +131,11 @@ namespace libmpdataxx
fill_halos_sclr(a, j);
}

void save_edge_vel(const arr_t &, const rng_t &) {}
void save_edge_vel(const arr_t &, const rng_t &) override {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &) override {}

void set_edge_pres(arr_t &, const rng_t &, int) {}
void set_edge_pres(arr_t &, const rng_t &, int) override {}

void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
{
Expand Down
4 changes: 4 additions & 0 deletions libmpdata++/bcond/cyclic_3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ namespace libmpdataxx

void save_edge_vel(const arr_t &, const rng_t &, const rng_t &) {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &, const rng_t &) {}

void set_edge_pres(arr_t &, const rng_t &, const rng_t &, int) {}

void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const rng_t &k, const bool ad = false)
Expand Down Expand Up @@ -130,6 +132,8 @@ namespace libmpdataxx

void save_edge_vel(const arr_t &, const rng_t &, const rng_t &) {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &, const rng_t &) {}

void set_edge_pres(arr_t &, const rng_t &, const rng_t &, int) {}

void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const rng_t &k, const bool ad = false)
Expand Down
12 changes: 11 additions & 1 deletion libmpdata++/bcond/detail/bcond_common.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ namespace libmpdataxx
{
using namespace arakawa_c;

enum bcond_e { null, cyclic, polar, open, rigid, remote, gndsky, custom };
enum bcond_e { null, cyclic, polar, open, rigid, remote, gndsky, fixed, custom };
enum drctn_e { left, rght };

template<
Expand Down Expand Up @@ -81,6 +81,11 @@ namespace libmpdataxx
assert(false && "bcond::save_edge_vel() called!");
};

virtual void save_edge_val(const arr_2d_t &edge, const arr_2d_t &val, const rng_t &rng)
{
assert(false && "bcond::save_edge_val() called!");
};

virtual void set_edge_pres(arr_2d_t &, const rng_t &, int)
{
assert(false && "bcond::set_edge() called!");
Expand Down Expand Up @@ -147,6 +152,11 @@ namespace libmpdataxx
assert(false && "bcond::save_edge_vel() called!");
};

virtual void save_edge_val(const arr_3d_t &, const arr_3d_t &, const rng_t &, const rng_t &)
{
assert(false && "bcond::save_edge_val() called!");
};

virtual void set_edge_pres(arr_3d_t &, const rng_t &, const rng_t &, int)
{
assert(false && "bcond::set_edge() called!");
Expand Down
107 changes: 107 additions & 0 deletions libmpdata++/bcond/fixed_2d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
// 2D fixed (Dirichlet) boundary conditions for libmpdata++
//
// licensing: GNU GPL v3
// copyright: University of Warsaw

#pragma once

#include <libmpdata++/bcond/detail/bcond_common.hpp>

namespace libmpdataxx
{
namespace bcond
{
template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
class bcond< real_t, halo, knd, dir, n_dims, d,
typename std::enable_if<
knd == fixed &&
dir == left &&
n_dims == 2
>::type
> : public bcond<real_t, halo, open, dir, n_dims, d> // TODO: in some cases (e.g. cloud chamber) we might want to inherit from rigid?
{
using parent_t = bcond<real_t, halo, open, dir, n_dims, d>;
using arr_t = blitz::Array<real_t, 2>;
using parent_t::parent_t; // inheriting ctor

// holds fixed wall values (equal to initial edge values)
std::unordered_map<const arr_t*, arr_t> edge_values;

public:

void fill_halos_sclr(arr_t &a, const rng_t &j, const bool deriv = false)
{
// #if !defined(NDEBUG)
// assert(edge_values.find(&a) != edge_values.end() && "fixed bcond: edge values not saved before filling halos");
// #endif

// if edge values are not saved, use open boundary (done to fix sgs eddy viscosity k_m)
if(edge_values.find(&a) == edge_values.end())
{
parent_t::fill_halos_sclr(a, j, deriv);
return;
}

using namespace idxperm;
for (int i = this->left_halo_sclr.first(); i <= this->left_halo_sclr.last(); ++i)
{
a(pi<d>(i, j)) = edge_values[&a](pi<d>(0, j));
}
}

void save_edge_val(const arr_t &a, const arr_t &val, const rng_t &j)
{
using namespace idxperm;
// assert(a.shape() == val.shape());
Copy link

Copilot AI Dec 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Commented-out assertion should be either enabled or removed. Line 48 contains a commented assertion to verify that array shapes match. If this check is needed, it should be enabled within an #if !defined(NDEBUG) block like the assertion on line 35. If the check is not needed, the commented code should be removed to reduce clutter and improve maintainability.

Copilot uses AI. Check for mistakes.
auto s = a.shape();
s[d] = 1;
edge_values.try_emplace(&a, arr_t(s));
if constexpr(d == 0) edge_values[&a].reindexSelf({0, a.lbound(1)});
else if constexpr(d == 1) edge_values[&a].reindexSelf({a.lbound(0), 0});
edge_values[&a](pi<d>(0, j)) = val(pi<d>(this->left_edge_sclr, j));
}
};

template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
class bcond< real_t, halo, knd, dir, n_dims, d,
typename std::enable_if<
knd == fixed &&
dir == rght &&
n_dims == 2
>::type
> : public bcond<real_t, halo, open, dir, n_dims, d>
{
using parent_t = bcond<real_t, halo, open, dir, n_dims, d>;
using arr_t = blitz::Array<real_t, 2>;
using parent_t::parent_t; // inheriting ctor

// holds fixed wall values (equal to initial edge values)
std::unordered_map<const arr_t*, arr_t> edge_values;

public:

// TODO

// void fill_halos_sclr(arr_t &a, const rng_t &j, const bool deriv = false)
// {
// using namespace idxperm;
// for (int i = this->rght_halo_sclr.first(); i <= this->rght_halo_sclr.last(); ++i)
// {
// a(pi<d>(i, j)) = edge_values[&a](pi<d>(0, j));
// }
// }

// void save_edge_val(const arr_t &a, const arr_t &val, const rng_t &j)
// {
// using namespace idxperm;
// // assert(a.shape() == val.shape());
// auto s = a.shape();
// s[d] = 1;
// edge_values.try_emplace(&a, arr_t(s));
// if constexpr (d == 0) edge_values[&a].reindexSelf({0, a.lbound(1)});
// else if constexpr(d == 1) edge_values[&a].reindexSelf({a.lbound(0), 0});
// edge_values[&a](pi<d>(0, j)) = val(pi<d>(this->rght_edge_sclr, j));
// }
};
} // namespace bcond
} // namespace libmpdataxx
105 changes: 105 additions & 0 deletions libmpdata++/bcond/fixed_3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
// 3D fixed (Dirichlet) boundary conditions for libmpdata++
//
// licensing: GNU GPL v3
// copyright: University of Warsaw

#pragma once

#include <libmpdata++/bcond/detail/bcond_common.hpp>

namespace libmpdataxx
{
namespace bcond
{
template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
class bcond< real_t, halo, knd, dir, n_dims, d,
typename std::enable_if<
knd == fixed &&
dir == left &&
n_dims == 3
>::type
> : public bcond<real_t, halo, open, dir, n_dims, d> // TODO: in some cases (e.g. cloud chamber) we might want to inherit from rigid?
{
using parent_t = bcond<real_t, halo, open, dir, n_dims, d>;
using arr_t = blitz::Array<real_t, 3>;
using parent_t::parent_t; // inheriting ctor

// holds fixed wall values (equal to initial edge values)
std::unordered_map<const arr_t*, arr_t> edge_values;

public:

void fill_halos_sclr(arr_t &a, const rng_t &j, const rng_t &k, const bool deriv = false)
{
#if !defined(NDEBUG)
assert(edge_values.find(&a) != edge_values.end() && "fixed bcond: edge values not saved before filling halos");
#endif

using namespace idxperm;
for (int i = this->left_halo_sclr.first(); i <= this->left_halo_sclr.last(); ++i)
{
a(pi<d>(i, j, k)) = edge_values[&a](pi<d>(0, j, k));
}
}

void save_edge_val(const arr_t &a, const arr_t &val, const rng_t &j, const rng_t &k)
{
using namespace idxperm;
// assert(a.shape() == val.shape());
Copy link

Copilot AI Dec 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Commented-out assertion should be either enabled or removed. Line 48 contains a commented assertion to verify that array shapes match. If this check is needed, it should be enabled within an #if !defined(NDEBUG) block like the assertion in the left direction implementation. If the check is not needed, the commented code should be removed to reduce clutter and improve maintainability.

Suggested change
// assert(a.shape() == val.shape());
#if !defined(NDEBUG)
assert(a.shape() == val.shape());
#endif

Copilot uses AI. Check for mistakes.
auto s = a.shape();
s[d] = 1;
edge_values.try_emplace(&a, arr_t(s));
if constexpr (d == 0) edge_values[&a].reindexSelf({0, a.lbound(1), a.lbound(2)});
else if constexpr (d == 1) edge_values[&a].reindexSelf({a.lbound(0), 0, a.lbound(2)});
else if constexpr (d == 2) edge_values[&a].reindexSelf({a.lbound(0), a.lbound(1), 0});
edge_values[&a](pi<d>(0, j, k)) = val(pi<d>(this->left_edge_sclr, j, k));
}
};


template <typename real_t, int halo, bcond_e knd, drctn_e dir, int n_dims, int d>
class bcond< real_t, halo, knd, dir, n_dims, d,
typename std::enable_if<
knd == fixed &&
dir == rght &&
n_dims == 3
>::type
> : public bcond<real_t, halo, open, dir, n_dims, d> // TODO: in some cases (e.g. cloud chamber) we might want to inherit from rigid?
{
using parent_t = bcond<real_t, halo, open, dir, n_dims, d>;
using arr_t = blitz::Array<real_t, 3>;
using parent_t::parent_t; // inheriting ctor

// holds fixed wall values (equal to initial edge values)
std::unordered_map<const arr_t*, arr_t> edge_values;

public:

void fill_halos_sclr(arr_t &a, const rng_t &j, const rng_t &k, const bool deriv = false)
{
#if !defined(NDEBUG)
assert(edge_values.find(&a) != edge_values.end() && "fixed bcond: edge values not saved before filling halos");
#endif

using namespace idxperm;
for (int i = this->rght_halo_sclr.first(); i <= this->rght_halo_sclr.last(); ++i)
{
a(pi<d>(i, j, k)) = edge_values[&a](pi<d>(0, j, k));
}
}

void save_edge_val(const arr_t &a, const arr_t &val, const rng_t &j, const rng_t &k)
{
using namespace idxperm;
// assert(a.shape() == val.shape());
Copy link

Copilot AI Dec 22, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Commented-out assertion should be either enabled or removed. Line 94 contains a commented assertion to verify that array shapes match. If this check is needed, it should be enabled within an #if !defined(NDEBUG) block like the assertion in the left direction implementation. If the check is not needed, the commented code should be removed to reduce clutter and improve maintainability.

Suggested change
// assert(a.shape() == val.shape());
#if !defined(NDEBUG)
assert(a.shape() == val.shape());
#endif

Copilot uses AI. Check for mistakes.
auto s = a.shape();
s[d] = 1;
edge_values.try_emplace(&a, arr_t(s));
if constexpr (d == 0) edge_values[&a].reindexSelf({0, a.lbound(1), a.lbound(2)});
else if constexpr (d == 1) edge_values[&a].reindexSelf({a.lbound(0), 0, a.lbound(2)});
else if constexpr (d == 2) edge_values[&a].reindexSelf({a.lbound(0), a.lbound(1), 0});
edge_values[&a](pi<d>(0, j, k)) = val(pi<d>(this->rght_edge_sclr, j, k));
}
};
} // namespace bcond
} // namespace libmpdataxx
4 changes: 4 additions & 0 deletions libmpdata++/bcond/open_2d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,8 @@ namespace libmpdataxx
edge_velocity(pi<d>(0, j)) = a(pi<d>(this->left_edge_sclr, j));
}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &) {}

void set_edge_pres(arr_t &a, const rng_t &j, int sign)
{
using namespace idxperm;
Expand Down Expand Up @@ -158,6 +160,8 @@ namespace libmpdataxx
edge_velocity(pi<d>(0, j)) = a(pi<d>(this->rght_edge_sclr, j));
}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &) {}

void set_edge_pres(arr_t &a, const rng_t &j, int sign)
{
using namespace idxperm;
Expand Down
4 changes: 4 additions & 0 deletions libmpdata++/bcond/open_3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,8 @@ namespace libmpdataxx
edge_velocity(pi<d>(0, j, k)) = a(pi<d>(this->left_edge_sclr, j, k));
}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &, const rng_t &) {}

void set_edge_pres(arr_t &a, const rng_t &j, const rng_t &k, int sign)
{
using namespace idxperm;
Expand Down Expand Up @@ -187,6 +189,8 @@ namespace libmpdataxx
edge_velocity(pi<d>(0, j, k)) = a(pi<d>(this->rght_edge_sclr, j, k));
}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &, const rng_t &) {}

void set_edge_pres(arr_t &a, const rng_t &j, const rng_t &k, int sign)
{
using namespace idxperm;
Expand Down
4 changes: 4 additions & 0 deletions libmpdata++/bcond/remote_2d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,8 @@ namespace libmpdataxx

void save_edge_vel(const arr_t &, const rng_t &) {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &) {}

void set_edge_pres(arr_t &, const rng_t &, int) {}

// we require that given process calculates its internal vectors + the nearest vector to the left, TODO: how to enforce this?
Expand Down Expand Up @@ -159,6 +161,8 @@ namespace libmpdataxx

void save_edge_vel(const arr_t &, const rng_t &) {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &) {}

void set_edge_pres(arr_t &, const rng_t &, int) {}

void fill_halos_vctr_alng(arrvec_t<arr_t> &av, const rng_t &j, const bool ad = false)
Expand Down
4 changes: 4 additions & 0 deletions libmpdata++/bcond/remote_3d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,8 @@ namespace libmpdataxx

void save_edge_vel(const arr_t &, const rng_t &, const rng_t &) {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &, const rng_t &) {}

void set_edge_pres(arr_t &, const rng_t &, const rng_t &, int) {}


Expand Down Expand Up @@ -220,6 +222,8 @@ namespace libmpdataxx

void save_edge_vel(const arr_t &, const rng_t &, const rng_t &) {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &, const rng_t &) {}

void set_edge_pres(arr_t &, const rng_t &, const rng_t &, int) {}


Expand Down
8 changes: 6 additions & 2 deletions libmpdata++/bcond/rigid_2d.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,11 @@ namespace libmpdataxx
}
}

void save_edge_vel(const arr_t &, const rng_t &) {}
void save_edge_vel(const arr_t &, const rng_t &) override {}

void set_edge_pres(arr_t &a, const rng_t &j, int)
void save_edge_val(const arr_t &, const arr_t &, const rng_t &) override {}

void set_edge_pres(arr_t &a, const rng_t &j, int) override
{
using namespace idxperm;
a(pi<d>(this->left_edge_sclr, j)) = 0;
Expand Down Expand Up @@ -145,6 +147,8 @@ namespace libmpdataxx

void save_edge_vel(const arr_t &, const rng_t &) {}

void save_edge_val(const arr_t &, const arr_t &, const rng_t &) {}

void set_edge_pres(arr_t &a, const rng_t &j, int)
{
using namespace idxperm;
Expand Down
Loading
Loading