astro_demos

demo/mars/CMakeLists.txt

@@ -0,0 +1,35 @@
+# -----------------------------------------------------------------------------
+#
+# 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.
+#
+# -----------------------------------------------------------------------------
+cmake_minimum_required(VERSION 3.19.3)
+project(mars)
+
+# BioDynaMo curretly uses the C++17 standard.
+set(CMAKE_CXX_STANDARD 17)
+
+# Use BioDynaMo in this project.
+find_package(BioDynaMo REQUIRED)
+
+# See UseBioDynaMo.cmake in your BioDynaMo build folder for details.
+# Note that BioDynaMo provides gtest header/libraries in its include/lib dir.
+include(${BDM_USE_FILE})
+
+# Consider all files in src/ for BioDynaMo simulation.
+include_directories("src")
+file(GLOB_RECURSE PROJECT_HEADERS src/*.h)
+file(GLOB_RECURSE PROJECT_SOURCES src/*.cc)
+
+bdm_add_executable(${CMAKE_PROJECT_NAME}
+                   HEADERS ${PROJECT_HEADERS}
+                   SOURCES ${PROJECT_SOURCES}
+                   LIBRARIES ${BDM_REQUIRED_LIBRARIES})

demo/mars/README.md

@@ -0,0 +1,73 @@
+# Mars
+
+This simulation models the gravitational interactions between Mars and its satellites, Phobos and Deimos, using realistic dimensions and masses. This demo consists of the planet Mars and its two satellites as agents, with gravitational interactions modeled using the gravity behavior.
+
+<div align="center">
+  <img src="Mars.png" alt="Simulation">
+</div>
+
+## Contents
+
+- **Celestial Object Agents**
+- **Gravity Behavior**
+
+## Assumptions
+
+To simplify the simulation, the following assumptions were made:
+ - All celestial objects (Mars, Phobos, and Deimos) are modeled as **spheres**.
+ - Celestial objects are treated as **point masses** for computational efficiency.
+ - Agents can't physically interact with each other.
+
+## Details
+
+- **Mars**: Radius ~3,389.5 km, Mass ~6.417 × 10^23 kg.
+- **Phobos**: Radius ~11.267 km, Mass ~1.0659 × 10^16 kg, Average orbital distance ~9,376 km.
+- **Deimos**: Radius ~6.2 km, Mass ~1.4762 × 10^15 kg, Average orbital distance ~23,463 km.
+
+## Newton's Law of Universal Gravitation
+
+Newton’s law of universal gravitation governs the movement of all celestial objects in the simulation, describing the force of attraction between two masses.
+
+The force of attraction between two objects is proportional to both their masses and inversely proportional to the square of the distance between them.
+
+$$F = \frac{G * (m_1 * m_2)}{r^2}$$
+
+where:
+-    $F$    is the gravitational force between the two masses,
+-    $G$    is the gravitational constant (6.674 × 10^-11 N·m²/kg²),
+-    $m_1$   and $m_2$ are the masses of the two objects, and
+-    $r$    is the distance between the centers of the two masses.
+
+By knowing the force of attraction and therefore the acceleration of a celestial object, the displacement can be predicted by solving the below differential equation:
+
+$$\frac{d^2x}{dt^2} = a$$
+
+---
+
+## Run the demos
+
+In order to run the demos, BioDynamo has to be correctly installed and sourced.
+
+```bash
+cd [path_to_biodynamo]demos/NewtonsLawTest
+bdm run
+bdm view
+```
+In order to clearly observe the satellites in relation to the planet, feel free to scale up the diameter or volume of the satellites in the simulation.
+
+## Verify the Results
+
+By running the following commands:
+
+```bash
+pip install vtk 
+```
+
+```bash
+cd [path_to_biodynamo]/demos/SolarSystem
+python3 check_results.py
+```
+
+You can verify the accuracy of the simulation by comparing the time it takes for a celestial object to complete a full rotation.
+
+The data were acquired from the [NASA Mars Moons Facts](https://science.nasa.gov/mars/moons/facts/), which provides real-world orbital and rotational period measurements.

demo/mars/bdm.toml

@@ -0,0 +1,11 @@
+[visualization]
+export = true
+interval = 1
+
+[simulation]
+time_step = 360
+
+[[visualize_agent]]
+name = "Planet"
+[[visualize_agent]]
+name = "Satellite"

demo/mars/check_results.py

@@ -0,0 +1,58 @@
+import vtk
+import os
+
+steps = 500
+ptuFiles = 8
+reader = vtk.vtkXMLUnstructuredGridReader()
+
+# first read all locations
+
+satellites = ["Phobos", "Deimos"]
+locations = [[], []]
+try:
+    # Read pvtu file for each time step and extract files
+    for i in range(steps):
+        k = 0
+        for j in range(ptuFiles):
+
+            # Find the output files
+            path = f"output/mars/Satellite-{i}_{j}.vtu"
+            if not os.path.exists(path):
+                raise Exception("Output file not found")
+            reader.SetFileName(path)
+            reader.Update()
+            output = reader.GetOutput()
+
+            # Get all points from each file
+            for l in range(output.GetNumberOfPoints()):
+                locations[k].append(output.GetPoints().GetPoint(l))
+                k = k + 1
+
+    # Scientific Data
+    # Since each step is 0.1 days
+    data = [76.6, 303.5]
+    stars = "**********************"
+    print(f'\n\n{stars}{stars}{stars}\n')
+    for satellite in satellites:
+        index = satellites.index(satellite)
+        # Find location after 1 full rotation
+        for step in range(1,steps):
+
+            A = locations[index][step][1] > 0
+            B = locations[index][step-1][1] < 0
+
+            if(A and B):
+                # Find the relative error
+                # Relative error = |x - x_0| / x
+                # where x is actual value
+                # and x is the measured value
+                relativeError = abs(step - data[index]) / data[index]
+                print(f"    The relative error of {satellite}'s full rotation is {relativeError*100:.2f}%")
+                break
+    print(f'\n{stars}{stars}{stars}') 
+    print(f"\n\nPrevious runs of this algorithm returned the following output:\n")
+    print(f"    The relative error of Phobo's full rotation is 0.52%")
+    print(f"    The relative error of Deimo's full rotation is 0.16%")
+
+except Exception as e:
+    print(f"error: {e}") 

demo/mars/src/CelestialObjects.h

@@ -0,0 +1,94 @@
+// -----------------------------------------------------------------------------
+//
+// 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 CELESTIALOBJECTS_H_
+#define CELESTIALOBJECTS_H_
+
+#include "biodynamo.h"
+
+namespace bdm {
+namespace astrophysics {
+
+// creating the custom agent that
+// defines celestial objects
+class CelestialObject : public Cell {
+  BDM_AGENT_HEADER(CelestialObject, Cell, 1);
+
+ public:
+  // constructors
+  CelestialObject() { UpdateVolume(); }
+
+  explicit CelestialObject(real_t diameter){
+    this->SetDiameter(diameter);
+    UpdateVolume();
+  }
+
+  explicit CelestialObject(const Real3& position){
+    this->SetPosition(position);
+    UpdateVolume();
+  }
+
+  // default destructor
+  ~CelestialObject() override = default;
+
+  Real3 GetVelocity() const { return velocity_; }
+
+  void SetVelocity(const Real3& vel) { velocity_ = vel; }
+
+ private:
+  Real3 velocity_ = {0, 0, 0};
+};
+
+// Planet subclass
+class Planet : public CelestialObject {
+  BDM_AGENT_HEADER(Planet, CelestialObject, 1);
+
+ public:
+  Planet() : CelestialObject() {}
+
+  explicit Planet(real_t diameter) : CelestialObject(diameter) {}
+
+  explicit Planet(const Real3 position) : CelestialObject(position) {}
+
+  ~Planet() override = default;
+};
+
+// Satellite subclass
+class Satellite : public CelestialObject {
+  BDM_AGENT_HEADER(Satellite, CelestialObject, 1);
+
+ public:
+  Satellite() : CelestialObject() {}
+
+  explicit Satellite(real_t diameter) : CelestialObject(diameter) {}
+
+  explicit Satellite(const Real3 position) : CelestialObject(position) {}
+
+  explicit Satellite(const Planet* main_planet, real_t diameter,
+                     real_t distance, real_t initial_speed) {
+    Real3 pos = main_planet->GetPosition();
+    pos[0] += distance;
+    Real3 initial_velocity = {0, initial_speed, 0};
+    this->SetVelocity(main_planet->GetVelocity() + initial_velocity);
+
+    this->SetPosition(pos);
+    this->SetDiameter(diameter);
+  }
+
+  ~Satellite() override = default;
+};
+
+}  // namespace astrophysics
+}  // namespace bdm
+
+#endif  // CELESTIALOBJECT_H_

demo/mars/src/Gravity.h

@@ -0,0 +1,116 @@
+// -----------------------------------------------------------------------------
+//
+// 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 GRAVITY_H_
+#define GRAVITY_H_
+
+#include "CelestialObjects.h"
+#include "biodynamo.h"
+
+namespace bdm {
+namespace astrophysics {
+
+// Gravitational constant in m^3/kg/s^2
+inline real_t gG = 6.674e-11;
+
+// set units of G constant
+inline void SetUnits(char units = 'O', int times = 1) {
+  gG = 6.674e-11;
+  switch (units) {
+    case 'k':
+      gG *= pow(times * 1e-3, 3);
+      break;
+    case 'M':
+      gG *= pow(times * 1e-6, 3);
+      break;
+    case 'G':
+      gG *= pow(times * 1e-9, 3);
+      break;
+    default:
+      gG *= pow(times, 3);
+      break;
+  }
+}
+
+// a behavior that applies Newtons Law for gravity and calculates the
+// displacement for any celestial object
+class Gravity : public Behavior {
+  Real3 acceleration_ = {0, 0, 0};
+  real_t dt_ = Simulation::GetActive()->GetParam()->simulation_time_step;
+
+ public:
+  Behavior* New() const override { return new Gravity(); }
+
+  Behavior* NewCopy() const override { return new Gravity(*this); }
+
+  void Run(Agent* agent) override {
+    CelestialObject* celestial = dynamic_cast<CelestialObject*>(agent);
+
+    Real3 celestial_pos = celestial->GetPosition();
+
+    auto* functor = new GravityFunctor(celestial, gG);
+    // for each agent calculate acceleration
+    Simulation::GetActive()->GetResourceManager()->ForEachAgentParallel(
+        *functor);
+
+    dt_ = Simulation::GetActive()->GetParam()->simulation_time_step;
+    Real3 velocity = celestial->GetVelocity();
+    Real3 acceleration = functor->acceleration;
+    Real3 position;
+
+    // explicit Euler's method
+    velocity += acceleration * dt_;
+    position = celestial_pos + velocity * dt_;
+
+    delete functor;
+
+    celestial->SetVelocity(velocity);
+    celestial->SetPosition(position);
+  }
+
+  class GravityFunctor : public Functor<void, Agent*, AgentHandle> {
+   public:
+    CelestialObject* celestial;
+    real_t g;
+    Real3 acceleration;
+
+    GravityFunctor(CelestialObject* celestial, real_t g)
+        : celestial(celestial), g(g), acceleration({0, 0, 0}) {}
+
+    // override the function that overloads () operator to
+    // calculate acceleration based on given agent
+    void operator()(Agent* other_agent, AgentHandle) override {
+      CelestialObject* other_celestial =
+          dynamic_cast<CelestialObject*>(other_agent);
+
+      if (!other_celestial || other_celestial == celestial) {
+        return;
+      }
+
+      Real3 dir = other_celestial->GetPosition() - celestial->GetPosition();
+      real_t magnitude = dir.Norm();
+
+      dir.Normalize();
+      real_t m1 = celestial->GetMass();
+      real_t m2 = other_celestial->GetMass();
+      real_t f = g * m1 * m2 / (magnitude * magnitude);
+
+      acceleration += (f / m1) * dir;
+    };
+  };
+};
+}  // namespace astrophysics
+
+}  // namespace bdm
+
+#endif  // GRAVITY_H_