Migration to github: Added tutorials.

Signed-off-by: Franz R. Sattler <[email protected]>

Tutorials/.gitignore

@@ -0,0 +1,4 @@
+*/build/
+*/TraceBuffer/
+*/diagrams/
+*/flows/

Tutorials/tut1/CMakeLists.txt

@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 3.26.4)
+
+project(tut1)
+
+find_package(DiFfRG REQUIRED HINTS /opt/DiFfRG)
+
+add_executable(tut1 tut1.cc)
+setup_application(tut1)
+
+execute_process(
+  COMMAND
+    ${CMAKE_COMMAND} "-E" "create_symlink"
+    "${CMAKE_CURRENT_SOURCE_DIR}/parameter.json"
+    "${CMAKE_CURRENT_BINARY_DIR}/parameter.json")

Tutorials/tut1/build.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+# ##############################################################################
+# Script setup
+# ##############################################################################
+
+threads='1'
+while getopts j: flag; do
+    case "${flag}" in
+    j) threads=${OPTARG} ;;
+    esac
+done
+scriptpath="$(
+    cd -- "$(dirname "$0")" >/dev/null 2>&1
+    pwd -P
+)"
+
+if [[ "$OSTYPE" =~ ^darwin ]]; then
+    export OpenMP_ROOT=$(brew --prefix)/opt/libomp
+fi
+
+# ##############################################################################
+# Build application
+# ##############################################################################
+
+mkdir -p build
+cd build
+cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS='-O3 -march=native -ffast-math -fno-finite-math-only' ..
+make -j"$threads"

Tutorials/tut1/model.hh

@@ -0,0 +1,58 @@
+#pragma once
+
+#include <DiFfRG/model/model.hh>
+
+using namespace DiFfRG;
+
+struct Parameters {
+  Parameters(const JSONValue &value)
+  {
+    try {
+      a = value.get_double("/physical/a");
+      b = value.get_double("/physical/b");
+      c = value.get_double("/physical/c");
+      d = value.get_double("/physical/d");
+    } catch (std::exception &e) {
+      std::cout << "Error in reading parameters: " << e.what() << std::endl;
+    }
+  }
+  double a, b, c, d;
+};
+
+using FEFunctionDesc = FEFunctionDescriptor<Scalar<"u">>;
+using Components = ComponentDescriptor<FEFunctionDesc>;
+constexpr auto idxf = FEFunctionDesc{};
+
+/**
+ * @brief This class implements the numerical model for Burgers' equation.
+ */
+class Tut1 : public def::AbstractModel<Tut1, Components>,
+             public def::fRG,                  // this handles the fRG time
+             public def::FlowBoundaries<Tut1>, // use Inflow/Outflow boundaries
+             public def::AD<Tut1>              // define all jacobians per AD
+{
+private:
+  const Parameters prm;
+
+public:
+  static constexpr uint dim = 1;
+
+  Tut1(const JSONValue &json) : def::fRG(json.get_double("/physical/Lambda")), prm(json) {}
+
+  template <typename Vector> void initial_condition(const Point<dim> &pos, Vector &values) const
+  {
+    const auto x = pos[0];
+    values[idxf("u")] = prm.a + prm.b * powr<1>(x) + prm.c * powr<2>(x) + prm.d * powr<3>(x);
+  }
+
+  template <typename NT, typename Solution>
+  void flux(std::array<Tensor<1, dim, NT>, Components::count_fe_functions(0)> &flux, const Point<dim> & /*pos*/,
+            const Solution &sol) const
+  {
+    const auto &fe_functions = get<"fe_functions">(sol);
+
+    const auto u = fe_functions[idxf("u")];
+
+    flux[idxf("u")][0] = 0.5 * powr<2>(u);
+  }
+};

Tutorials/tut1/parameter.json

@@ -0,0 +1,56 @@
+{
+  "physical": {
+    "Lambda" : 1.0,
+    "a" : 0.0,
+    "b" : 1.0,
+    "c" : 0.0,
+    "d" : 0.0
+  },
+  "discretization": {
+    "fe_order": 3,
+    "threads": 8,
+    "batch_size": 64,
+    "overintegration": 0,
+    "output_subdivisions": 2,
+
+    "EoM_abs_tol": 1e-10,
+    "EoM_max_iter": 200,
+
+    "grid": {
+      "x_grid": "0:1e-2:1",
+      "y_grid": "0:0.1:1",
+      "z_grid": "0:0.1:1",
+      "refine": 0
+    },
+    "adaptivity": {
+      "start_adapt_at": 0E0,
+      "adapt_dt": 1E-1,
+      "level": 0,
+      "refine_percent": 1E-1,
+      "coarsen_percent": 5E-2
+    }
+  },
+  "timestepping": {
+    "final_time": 10.0,
+    "output_dt": 1E-1,
+    "explicit": {
+      "dt": 1E-4,
+      "minimal_dt": 1E-6,
+      "maximal_dt": 1E-1,
+      "abs_tol": 1E-12,
+      "rel_tol": 1E-4
+    },
+    "implicit": {
+      "dt": 1E-4,
+      "minimal_dt": 1E-6,
+      "maximal_dt": 1E-1,
+      "abs_tol": 1E-13,
+      "rel_tol": 1E-7
+    }
+  },
+  "output": {
+    "verbosity": 1,
+    "folder": "./",
+    "name": "output"
+  }
+}

Tutorials/tut1/tut1.cc

@@ -0,0 +1,50 @@
+#include <DiFfRG/common/configuration_helper.hh>
+#include <DiFfRG/common/utils.hh>
+#include <DiFfRG/discretization/discretization.hh>
+#include <DiFfRG/timestepping/timestepping.hh>
+
+#include "model.hh"
+
+using namespace DiFfRG;
+
+int main(int argc, char *argv[])
+{
+  // get all needed parameters and parse from the CLI
+  ConfigurationHelper config_helper(argc, argv);
+  const auto json = config_helper.get_json();
+
+  // Choices for types
+  using Model = Tut1;
+  constexpr uint dim = Model::dim;
+  using Discretization = CG::Discretization<Model::Components, double, RectangularMesh<dim>>;
+  using VectorType = typename Discretization::VectorType;
+  using SparseMatrixType = typename Discretization::SparseMatrixType;
+  using Assembler = CG::Assembler<Discretization, Model>;
+  using TimeStepper = TimeStepperSUNDIALS_IDA<VectorType, SparseMatrixType, dim, UMFPack>;
+
+  // Define the objects needed to run the simulation
+  Model model(json);
+  RectangularMesh<dim> mesh(json);
+  Discretization discretization(mesh, json);
+  Assembler assembler(discretization, model, json);
+  TimeStepper time_stepper(json, &assembler);
+
+  // Set up the initial condition
+  FlowingVariables initial_condition(discretization);
+  initial_condition.interpolate(model);
+
+  // Now we start the timestepping
+  Timer timer;
+  try {
+    time_stepper.run(&initial_condition, 0., json.get_double("/timestepping/final_time"));
+  } catch (std::exception &e) {
+    spdlog::get("log")->error("Simulation finished with exception {}", e.what());
+    return -1;
+  }
+
+  // We print a bit of exit information.
+  assembler.log("log");
+  auto time = timer.wall_time();
+  spdlog::get("log")->info("Simulation finished after " + time_format(time));
+  return 0;
+}

Tutorials/tut2/CMakeLists.txt

@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 3.26.4)
+
+project(tut2)
+
+find_package(DiFfRG REQUIRED HINTS /opt/DiFfRG)
+
+add_executable(tut2 tut2.cc)
+setup_application(tut2)
+
+execute_process(
+  COMMAND
+    ${CMAKE_COMMAND} "-E" "create_symlink"
+    "${CMAKE_CURRENT_SOURCE_DIR}/parameter.json"
+    "${CMAKE_CURRENT_BINARY_DIR}/parameter.json")

Tutorials/tut2/build.sh

@@ -0,0 +1,29 @@
+#!/bin/bash
+
+# ##############################################################################
+# Script setup
+# ##############################################################################
+
+threads='1'
+while getopts j: flag; do
+    case "${flag}" in
+    j) threads=${OPTARG} ;;
+    esac
+done
+scriptpath="$(
+    cd -- "$(dirname "$0")" >/dev/null 2>&1
+    pwd -P
+)"
+
+if [[ "$OSTYPE" =~ ^darwin ]]; then
+    export OpenMP_ROOT=$(brew --prefix)/opt/libomp
+fi
+
+# ##############################################################################
+# Build application
+# ##############################################################################
+
+mkdir -p build
+cd build
+cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS='-O3 -march=native -ffast-math -fno-finite-math-only' ..
+make -j"$threads"

Tutorials/tut2/model.hh

@@ -0,0 +1,84 @@
+#pragma once
+
+#include <DiFfRG/model/model.hh>
+#include <DiFfRG/physics/physics.hh>
+
+using namespace DiFfRG;
+
+struct Parameters {
+  Parameters(const JSONValue &value)
+  {
+    try {
+      Nf = value.get_double("/physical/Nf");
+      Nc = value.get_double("/physical/Nc");
+
+      T = value.get_double("/physical/T");
+      muq = value.get_double("/physical/muq");
+
+      m2Phi = value.get_double("/physical/m2Phi");
+      lambdaPhi = value.get_double("/physical/lambdaPhi");
+      hPhi = value.get_double("/physical/hPhi");
+    } catch (std::exception &e) {
+      std::cout << "Error in reading parameters: " << e.what() << std::endl;
+    }
+  }
+  double Nf, Nc, T, muq, m2Phi, lambdaPhi, hPhi;
+};
+
+using FEFunctionDesc = FEFunctionDescriptor<Scalar<"u">>;
+using Components = ComponentDescriptor<FEFunctionDesc>;
+constexpr auto idxf = FEFunctionDesc{};
+
+namespace tF = fRG::TFLitimSpatial; // use Litim threshold functions
+
+/**
+ * @brief This class implements the numerical model for a Quark-Meson model.
+ */
+class Tut2 : public def::AbstractModel<Tut2, Components>,
+             public def::fRG,                  // this handles the fRG time
+             public def::FlowBoundaries<Tut2>, // use Inflow/Outflow boundaries
+             public def::AD<Tut2>              // define all jacobians per AD
+{
+private:
+  const Parameters prm;
+
+public:
+  static constexpr uint dim = 1;
+
+  Tut2(const JSONValue &json) : def::fRG(json.get_double("/physical/Lambda")), prm(json) {}
+
+  template <typename Vector> void initial_condition(const Point<dim> &pos, Vector &values) const
+  {
+    const auto rhoPhi = pos[0];
+    values[idxf("u")] = prm.m2Phi + prm.lambdaPhi / 2. * rhoPhi;
+  }
+
+  template <typename NT, typename Solution>
+  void flux(std::array<Tensor<1, dim, NT>, Components::count_fe_functions(0)> &flux, const Point<dim> &pos,
+            const Solution &sol) const
+  {
+    const auto rhoPhi = pos[0];
+    const auto &fe_functions = get<"fe_functions">(sol);
+    const auto &fe_derivatives = get<"fe_derivatives">(sol);
+
+    const auto m2Quark = powr<2>(prm.hPhi) * rhoPhi / (prm.Nf);
+    const auto m2Pion = fe_functions[idxf("u")];
+    const auto m2Sigma = m2Pion + 2. * rhoPhi * fe_derivatives[idxf("u")][0];
+
+    flux[idxf("u")][0] = fluxPion(m2Pion) + fluxSigma(m2Sigma) + fluxQuark(m2Quark);
+  }
+
+private:
+  template <typename NT> NT fluxSigma(const NT m2Sigma) const
+  {
+    return (k4 * tF::lB<0, NT>(m2Sigma / k2, k, prm.T, 4)) / (4. * powr<2>(M_PI));
+  }
+  template <typename NT> NT fluxPion(const NT m2Pion) const
+  {
+    return (k4 * (powr<2>(prm.Nf) - 1.) * tF::lB<0, NT>(m2Pion / k2, k, prm.T, 4)) / (4. * powr<2>(M_PI));
+  }
+  template <typename NT> NT fluxQuark(const NT m2Quark) const
+  {
+    return -((k4 * prm.Nc * prm.Nf * tF::lF<0, NT>(m2Quark / k2, k, prm.T, prm.muq, 4)) / powr<2>(M_PI));
+  }
+};