diff --git a/cli/new_command.py b/cli/new_command.py index 90df722b2..ab06091a4 100755 --- a/cli/new_command.py +++ b/cli/new_command.py @@ -24,7 +24,7 @@ def ValidateSimName(sim_name): - pattern = re.compile("^[a-zA-Z]+[a-zA-Z0-9\-_]+$") + pattern = re.compile(r"^[a-zA-Z]+[a-zA-Z0-9\-_]+$") if not pattern.match(sim_name): Print.error('Error: simulation name "{}" is not valid.'.format(sim_name)) Print.error(" Allowed characters are a-z A-Z 0-9 - and _") diff --git a/cli/util.py b/cli/util.py index 5c9762178..53095617a 100755 --- a/cli/util.py +++ b/cli/util.py @@ -23,5 +23,5 @@ def GetBinaryName(): with open("CMakeLists.txt") as f: content = f.read() - return re.search("project\((.*)\)", content).group(1) + return re.search(r"project\((.*)\)", content).group(1) diff --git a/cmake/external/ROOT.cmake b/cmake/external/ROOT.cmake index 4be6d2b11..a137b295f 100644 --- a/cmake/external/ROOT.cmake +++ b/cmake/external/ROOT.cmake @@ -13,7 +13,11 @@ if(APPLE) "${DETECTED_OS_VERS}" MATCHES "^osx-11.7") execute_process(COMMAND bash "-c" "xcodebuild -version | sed -En 's/Xcode[[:space:]]+([0-9\.]*)/\\1/p'" OUTPUT_VARIABLE XCODE_VERS) message(STATUS "##### XCODE version: ${XCODE_VERS}") - if("${XCODE_VERS}" GREATER_EQUAL "15.0") + if("${XCODE_VERS}" GREATER_EQUAL "15.3") + message(STATUS "##### Using ROOT builds for XCODE 15.3") + set(ROOT_TAR_FILE root_v6.30.06_cxx17_python3.9_osx-xcode-15.3-${DETECTED_ARCH}.tar.gz) + set(ROOT_SHA_KEY osx-xcode-15.3-${DETECTED_ARCH}-ROOT) + elseif("${XCODE_VERS}" GREATER_EQUAL "15.0") message(STATUS "##### Using ROOT builds for XCODE 15.2") set(ROOT_TAR_FILE root_v6.30.02_cxx17_python3.9_osx-xcode-15.2-${DETECTED_ARCH}.tar.gz) set(ROOT_SHA_KEY osx-xcode-15.2-${DETECTED_ARCH}-ROOT) diff --git a/cmake/external/SHA256Digests.cmake b/cmake/external/SHA256Digests.cmake index a9274759c..3dadffbd1 100644 --- a/cmake/external/SHA256Digests.cmake +++ b/cmake/external/SHA256Digests.cmake @@ -58,6 +58,10 @@ SET(osx-14.2-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e8055 SET(osx-14.2-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4) SET(osx-14.3-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e80554f7fa617ccb) SET(osx-14.3-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4) +SET(osx-14.4-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e80554f7fa617ccb) +SET(osx-14.4-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4) +SET(osx-14.5-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e80554f7fa617ccb) +SET(osx-14.5-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4) SET(osx-xcode-13.1-i386-ROOT be97dd72022c8d082fbe4394f18b55c4920f20b138cfff1b5fc2b41d397ac203) SET(osx-xcode-13.1-arm64-ROOT 0a55b91c6df42d152b7943912e134f05c3872a73e73fcc129ee87fd847240ec8) SET(osx-xcode-14.1-i386-ROOT 001311608512b24535bb9710b8baf006bd00e9b0595fd6bdf900d28b1e22c395) @@ -68,6 +72,8 @@ SET(osx-xcode-14.3-i386-ROOT c4f4c7a9848121ae19f3fe3c40c5812813ed42d9dfa7f1948f6 SET(osx-xcode-14.3-arm64-ROOT 6ffcfaadd4de627cffde2e42fded9175b730769071396d7abdabeefc2714839a) SET(osx-xcode-15.2-i386-ROOT 3dea500451dc2c84d150a7d5f4b27d2cba3e867d7c775b744ac898f9f6e51249) SET(osx-xcode-15.2-arm64-ROOT 677429937f13a44a088b71851f0af111509bf81a9a25a299c905b6c9f7249587) +SET(osx-xcode-15.3-i386-ROOT 52cef545acf679c3e491c9ef14e18edee454c90309b7fcf5f36fc2e659bdbde5) +SET(osx-xcode-15.3-arm64-ROOT 3ee0885c329fab9a965c278cc491db82dcbde3999a96ac80d6f5145fb6794d1e) SET(ubuntu-18.04-Libroadrunner bf9293c1c95d0b65227bd7e08c0048116ba851bcec5028ef72ea13762ef79276) SET(ubuntu-18.04-ParaView e3fd74b13e9a4086988c5104c6b8d95c56365d25b491706a8e72018d0e5c76f1) diff --git a/src/core/diffusion/diffusion_grid.cc b/src/core/diffusion/diffusion_grid.cc index 065c173e2..bacec5135 100644 --- a/src/core/diffusion/diffusion_grid.cc +++ b/src/core/diffusion/diffusion_grid.cc @@ -35,6 +35,8 @@ std::string BoundaryTypeToString(const BoundaryConditionType& type) { return "Dirichlet"; case BoundaryConditionType::kPeriodic: return "Periodic"; + case BoundaryConditionType::kComplex: + return "Complex"; default: return "unknown"; } @@ -52,6 +54,8 @@ BoundaryConditionType StringToBoundaryType(const std::string& type) { return BoundaryConditionType::kDirichlet; } else if (type == "Periodic") { return BoundaryConditionType::kPeriodic; + } else if (type == "Complex") { + return BoundaryConditionType::kComplex; } else { Log::Fatal("StringToBoundaryType", "Unknown boundary type: ", type); return BoundaryConditionType::kNeumann; @@ -64,7 +68,9 @@ DiffusionGrid::DiffusionGrid(int substance_id, : dc_({{1 - dc, dc / 6, dc / 6, dc / 6, dc / 6, dc / 6, dc / 6}}), mu_(mu), resolution_(resolution), - boundary_condition_(std::make_unique(0.0)) { + boundary_condition_(std::make_unique(0.0)), + boundary_condition_neumann_(std::make_unique(0.0)), + boundary_condition_dirichlet_(std::make_unique(0.0)) { // Compatibility with new abstract interface SetContinuumId(substance_id); SetContinuumName(substance_name); @@ -81,6 +87,24 @@ void DiffusionGrid::Initialize() { "')"); } + if (bc_type_ == BoundaryConditionType::kComplex) { + auto* param = Simulation::GetActive()->GetParam(); + complex_bct_ = {param->diffusion_bc_x_min, + param->diffusion_bc_x_max, + param->diffusion_bc_y_min, + param->diffusion_bc_y_max, + param->diffusion_bc_z_min, + param->diffusion_bc_z_max}; + for (int i=0; i<6; i++) { + if (complex_bct_[i] != "Neumann" && complex_bct_[i] != "Dirichlet" && + complex_bct_[i] != "Periodic") { + Log::Fatal("DiffusionGrid::Initialize", + "Complex boundaries: unknown boundary condition ", complex_bct_[i], + ". Select 'Neumann', 'Dirichlet' or 'Periodic' for each side of boundary."); + } + } + } + // Get neighbor grid dimensions auto* env = Simulation::GetActive()->GetEnvironment(); auto bounds = env->GetDimensionThresholds(); @@ -145,6 +169,8 @@ void DiffusionGrid::Diffuse(real_t dt) { DiffuseWithNeumann(dt); } else if (bc_type_ == BoundaryConditionType::kPeriodic) { DiffuseWithPeriodic(dt); + } else if (bc_type_ == BoundaryConditionType::kComplex) { + DiffuseWithComplex(dt); } else { Log::Error( "EulerGrid::Diffuse", "Boundary condition of type '", @@ -533,6 +559,16 @@ void DiffusionGrid::SetBoundaryCondition( boundary_condition_ = std::move(bc); } +void DiffusionGrid::SetBoundaryCondition_neumann( + std::unique_ptr bc) { + boundary_condition_neumann_ = std::move(bc); +} + +void DiffusionGrid::SetBoundaryCondition_dirichlet( + std::unique_ptr bc) { + boundary_condition_dirichlet_ = std::move(bc); +} + BoundaryCondition* DiffusionGrid::GetBoundaryCondition() const { return boundary_condition_.get(); } diff --git a/src/core/diffusion/diffusion_grid.h b/src/core/diffusion/diffusion_grid.h index 0b5e74c18..4a83fd79f 100644 --- a/src/core/diffusion/diffusion_grid.h +++ b/src/core/diffusion/diffusion_grid.h @@ -37,7 +37,8 @@ enum class BoundaryConditionType { kNeumann, kOpenBoundaries, kClosedBoundaries, - kPeriodic + kPeriodic, + kComplex }; enum class InteractionMode { kAdditive = 0, kExponential = 1, kLogistic = 2 }; @@ -121,6 +122,11 @@ class DiffusionGrid : public ScalarField { /// Compute the diffusion of the substance with Periodic boundary conditions /// (u_{-1} = {u_N-1}, u_{N} = u_{0}) for a given time step dt. virtual void DiffuseWithPeriodic(real_t dt) = 0; + /// Compute the diffusion of the substance with Complex boundary conditions, + /// where the boundary type can be Dirichlet, Neumann or Periodic depending on + /// the side of the environment. + virtual void DiffuseWithComplex(real_t dt) = 0; + /// Calculates the gradient for each box in the diffusion grid. /// The gradient is calculated in each direction (x, y, z) as following: @@ -323,6 +329,10 @@ class DiffusionGrid : public ScalarField { /// `ConstantBoundaryCondition` void SetBoundaryCondition(std::unique_ptr bc); + void SetBoundaryCondition_neumann(std::unique_ptr bc); + + void SetBoundaryCondition_dirichlet(std::unique_ptr bc); + /// @brief Returns the boundary condition. Does not transfer ownership. /// @return const Pointer to the boundary condition BoundaryCondition* GetBoundaryCondition() const; @@ -415,6 +425,14 @@ class DiffusionGrid : public ScalarField { BoundaryConditionType bc_type_ = BoundaryConditionType::kDirichlet; /// Object that implements the boundary conditions. std::unique_ptr boundary_condition_ = nullptr; + /// Object that implements the boundary condition for Neumann type when + /// using complex boundaries. + std::unique_ptr boundary_condition_neumann_ = nullptr; + /// Object that implements the boundary condition for Dirichlet type when + /// using complex boundaries. + std::unique_ptr boundary_condition_dirichlet_ = nullptr; + /// Stores the boundary condition types when using complex boundaries. + std::array complex_bct_ = {{""}}; /// Flag to indicate if the grid is initialized bool initialized_ = false; /// Flag to indicate if we want to print information about the grid after diff --git a/src/core/diffusion/euler_grid.cc b/src/core/diffusion/euler_grid.cc index f8ac1b722..e4e2e2fe9 100644 --- a/src/core/diffusion/euler_grid.cc +++ b/src/core/diffusion/euler_grid.cc @@ -366,4 +366,150 @@ void EulerGrid::DiffuseWithPeriodic(real_t dt) { c1_.swap(c2_); } + + + +void EulerGrid::DiffuseWithComplex(real_t dt) { + + const auto nx = resolution_; + const auto ny = resolution_; + const auto nz = resolution_; + const auto num_boxes = nx * ny * nz; + + const real_t d = 1 - dc_[0]; + const real_t dx = box_length_; + + const auto sim_time = GetSimulatedTime(); + + constexpr size_t YBF = 16; +#pragma omp parallel for collapse(2) + for (size_t yy = 0; yy < ny; yy += YBF) { + for (size_t z = 0; z < nz; z++) { + size_t ymax = yy + YBF; + if (ymax >= ny) { + ymax = ny; + } + for (size_t y = yy; y < ymax; y++) { + size_t x{0}; + size_t c{0}; + size_t l{0}; + size_t r{0}; + size_t n{0}; + size_t s{0}; + size_t b{0}; + size_t t{0}; + c = x + y * nx + z * nx * ny; +#pragma omp simd + for (x = 0; x < nx; x++) { + l = c - 1; + r = c + 1; + n = c - nx; + s = c + nx; + b = c - nx * ny; + t = c + nx * ny; + + // Clamp to avoid out of bounds access. Clamped values are initialized + // to a wrong value but will be overwritten by the boundary condition + // evaluation. All other values are correct. + real_t left{c1_[std::clamp(l, size_t{0}, num_boxes - 1)]}; + real_t right{c1_[std::clamp(r, size_t{0}, num_boxes - 1)]}; + real_t bottom{c1_[std::clamp(b, size_t{0}, num_boxes - 1)]}; + real_t top{c1_[std::clamp(t, size_t{0}, num_boxes - 1)]}; + real_t north{c1_[std::clamp(n, size_t{0}, num_boxes - 1)]}; + real_t south{c1_[std::clamp(s, size_t{0}, num_boxes - 1)]}; + + real_t CF_x{2.0}; + real_t CF_y{2.0}; + real_t CF_z{2.0}; + + if (x == 0 || x == (nx - 1) || y == 0 || y == (ny - 1) || z == 0 || + z == (nz - 1)) { + + real_t real_x = grid_dimensions_[0] + x * box_length_; + real_t real_y = grid_dimensions_[2] + y * box_length_; + real_t real_z = grid_dimensions_[4] + z * box_length_; + + real_t boundary_value_neumann = boundary_condition_neumann_->Evaluate(real_x, real_y, real_z, sim_time); + real_t boundary_value_dirichlet = boundary_condition_dirichlet_->Evaluate(real_x, real_y, real_z, sim_time); + + if (x == 0) { + if (complex_bct_[0] == "Dirichlet") { + left = boundary_value_dirichlet; + } else if (complex_bct_[0] == "Neumann") { + left = - boundary_value_neumann * box_length_; + CF_x -= 1.0; + } else if (complex_bct_[0] == "Periodic") { + l = nx-1 + y * nx + z * nx * ny; + left = c1_[l]; + } + } else if (x == (nx - 1)) { + if (complex_bct_[1] == "Dirichlet") { + right = boundary_value_dirichlet; + } else if (complex_bct_[1] == "Neumann") { + right = - boundary_value_neumann * box_length_; + CF_x -= 1.0; + } else if (complex_bct_[1] == "Periodic") { + r = 0 + y * nx + z * nx * ny; + right = c1_[r]; + } + } + + if (y == 0) { + if (complex_bct_[2] == "Dirichlet") { + north = boundary_value_dirichlet; + } else if (complex_bct_[2] == "Neumann") { + north = - boundary_value_neumann * box_length_; + CF_y -= 1.0; + } else if (complex_bct_[2] == "Periodic") { + n = x + (ny-1) * nx + z * nx * ny; + north = c1_[n]; + } + } else if (y == (ny - 1)) { + if (complex_bct_[3] == "Dirichlet") { + south = boundary_value_dirichlet; + } else if (complex_bct_[3] == "Neumann") { + south = - boundary_value_neumann * box_length_; + CF_y -= 1.0; + } else if (complex_bct_[3] == "Periodic") { + s = x + 0 * nx + z * nx * ny; + south = c1_[s]; + } + } + + if (z == 0) { + if (complex_bct_[4] == "Dirichlet") { + bottom = boundary_value_dirichlet; + } else if (complex_bct_[4] == "Neumann") { + bottom = - boundary_value_neumann * box_length_; + CF_z -= 1.0; + } else if (complex_bct_[4] == "Periodic") { + b = x + y * nx + (nz-1) * nx * ny; + bottom = c1_[b]; + } + } else if (z == (nz - 1)) { + if (complex_bct_[5] == "Dirichlet") { + top = boundary_value_dirichlet; + } else if (complex_bct_[5] == "Neumann") { + top = - boundary_value_neumann * box_length_; + CF_z -= 1.0; + } else if (complex_bct_[5] == "Periodic") { + t = x + y * nx + 0 * nx * ny; + top = c1_[t]; + } + } + } + + c2_[c] = c1_[c] * (1 - (mu_ * dt)) + ( (d * dt / (dx*dx) ) * (left + right - CF_x*c1_[c]) ) + + ( (d * dt / (dx*dx) ) * (south + north - CF_y*c1_[c]) ) + + ( (d * dt / (dx*dx) ) * (top + bottom - CF_z*c1_[c]) ); + + ++c; + } + } // tile ny + } // tile nz + } // block ny + c1_.swap(c2_); +} + + } // namespace bdm diff --git a/src/core/diffusion/euler_grid.h b/src/core/diffusion/euler_grid.h index 5052ae49d..5cd5aa6ee 100644 --- a/src/core/diffusion/euler_grid.h +++ b/src/core/diffusion/euler_grid.h @@ -53,6 +53,7 @@ class EulerGrid : public DiffusionGrid { void DiffuseWithDirichlet(real_t dt) override; void DiffuseWithNeumann(real_t dt) override; void DiffuseWithPeriodic(real_t dt) override; + void DiffuseWithComplex(real_t dt) override; private: BDM_CLASS_DEF_OVERRIDE(EulerGrid, 1); diff --git a/src/core/param/param.cc b/src/core/param/param.cc index bdcecb7dc..b2a3b8228 100644 --- a/src/core/param/param.cc +++ b/src/core/param/param.cc @@ -239,6 +239,12 @@ void Param::AssignFromConfig(const std::shared_ptr& config) { BDM_ASSIGN_CONFIG_VALUE(max_bound, "simulation.max_bound"); BDM_ASSIGN_CONFIG_VALUE(diffusion_boundary_condition, "simulation.diffusion_boundary_condition"); + BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_x_min, "simulation.diffusion_bc_x_min"); + BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_x_max, "simulation.diffusion_bc_x_max"); + BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_y_min, "simulation.diffusion_bc_y_min"); + BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_y_max, "simulation.diffusion_bc_y_max"); + BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_z_min, "simulation.diffusion_bc_z_min"); + BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_z_max, "simulation.diffusion_bc_z_max"); BDM_ASSIGN_CONFIG_VALUE(diffusion_method, "simulation.diffusion_method"); BDM_ASSIGN_CONFIG_VALUE(calculate_gradients, "simulation.calculate_gradients"); diff --git a/src/core/param/param.h b/src/core/param/param.h index f1d394ede..1cba5bc9b 100644 --- a/src/core/param/param.h +++ b/src/core/param/param.h @@ -243,6 +243,22 @@ struct Param { /// diffusion_boundary_condition = "Neumann" std::string diffusion_boundary_condition = "Neumann"; + /// Define the boundary condition of each side of the simulation space + /// when using complex boundaries [Neumann, Dirichlet, Periodic]. + /// Default value: `"Neumann"`\n TOML config file: + /// + /// [simulation] + /// diffusion_bc_x_min = "Neumann" + /// diffusion_bc_x_max = "Neumann" + /// .... + std::string diffusion_bc_x_min = "Neumann"; + std::string diffusion_bc_x_max = "Neumann"; + std::string diffusion_bc_y_min = "Neumann"; + std::string diffusion_bc_y_max = "Neumann"; + std::string diffusion_bc_z_min = "Neumann"; + std::string diffusion_bc_z_max = "Neumann"; + + /// A string for determining diffusion type within the simulation space. /// Currently, only the method "euler" implementing a FTCS scheme is /// supported. See for instance here: diff --git a/test/unit/core/diffusion_init_test.h b/test/unit/core/diffusion_init_test.h index 7ac3bec88..bc81cc43b 100644 --- a/test/unit/core/diffusion_init_test.h +++ b/test/unit/core/diffusion_init_test.h @@ -32,6 +32,7 @@ class TestGrid : public DiffusionGrid { void DiffuseWithDirichlet(real_t dt) override { return; }; void DiffuseWithNeumann(real_t dt) override { return; }; void DiffuseWithPeriodic(real_t dt) override { return; }; + void DiffuseWithComplex(real_t dt) override { return; }; void Swap() { c1_.swap(c2_); } diff --git a/test/unit/core/diffusion_test.cc b/test/unit/core/diffusion_test.cc index 7a9d78fe8..4ad6c89ab 100644 --- a/test/unit/core/diffusion_test.cc +++ b/test/unit/core/diffusion_test.cc @@ -970,6 +970,215 @@ TEST(DiffusionTest, DepletionMissmatch) { "that of the depleted one (ECM)*"); } +TEST(DiffusionTest, EulerComplexBoundariesUnknownBounds) { + ASSERT_DEATH( + { + auto set_param = [&](Param* param) { + param->bound_space = Param::BoundSpaceMode::kClosed; + param->diffusion_boundary_condition = "Complex"; + // Initialisation should fail when given boundary that is + // not 'Neumann', 'Dirichlet' or 'Periodic'. + param->diffusion_bc_z_max = "Neum"; + }; + Simulation simulation(TEST_NAME, set_param); + simulation.GetEnvironment()->Update(); + + auto* dgrid = new EulerGrid(0, "Kalium", 10, 0.0, 10); + dgrid->Initialize(); + }, + ".*Complex boundaries: unknown boundary condition*"); +} + +TEST(DiffusionTest, EulerComplexBoundariesCase1) { + double simulation_time_step{0.1}; + auto set_param = [&](Param* param) { + param->bound_space = Param::BoundSpaceMode::kClosed; + param->min_bound = -50; + param->max_bound = 50; + param->diffusion_boundary_condition = "Complex"; + param->diffusion_bc_x_min = "Dirichlet"; + param->diffusion_bc_x_max = "Neumann"; + param->diffusion_bc_y_min = "Periodic"; + param->diffusion_bc_y_max = "Periodic"; + param->diffusion_bc_z_min = "Periodic"; + param->diffusion_bc_z_max = "Periodic"; + }; + Simulation simulation(TEST_NAME, set_param); + simulation.GetEnvironment()->Update(); + auto* rm = simulation.GetResourceManager(); + + double decay_coef = 0.0; + double diff_coef = 400.0; + int res = 5; + double init = 1e5; + + // Dirichlet boundaries on one side and closed Neumann boundaries on the opposite. + auto* dgrid = new EulerGrid(0, "Kalium", diff_coef, decay_coef, res); + dgrid->Initialize(); + dgrid->SetBoundaryCondition_neumann(std::make_unique(0.0)); + dgrid->SetBoundaryCondition_dirichlet(std::make_unique(init)); + dgrid->SetUpperThreshold(1e15); + rm->AddContinuum(dgrid); + + double close = 0.0; + double far = 0.0; + + int tot = 2500; + for (int i=0; iDiffuse(simulation_time_step); + // We should a higher concentration closer to the dirichlet boundary than the + // neumann boundary at the start of the simulation. + if (i < 50) { + close = std::abs(dgrid->GetValue({-40, 0, 0})); + far = std::abs(dgrid->GetValue({40, 0, 0})); + EXPECT_LT(far, close); + } + } + + // After some time, the concentration across the space should become homogenous. + double final_close = std::abs(dgrid->GetValue({-40, 0, 0})); + double final_far = std::abs(dgrid->GetValue({40, 0, 0})); + EXPECT_LT((final_close - final_far) / final_close, 0.0001); + + rm->RemoveContinuum(0); +} + +TEST(DiffusionTest, EulerComplexBoundariesCase2) { + double simulation_time_step{0.1}; + auto set_param = [&](Param* param) { + param->bound_space = Param::BoundSpaceMode::kClosed; + param->min_bound = -50; + param->max_bound = 50; + param->diffusion_boundary_condition = "Complex"; + param->diffusion_bc_x_min = "Neumann"; + param->diffusion_bc_x_max = "Neumann"; + param->diffusion_bc_y_min = "Periodic"; + param->diffusion_bc_y_max = "Periodic"; + param->diffusion_bc_z_min = "Dirichlet"; + param->diffusion_bc_z_max = "Dirichlet"; + }; + Simulation simulation(TEST_NAME, set_param); + simulation.GetEnvironment()->Update(); + auto* rm = simulation.GetResourceManager(); + + double decay_coef = 0.0; + double diff_coef = 100.0; + int res = 5; + double init = 1e5; + + // Test 6 sources; one at each boundary of the simulation. + std::vector sources; + sources.push_back({40, 0, 0}); + sources.push_back({-40, 0, 0}); + sources.push_back({0, 40, 0}); + sources.push_back({0, -40, 0}); + sources.push_back({0, 0, 40}); + sources.push_back({0, 0, -40}); + + // Measure box adjacent to source within the simulation. + std::vector inside; + inside.push_back({20, 0, 0}); + inside.push_back({-20, 0, 0}); + inside.push_back({0, 20, 0}); + inside.push_back({0, -20, 0}); + inside.push_back({0, 0, 20}); + inside.push_back({0, 0, -20}); + + // Measure box adjacent to source on opposite side of simulation. + std::vector outside; + outside.push_back({-40, 0, 0}); + outside.push_back({40, 0, 0}); + outside.push_back({0, -40, 0}); + outside.push_back({0, 40, 0}); + outside.push_back({0, 0, -40}); + outside.push_back({0, 0, 40}); + + for (size_t s = 0; s < sources.size(); s++) { + auto* dgrid = new EulerGrid(0, "Kalium", diff_coef, decay_coef, res); + dgrid->Initialize(); + dgrid->SetUpperThreshold(1e15); + dgrid->ChangeConcentrationBy(sources[s], init); + rm->AddContinuum(dgrid); + + int tot = 20; + for (int t = 0; t < tot; t++) { + dgrid->Diffuse(simulation_time_step); + if (s == 2 || s == 3) { + // Source should pass through periodic boundaries. + EXPECT_FLOAT_EQ(dgrid->GetValue(inside[s]), dgrid->GetValue(outside[s])); + } + else { + // Source should not pass through dirichlet and neumann boundaries. + EXPECT_GT(dgrid->GetValue(inside[s]), dgrid->GetValue(outside[s])); + } + } + rm->RemoveContinuum(0); + } +} + +TEST(DiffusionTest, EulerComplexBoundariesCase3) { + double simulation_time_step{0.1}; + auto set_param = [&](Param* param) { + param->bound_space = Param::BoundSpaceMode::kClosed; + param->min_bound = -50; + param->max_bound = 50; + param->diffusion_boundary_condition = "Complex"; + param->diffusion_bc_x_min = "Periodic"; + param->diffusion_bc_x_max = "Periodic"; + param->diffusion_bc_y_min = "Neumann"; + param->diffusion_bc_y_max = "Dirichlet"; + param->diffusion_bc_z_min = "Periodic"; + param->diffusion_bc_z_max = "Periodic"; + }; + Simulation simulation(TEST_NAME, set_param); + simulation.GetEnvironment()->Update(); + auto* rm = simulation.GetResourceManager(); + + double decay_coef = 0.0; + double diff_coef = 400.0; + int res = 5; + + // Constant Incoming flux at neumann boundary on one side of simulation and zero conc dirichlet + // boundary on the other side. + auto* dgrid = new EulerGrid(0, "Kalium", diff_coef, decay_coef, res); + dgrid->Initialize(); + dgrid->SetBoundaryCondition_neumann(std::make_unique(-1.0)); + dgrid->SetBoundaryCondition_dirichlet(std::make_unique(0.0)); + dgrid->SetUpperThreshold(1e15); + rm->AddContinuum(dgrid); + + double total_conc = 0.0; + double total_conc_backup = 0.0; + double close = 0.0; + double far = 0.0; + + int tot = 2500; + for (int t=0; tDiffuse(simulation_time_step); + auto conc = dgrid->GetAllConcentrations(); + total_conc = 0.0; + for (size_t i = 0; i < dgrid->GetNumBoxes(); i++) { + total_conc += conc[i]; + } + // We should see a higher concentration at neumann boundary compared to dirichlet + // boundary at any point in the simulation. + close = std::abs(dgrid->GetValue({0, -40, 0})); + far = std::abs(dgrid->GetValue({0, 40, 0})); + EXPECT_LT(far, close); + // The total concentration in the system increases at the start + if (t < 100) { + EXPECT_GT(total_conc, total_conc_backup); + } + total_conc_backup = total_conc; + } + // As the concenentration at the dirichlet boundary increases, so does the flux out of the + // boundary. Eventually this flux will equal the flux concentration coming in from the neumann + // boundary. The total concentation reaches a steady state. + EXPECT_LT((total_conc - total_conc_backup) / total_conc, 0.0001); + + rm->RemoveContinuum(0); +} + TEST(DiffusionTest, EulerDirichletBoundaries) { double simulation_time_step{0.1}; auto set_param = [](auto* param) { diff --git a/util/version/generate_version_files.py b/util/version/generate_version_files.py index 5fd71373e..b66b67330 100755 --- a/util/version/generate_version_files.py +++ b/util/version/generate_version_files.py @@ -84,13 +84,13 @@ def UpdateVersionInfo(version_file, version_string): version = GetGitDescribeString(sys.argv[4], sys.argv[3]) # extract information - search = re.search('v([0-9]+)\.([0-9]+)\.([0-9]+)\-(.*)', version) + search = re.search(r'v([0-9]+)\.([0-9]+)\.([0-9]+)\-(.*)', version) if search != None: major = search.group(1) minor = search.group(2) patch = search.group(3) else: - search = re.search('v([0-9]+)\.([0-9]+)', version) + search = re.search(r'v([0-9]+)\.([0-9]+)', version) major = search.group(1) minor = search.group(2) patch = 0