model_initializer.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
// BioDynaMo collaboration. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
//
// See the LICENSE file distributed with this work for details.
// See the NOTICE file distributed with this work for additional information
// regarding copyright ownership.
//
// -----------------------------------------------------------------------------
 
#ifndef CORE_MODEL_INITIALIZER_H_
#define CORE_MODEL_INITIALIZER_H_
 
#include <Math/DistFunc.h>
#include <omp.h>
#include <ctime>
#include <string>
#include <vector>
 
#include "core/container/math_array.h"
#include "core/diffusion/diffusion_grid.h"
#include "core/resource_manager.h"
#include "core/simulation.h"
#include "core/util/random.h"
 
class EulerGrid;
class EulerDepletionGrid;
 
namespace bdm {
 
struct ModelInitializer {
  template <typename Function>
  static void Grid3D(size_t agents_per_dim, real_t space,
                     Function agent_builder) {
#pragma omp parallel
    {
      auto* sim = Simulation::GetActive();
      auto* ctxt = sim->GetExecutionContext();
 
#pragma omp for
      for (size_t x = 0; x < agents_per_dim; x++) {
        auto x_pos = x * space;
        for (size_t y = 0; y < agents_per_dim; y++) {
          auto y_pos = y * space;
          for (size_t z = 0; z < agents_per_dim; z++) {
            auto* new_agent = agent_builder({x_pos, y_pos, z * space});
            ctxt->AddAgent(new_agent);
          }
        }
      }
    }
  }
 
  template <typename Function>
  static void Grid3D(const std::array<size_t, 3>& agents_per_dim, real_t space,
                     Function agent_builder) {
#pragma omp parallel
    {
      auto* sim = Simulation::GetActive();
      auto* ctxt = sim->GetExecutionContext();
 
#pragma omp for
      for (size_t x = 0; x < agents_per_dim[0]; x++) {
        auto x_pos = x * space;
        for (size_t y = 0; y < agents_per_dim[1]; y++) {
          auto y_pos = y * space;
          for (size_t z = 0; z < agents_per_dim[2]; z++) {
            auto* new_agent = agent_builder({x_pos, y_pos, z * space});
            ctxt->AddAgent(new_agent);
          }
        }
      }
    }
  }
 
  template <typename Function>
  static void CreateAgents(const std::vector<Real3>& positions,
                           Function agent_builder) {
#pragma omp parallel
    {
      auto* sim = Simulation::GetActive();
      auto* ctxt = sim->GetExecutionContext();
 
#pragma omp for
      for (size_t i = 0; i < positions.size(); i++) {
        auto* new_agent =
            agent_builder({positions[i][0], positions[i][1], positions[i][2]});
        ctxt->AddAgent(new_agent);
      }
    }
  }
 
  template <typename Function>
  static void CreateAgentsRandom(real_t min, real_t max, uint64_t num_agents,
                                 Function agent_builder,
                                 DistributionRng<real_t>* rng = nullptr) {
#pragma omp parallel
    {
      auto* sim = Simulation::GetActive();
      auto* ctxt = sim->GetExecutionContext();
      auto* random = sim->GetRandom();
 
#pragma omp for
      for (uint64_t i = 0; i < num_agents; i++) {
        if (rng != nullptr) {
          Real3 pos;
          bool in_range = false;
          do {
            pos = rng->Sample3();
            in_range = (pos[0] >= min) && (pos[0] <= max) && (pos[1] >= min) &&
                       (pos[1] <= max) && (pos[2] >= min) && (pos[2] <= max);
          } while (!in_range);
          auto* new_agent = agent_builder(pos);
          ctxt->AddAgent(new_agent);
        } else {
          auto* new_agent = agent_builder(random->UniformArray<3>(min, max));
          ctxt->AddAgent(new_agent);
        }
      }
    }
  }
 
  template <typename Function>
  static void CreateAgentsOnSurface(
      real_t (*f)(const real_t*, const real_t*),
      const FixedSizeVector<real_t, 10>& fn_params, real_t xmin, real_t xmax,
      real_t deltax, real_t ymin, real_t ymax, real_t deltay,
      Function agent_builder) {
#pragma omp parallel
    {
      auto* sim = Simulation::GetActive();
      auto* ctxt = sim->GetExecutionContext();
 
      auto xiterations =
          static_cast<uint64_t>(std::floor((xmax - xmin) / deltax));
      auto yiterations =
          static_cast<uint64_t>(std::floor((ymax - ymin) / deltay));
 
#pragma omp for
      for (uint64_t xit = 0; xit < xiterations; ++xit) {
        real_t x = xmin + xit * deltax;
        for (uint64_t yit = 0; yit < yiterations; ++yit) {
          real_t y = ymin + yit * deltay;
          Real3 pos = {x, y};
          pos[2] = f(pos.data(), fn_params.data());
          ctxt->AddAgent(agent_builder(pos));
        }
      }
    }
  }
 
  template <typename Function>
  static void CreateAgentsOnSurfaceRndm(
      real_t (*f)(const real_t*, const real_t*),
      const FixedSizeVector<real_t, 10>& fn_params, real_t xmin, real_t xmax,
      real_t ymin, real_t ymax, uint64_t num_agents, Function agent_builder) {
#pragma omp parallel
    {
      auto* sim = Simulation::GetActive();
      auto* ctxt = sim->GetExecutionContext();
      auto* random = sim->GetRandom();
 
#pragma omp for
      for (uint64_t i = 0; i < num_agents; ++i) {
        Real3 pos = {random->Uniform(xmin, xmax), random->Uniform(ymin, ymax)};
        pos[2] = f(pos.data(), fn_params.data());
        ctxt->AddAgent(agent_builder(pos));
      }
    }
  }
 
  template <typename Function>
  static void CreateAgentsOnSphereRndm(const Real3& center, real_t radius,
                                       uint64_t num_agents,
                                       Function agent_builder) {
#pragma omp parallel
    {
      auto* sim = Simulation::GetActive();
      auto* ctxt = sim->GetExecutionContext();
      auto* random = sim->GetRandom();
 
#pragma omp for
      for (uint64_t i = 0; i < num_agents; i++) {
        auto pos = random->Sphere(radius) + center;
        auto* new_agent = agent_builder(pos);
        ctxt->AddAgent(new_agent);
      }
    }
  }
 
  template <typename Function>
  static void CreateAgentsInSphereRndm(const Real3& center, real_t radius,
                                       uint64_t num_agents,
                                       Function agent_builder) {
    // We use a probability density function (PDF) to model the probability of
    // an agent to occur at a distance `r>=0` of the center. As the surface of
    // a sphere scales as `r^2`, the PDF does as well. Thus
    // `p(r)=a*r^2*\Theta(R-r)`, where `\Theta` is a heavyside function and R is
    // largest allowed radius (interpretation: no agents outside the sphere). We
    // can fix `a` by requiring `\int_0^\inf p(r') dr' = 1` and obtain
    // `a=3/R^3`.
    auto radial_pdf_sphere = [](const double* x, const double* params) {
      double R{params[0]};
      double r{x[0]};
      if (r > 0.0 && r <= R) {
        return 3.0 * std::pow(r, 2.0) / std::pow(R, 3.0);
      } else {
        return 0.0;
      }
    };
 
    // Get a random number generator to sample from our PDF.
    auto* random = Simulation::GetActive()->GetRandom();
    auto rng =
        random->GetUserDefinedDistRng1D(radial_pdf_sphere, {radius}, 0, radius);
 
    // Create a random radius for each of the agents. Note: this is done
    // serially because we GetUserDefinedDistRng1D does not work in parallel
    // regions at the moment.
    std::vector<real_t> random_radius;
    random_radius.resize(num_agents);
    for (size_t i = 0; i < num_agents; i++) {
      random_radius[i] = rng.Sample();
    }
#pragma omp parallel shared(random_radius)
    {
      auto* ctxt_tl = Simulation::GetActive()->GetExecutionContext();
#pragma omp for schedule(static)
      for (uint64_t i = 0; i < num_agents; i++) {
        auto pos = random->Sphere(random_radius[i]) + center;
        auto* new_agent = agent_builder(pos);
        ctxt_tl->AddAgent(new_agent);
      }
    }
  }
 
  static void DefineSubstance(
      int substance_id, const std::string& substance_name,
      real_t diffusion_coeff, real_t decay_constant, int resolution = 10,
      const std::vector<real_t>& binding_coefficients = {},
      const std::vector<int>& binding_substances = {});
 
  template <typename F>
  static void InitializeSubstance(int substance_id, F function) {
    auto* sim = Simulation::GetActive();
    auto* rm = sim->GetResourceManager();
    auto diffusion_grid = rm->GetDiffusionGrid(substance_id);
    diffusion_grid->AddInitializer(function);
  }
 
  static void AddBoundaryConditions(int substance_id,
                                    BoundaryConditionType bc_type,
                                    std::unique_ptr<BoundaryCondition> bc) {
    auto* sim = Simulation::GetActive();
    const auto* rm = sim->GetResourceManager();
    auto diffusion_grid = rm->GetDiffusionGrid(substance_id);
    diffusion_grid->SetBoundaryConditionType(bc_type);
    diffusion_grid->SetBoundaryCondition(std::move(bc));
  }
};
 
}  // namespace bdm
 
#endif  // CORE_MODEL_INITIALIZER_H_