diff --git a/CMakeLists.txt b/CMakeLists.txt
index 61e560c..66b8001 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -105,6 +105,7 @@ add_subdirectory(src/core)
if (BUILD_SAMPLES)
add_subdirectory(samples/standalone)
+ add_subdirectory(samples/standalone_mpi)
add_subdirectory(samples/ctest)
add_subdirectory(samples/cpptest)
add_subdirectory(samples/mpitest)
diff --git a/samples/standalone_mpi/CMakeLists.txt b/samples/standalone_mpi/CMakeLists.txt
new file mode 100644
index 0000000..979b4d4
--- /dev/null
+++ b/samples/standalone_mpi/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_executable(ac_run_mpi main.cc host_memory.cc host_forcing.cc config_loader.cc)
+target_link_libraries(ac_run_mpi astaroth_utils astaroth_core)
diff --git a/samples/standalone_mpi/config_loader.cc b/samples/standalone_mpi/config_loader.cc
new file mode 100644
index 0000000..6577707
--- /dev/null
+++ b/samples/standalone_mpi/config_loader.cc
@@ -0,0 +1,184 @@
+/*
+ Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+
+/**
+ * @file
+ * \brief Brief info.
+ *
+ * Detailed info.
+ *
+ */
+#include "config_loader.h"
+
+#include // UINT_MAX
+#include // uint8_t, uint32_t
+#include // print
+#include // memset
+
+#include "errchk.h"
+#include "math_utils.h"
+
+/**
+ \brief Find the index of the keyword in names
+ \return Index in range 0...n if the keyword is in names. -1 if the keyword was
+ not found.
+ */
+static int
+find_str(const char keyword[], const char* names[], const int& n)
+{
+ for (int i = 0; i < n; ++i)
+ if (!strcmp(keyword, names[i]))
+ return i;
+
+ return -1;
+}
+
+static void
+parse_config(const char* path, AcMeshInfo* config)
+{
+ FILE* fp;
+ fp = fopen(path, "r");
+ // For knowing which .conf file will be used
+ printf("Config file path: \n %s \n ", path);
+ ERRCHK(fp != NULL);
+
+ const size_t BUF_SIZE = 128;
+ char keyword[BUF_SIZE];
+ char value[BUF_SIZE];
+ int items_matched;
+ while ((items_matched = fscanf(fp, "%s = %s", keyword, value)) != EOF) {
+
+ if (items_matched < 2)
+ continue;
+
+ int idx = -1;
+ if ((idx = find_str(keyword, intparam_names, NUM_INT_PARAMS)) >= 0)
+ config->int_params[idx] = atoi(value);
+ else if ((idx = find_str(keyword, realparam_names, NUM_REAL_PARAMS)) >= 0)
+ config->real_params[idx] = AcReal(atof(value));
+ }
+
+ fclose(fp);
+}
+
+void
+update_config(AcMeshInfo* config)
+{
+ config->int_params[AC_mx] = config->int_params[AC_nx] + STENCIL_ORDER;
+ ///////////// PAD TEST
+ // config->int_params[AC_mx] = config->int_params[AC_nx] + STENCIL_ORDER + PAD_SIZE;
+ ///////////// PAD TEST
+ config->int_params[AC_my] = config->int_params[AC_ny] + STENCIL_ORDER;
+ config->int_params[AC_mz] = config->int_params[AC_nz] + STENCIL_ORDER;
+
+ // Bounds for the computational domain, i.e. nx_min <= i < nx_max
+ config->int_params[AC_nx_min] = STENCIL_ORDER / 2;
+ config->int_params[AC_nx_max] = config->int_params[AC_nx_min] + config->int_params[AC_nx];
+ config->int_params[AC_ny_min] = STENCIL_ORDER / 2;
+ config->int_params[AC_ny_max] = config->int_params[AC_ny] + STENCIL_ORDER / 2;
+ config->int_params[AC_nz_min] = STENCIL_ORDER / 2;
+ config->int_params[AC_nz_max] = config->int_params[AC_nz] + STENCIL_ORDER / 2;
+
+ // Spacing
+ config->real_params[AC_inv_dsx] = AcReal(1.) / config->real_params[AC_dsx];
+ config->real_params[AC_inv_dsy] = AcReal(1.) / config->real_params[AC_dsy];
+ config->real_params[AC_inv_dsz] = AcReal(1.) / config->real_params[AC_dsz];
+ config->real_params[AC_dsmin] = min(
+ config->real_params[AC_dsx], min(config->real_params[AC_dsy], config->real_params[AC_dsz]));
+
+ // Real grid coordanates (DEFINE FOR GRID WITH THE GHOST ZONES)
+ config->real_params[AC_xlen] = config->real_params[AC_dsx] * config->int_params[AC_mx];
+ config->real_params[AC_ylen] = config->real_params[AC_dsy] * config->int_params[AC_my];
+ config->real_params[AC_zlen] = config->real_params[AC_dsz] * config->int_params[AC_mz];
+
+ config->real_params[AC_xorig] = AcReal(.5) * config->real_params[AC_xlen];
+ config->real_params[AC_yorig] = AcReal(.5) * config->real_params[AC_ylen];
+ config->real_params[AC_zorig] = AcReal(.5) * config->real_params[AC_zlen];
+
+ /* Additional helper params */
+ // Int helpers
+ config->int_params[AC_mxy] = config->int_params[AC_mx] * config->int_params[AC_my];
+ config->int_params[AC_nxy] = config->int_params[AC_nx] * config->int_params[AC_ny];
+ config->int_params[AC_nxyz] = config->int_params[AC_nxy] * config->int_params[AC_nz];
+
+ // Real helpers
+ config->real_params[AC_cs2_sound] = config->real_params[AC_cs_sound] *
+ config->real_params[AC_cs_sound];
+
+ config->real_params[AC_cv_sound] = config->real_params[AC_cp_sound] /
+ config->real_params[AC_gamma];
+
+ AcReal G_CONST_CGS = AcReal(
+ 6.674e-8); // cm^3/(g*s^2) GGS definition //TODO define in a separate module
+ AcReal M_sun = AcReal(1.989e33); // g solar mass
+
+ config->real_params[AC_unit_mass] = (config->real_params[AC_unit_length] *
+ config->real_params[AC_unit_length] *
+ config->real_params[AC_unit_length]) *
+ config->real_params[AC_unit_density];
+
+ config->real_params[AC_M_sink] = config->real_params[AC_M_sink_Msun] * M_sun /
+ config->real_params[AC_unit_mass];
+ config->real_params[AC_M_sink_init] = config->real_params[AC_M_sink_Msun] * M_sun /
+ config->real_params[AC_unit_mass];
+
+ config->real_params[AC_G_const] = G_CONST_CGS / ((config->real_params[AC_unit_velocity] *
+ config->real_params[AC_unit_velocity]) /
+ (config->real_params[AC_unit_density] *
+ config->real_params[AC_unit_length] *
+ config->real_params[AC_unit_length]));
+
+ config->real_params[AC_sq2GM_star] = AcReal(sqrt(AcReal(2) * config->real_params[AC_GM_star]));
+
+#if VERBOSE_PRINTING // Defined in astaroth.h
+ printf("###############################################################\n");
+ printf("Config dimensions recalculated:\n");
+ acPrintMeshInfo(*config);
+ printf("###############################################################\n");
+#endif
+}
+
+/**
+\brief Loads data from astaroth.conf into a config struct.
+\return 0 on success, -1 if there are potentially uninitialized values.
+*/
+int
+load_config(const char* config_path, AcMeshInfo* config)
+{
+ int retval = 0;
+ ERRCHK(config_path);
+
+ // memset reads the second parameter as a byte even though it says int in
+ // the function declaration
+ memset(config, (uint8_t)0xFF, sizeof(*config));
+
+ parse_config(config_path, config);
+ update_config(config);
+
+ // sizeof(config) must be a multiple of 4 bytes for this to work
+ ERRCHK(sizeof(*config) % sizeof(uint32_t) == 0);
+ for (size_t i = 0; i < sizeof(*config) / sizeof(uint32_t); ++i) {
+ if (((uint32_t*)config)[i] == (uint32_t)0xFFFFFFFF) {
+ WARNING("Some config values may be uninitialized. "
+ "See that all are defined in astaroth.conf\n");
+ retval = -1;
+ }
+ }
+ return retval;
+}
diff --git a/samples/standalone_mpi/config_loader.h b/samples/standalone_mpi/config_loader.h
new file mode 100644
index 0000000..5dae73f
--- /dev/null
+++ b/samples/standalone_mpi/config_loader.h
@@ -0,0 +1,34 @@
+/*
+ Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+
+/**
+ * @file
+ * \brief Functions for loading and updating AcMeshInfo.
+ *
+ */
+#pragma once
+#include "astaroth.h"
+
+/** Loads data from the config file */
+int load_config(const char* config_path, AcMeshInfo* config);
+
+/** Recalculates the portion of int parameters which get their values from nx,
+ * ny and nz. Must be called after modifying the config struct or otherwise
+ * contents of the struct will be incorrect */
+void update_config(AcMeshInfo* config);
diff --git a/samples/standalone_mpi/host_forcing.cc b/samples/standalone_mpi/host_forcing.cc
new file mode 100644
index 0000000..300f4ad
--- /dev/null
+++ b/samples/standalone_mpi/host_forcing.cc
@@ -0,0 +1,294 @@
+/*
+ Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+
+/**
+ * @file
+ * \brief Brief info.
+ *
+ * Detailed info.
+ *
+ */
+#include "host_forcing.h"
+
+// #include "math_utils.h"
+#include
+using namespace std;
+
+// The is a wrapper for genering random numbers with a chosen system.
+AcReal
+get_random_number_01()
+{
+ // TODO: Implement better randon number generator http://www.cplusplus.com/reference/random/
+ return AcReal(rand()) / AcReal(RAND_MAX);
+}
+
+static AcReal3
+cross(const AcReal3& a, const AcReal3& b)
+{
+ AcReal3 c;
+
+ c.x = a.y * b.z - a.z * b.y;
+ c.y = a.z * b.x - a.x * b.z;
+ c.z = a.x * b.y - a.y * b.x;
+
+ return c;
+}
+
+static AcReal
+dot(const AcReal3& a, const AcReal3& b)
+{
+ return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+static AcReal3
+vec_norm(const AcReal3& a)
+{
+ AcReal3 c;
+ AcReal norm = dot(a, a);
+
+ c.x = a.x / sqrt(norm);
+ c.y = a.y / sqrt(norm);
+ c.z = a.z / sqrt(norm);
+
+ return c;
+}
+
+static AcReal3
+vec_multi_scal(const AcReal scal, const AcReal3& a)
+{
+ AcReal3 c;
+
+ c.x = a.x * scal;
+ c.y = a.y * scal;
+ c.z = a.z * scal;
+
+ return c;
+}
+
+// Generate forcing wave vector k_force
+AcReal3
+helical_forcing_k_generator(const AcReal kmax, const AcReal kmin)
+{
+ AcReal phi, theta, kk; // Spherical direction coordinates
+ AcReal3 k_force; // forcing wave vector
+
+ AcReal delta_k = kmax - kmin;
+
+ // Generate vector in spherical coordinates
+ phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
+ theta = AcReal(M_PI) * get_random_number_01();
+ kk = delta_k * get_random_number_01() + kmin;
+
+ // Cast into Cartesian form
+ k_force = (AcReal3){AcReal(kk * sin(theta) * cos(phi)), //
+ AcReal(kk * sin(theta) * sin(phi)), //
+ AcReal(kk * cos(theta))};
+
+ // printf("k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x, k_force.y, k_force.z);
+
+ // Round the numbers. In that way k(x/y/z) will get complete waves.
+ k_force.x = round(k_force.x);
+ k_force.y = round(k_force.y);
+ k_force.z = round(k_force.z);
+
+ // printf("After rounding --> k_force.x %f, k_force.y %f, k_force.z %f \n", k_force.x,
+ // k_force.y, k_force.z);
+
+ return k_force;
+}
+
+// Generate the unit perpendicular unit vector e required for helical forcing
+// Addapted from Pencil code forcing.f90 hel_vec() subroutine.
+void
+helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force)
+{
+
+ AcReal3 k_cross_e = cross(k_force, *e_force);
+ k_cross_e = vec_norm(k_cross_e);
+ AcReal3 k_cross_k_cross_e = cross(k_force, k_cross_e);
+ k_cross_k_cross_e = vec_norm(k_cross_k_cross_e);
+ AcReal phi = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
+ AcReal3 ee_tmp1 = vec_multi_scal(cos(phi), k_cross_e);
+ AcReal3 ee_tmp2 = vec_multi_scal(sin(phi), k_cross_k_cross_e);
+
+ *e_force = (AcReal3){ee_tmp1.x + ee_tmp2.x, ee_tmp1.y + ee_tmp2.y, ee_tmp1.z + ee_tmp2.z};
+}
+
+// PC Manual Eq. 223
+void
+helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force,
+ const AcReal3 e_force, const AcReal relhel)
+{
+
+ // k dot e
+ AcReal3 kdote;
+ kdote.x = k_force.x * e_force.x;
+ kdote.y = k_force.y * e_force.y;
+ kdote.z = k_force.z * e_force.z;
+
+ // k cross e
+ AcReal3 k_cross_e;
+ k_cross_e.x = k_force.y * e_force.z - k_force.z * e_force.y;
+ k_cross_e.y = k_force.z * e_force.x - k_force.x * e_force.z;
+ k_cross_e.z = k_force.x * e_force.y - k_force.y * e_force.x;
+
+ // k cross k cross e
+ AcReal3 k_cross_k_cross_e;
+ k_cross_k_cross_e.x = k_force.y * k_cross_e.z - k_force.z * k_cross_e.y;
+ k_cross_k_cross_e.y = k_force.z * k_cross_e.x - k_force.x * k_cross_e.z;
+ k_cross_k_cross_e.z = k_force.x * k_cross_e.y - k_force.y * k_cross_e.x;
+
+ // abs(k)
+ AcReal kabs = sqrt(k_force.x * k_force.x + k_force.y * k_force.y + k_force.z * k_force.z);
+
+ AcReal denominator = sqrt(AcReal(1.0) + relhel * relhel) * kabs *
+ sqrt(kabs * kabs -
+ (kdote.x * kdote.x + kdote.y * kdote.y + kdote.z * kdote.z));
+
+ // MV: I suspect there is a typo in the Pencil Code manual!
+ //*ff_hel_re = (AcReal3){-relhel*kabs*k_cross_e.x/denominator,
+ // -relhel*kabs*k_cross_e.y/denominator,
+ // -relhel*kabs*k_cross_e.z/denominator};
+
+ //*ff_hel_im = (AcReal3){k_cross_k_cross_e.x/denominator,
+ // k_cross_k_cross_e.y/denominator,
+ // k_cross_k_cross_e.z/denominator};
+
+ // See PC forcing.f90 forcing_hel_both()
+ *ff_hel_im = (AcReal3){kabs * k_cross_e.x / denominator, kabs * k_cross_e.y / denominator,
+ kabs * k_cross_e.z / denominator};
+
+ *ff_hel_re = (AcReal3){relhel * k_cross_k_cross_e.x / denominator,
+ relhel * k_cross_k_cross_e.y / denominator,
+ relhel * k_cross_k_cross_e.z / denominator};
+}
+
+// Tool for loading forcing vector information into the device memory
+// %JP: Added a generic function for loading device constants (acLoadDeviceConstant).
+// This acForcingVec should go outside the core library since it references user-defined
+// parameters such as AC_forcing_magnitude which may not be defined in all projects.
+// host_forcing.cc is probably a good place for this.
+// %JP update: moved this here out of astaroth.cu. Should be renamed such that it cannot be
+// confused with actual interface functions.
+// %JP update 2: deprecated acForcingVec: use loadForcingParams instead
+void
+DEPRECATED_acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force,
+ const AcReal3 ff_hel_re, const AcReal3 ff_hel_im,
+ const AcReal forcing_phase, const AcReal kaver)
+{
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_forcing_magnitude, forcing_magnitude);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_forcing_phase, forcing_phase);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_k_forcex, k_force.x);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_k_forcey, k_force.y);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_k_forcez, k_force.z);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_rex, ff_hel_re.x);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_rey, ff_hel_re.y);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_rez, ff_hel_re.z);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_imx, ff_hel_im.x);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_imy, ff_hel_im.y);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_imz, ff_hel_im.z);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_kaver, kaver);
+}
+
+void
+loadForcingParamsToGrid(const ForcingParams& forcing_params)
+{
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_forcing_magnitude, forcing_params.magnitude);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_forcing_phase, forcing_params.phase);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_k_forcex, forcing_params.k_force.x);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_k_forcey, forcing_params.k_force.y);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_k_forcez, forcing_params.k_force.z);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_rex, forcing_params.ff_hel_re.x);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_rey, forcing_params.ff_hel_re.y);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_rez, forcing_params.ff_hel_re.z);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_imx, forcing_params.ff_hel_im.x);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_imy, forcing_params.ff_hel_im.y);
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_ff_hel_imz, forcing_params.ff_hel_im.z);
+
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_kaver, forcing_params.kaver);
+ acGridSynchronizeStream(STREAM_ALL);
+}
+
+/** This function would be used in autotesting to update the forcing params of the host
+ configuration. */
+void
+loadForcingParamsToHost(const ForcingParams& forcing_params, AcMesh* mesh)
+{
+ // %JP: Left some regex magic here in case we need to modify the ForcingParams struct
+ // acLoadDeviceConstant\(([A-Za-z_]*), ([a-z_.]*)\);
+ // mesh->info.real_params[$1] = $2;
+ mesh->info.real_params[AC_forcing_magnitude] = forcing_params.magnitude;
+ mesh->info.real_params[AC_forcing_phase] = forcing_params.phase;
+
+ mesh->info.real_params[AC_k_forcex] = forcing_params.k_force.x;
+ mesh->info.real_params[AC_k_forcey] = forcing_params.k_force.y;
+ mesh->info.real_params[AC_k_forcez] = forcing_params.k_force.z;
+
+ mesh->info.real_params[AC_ff_hel_rex] = forcing_params.ff_hel_re.x;
+ mesh->info.real_params[AC_ff_hel_rey] = forcing_params.ff_hel_re.y;
+ mesh->info.real_params[AC_ff_hel_rez] = forcing_params.ff_hel_re.z;
+
+ mesh->info.real_params[AC_ff_hel_imx] = forcing_params.ff_hel_im.x;
+ mesh->info.real_params[AC_ff_hel_imy] = forcing_params.ff_hel_im.y;
+ mesh->info.real_params[AC_ff_hel_imz] = forcing_params.ff_hel_im.z;
+
+ mesh->info.real_params[AC_kaver] = forcing_params.kaver;
+}
+
+ForcingParams
+generateForcingParams(const AcMeshInfo& mesh_info)
+{
+ ForcingParams params = {};
+
+ // Forcing properties
+ AcReal relhel = mesh_info.real_params[AC_relhel];
+ params.magnitude = mesh_info.real_params[AC_forcing_magnitude];
+ AcReal kmin = mesh_info.real_params[AC_kmin];
+ AcReal kmax = mesh_info.real_params[AC_kmax];
+
+ params.kaver = (kmax - kmin) / AcReal(2.0);
+
+ // Generate forcing wave vector k_force
+ params.k_force = helical_forcing_k_generator(kmax, kmin);
+
+ // Randomize the phase
+ params.phase = AcReal(2.0) * AcReal(M_PI) * get_random_number_01();
+
+ // Generate e for k. Needed for the sake of isotrophy.
+ AcReal3 e_force;
+ if ((params.k_force.y == AcReal(0.0)) && (params.k_force.z == AcReal(0.0))) {
+ e_force = (AcReal3){0.0, 1.0, 0.0};
+ }
+ else {
+ e_force = (AcReal3){1.0, 0.0, 0.0};
+ }
+ helical_forcing_e_generator(&e_force, params.k_force);
+
+ helical_forcing_special_vector(¶ms.ff_hel_re, ¶ms.ff_hel_im, params.k_force, e_force,
+ relhel);
+
+ return params;
+}
diff --git a/samples/standalone_mpi/host_forcing.h b/samples/standalone_mpi/host_forcing.h
new file mode 100644
index 0000000..14b77aa
--- /dev/null
+++ b/samples/standalone_mpi/host_forcing.h
@@ -0,0 +1,60 @@
+
+/*
+ Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+
+/**
+ * @file
+ * \brief Brief info.
+ *
+ * Detailed info.
+ *
+ */
+#pragma once
+#include "astaroth.h"
+
+AcReal get_random_number_01();
+
+AcReal3 helical_forcing_k_generator(const AcReal kmax, const AcReal kmin);
+
+void helical_forcing_e_generator(AcReal3* e_force, const AcReal3 k_force);
+
+void helical_forcing_special_vector(AcReal3* ff_hel_re, AcReal3* ff_hel_im, const AcReal3 k_force,
+ const AcReal3 e_force, const AcReal relhel);
+
+/** Tool for loading forcing vector information into the device memory
+ // DEPRECATED in favour of loadForcingParams
+ */
+void DEPRECATED_acForcingVec(const AcReal forcing_magnitude, const AcReal3 k_force,
+ const AcReal3 ff_hel_re, const AcReal3 ff_hel_im,
+ const AcReal forcing_phase, const AcReal kaver);
+
+typedef struct {
+ AcReal magnitude;
+ AcReal3 k_force;
+ AcReal3 ff_hel_re;
+ AcReal3 ff_hel_im;
+ AcReal phase;
+ AcReal kaver;
+} ForcingParams;
+
+void loadForcingParamsToGrid(const ForcingParams& forcing_params);
+
+void loadForcingParamsToHost(const ForcingParams& forcing_params, AcMesh* mesh);
+
+ForcingParams generateForcingParams(const AcMeshInfo& mesh_info);
diff --git a/samples/standalone_mpi/host_memory.cc b/samples/standalone_mpi/host_memory.cc
new file mode 100644
index 0000000..1e809df
--- /dev/null
+++ b/samples/standalone_mpi/host_memory.cc
@@ -0,0 +1,718 @@
+/*
+ Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+
+/**
+ * @file
+ * \brief Brief info.
+ *
+ * Detailed info.
+ *
+ */
+#include "host_memory.h"
+
+#include
+
+#include "astaroth_utils.h"
+#include "errchk.h"
+
+#define AC_GEN_STR(X) #X
+const char* init_type_names[] = {AC_FOR_INIT_TYPES(AC_GEN_STR)};
+#undef AC_GEN_STR
+
+#define XORIG (AcReal(.5) * mesh->info.int_params[AC_nx] * mesh->info.real_params[AC_dsx])
+#define YORIG (AcReal(.5) * mesh->info.int_params[AC_ny] * mesh->info.real_params[AC_dsy])
+#define ZORIG (AcReal(.5) * mesh->info.int_params[AC_nz] * mesh->info.real_params[AC_dsz])
+
+static AcReal
+randr(void)
+{
+ return AcReal(rand()) / AcReal(RAND_MAX);
+}
+
+void
+lnrho_step(AcMesh* mesh)
+{
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+ const int mz = mesh->info.int_params[AC_mz];
+
+ // const int nx_min = mesh->info.int_params[AC_nx_min];
+ // const int nx_max = mesh->info.int_params[AC_nx_max];
+ // const int ny_min = mesh->info.int_params[AC_ny_min];
+ // const int ny_max = mesh->info.int_params[AC_ny_max];
+ // const int nz_min = mesh->info.int_params[AC_nz_min];
+ // const int nz_max = mesh->info.int_params[AC_nz_max];
+
+ // const AcReal DX = mesh->info.real_params[AC_dsx];
+ // const AcReal DY = mesh->info.real_params[AC_dsy];
+ // const AcReal DZ = mesh->info.real_params[AC_dsz];
+ // const AcReal xmax = DX * (nx_max - nx_min) ;
+ // const AcReal zmax = DZ * (nz_max - nz_min) ;
+
+ // const AcReal lnrho1 = (AcReal) -1.0; // TODO mesh->info.real_params[AC_lnrho1];
+ const AcReal lnrho2 = (AcReal)0.0; // TODO mesh->info.real_params[AC_lnrho2];
+ // const AcReal rho1 = (AcReal) exp(lnrho1);
+ // const AcReal rho2 = (AcReal) exp(lnrho2);
+
+ // const AcReal k_pert = (AcReal) 1.0; //mesh->info.real_params[AC_k_pert]; //Wamenumber of
+ // the perturbation const AcReal k_pert = 4.0; //mesh->info.real_params[AC_k_pert];
+ // //Wamenumber of the perturbation
+ // const AcReal ampl_pert = xmax/10.0; // xmax/mesh->info.real_params[AC_pert]; //Amplitude of
+ // the perturbation
+ // const AcReal ampl_pert = (AcReal) 0.0;//xmax/20.0; // xmax/mesh->info.real_params[AC_pert];
+ // //Amplitude of the perturbation const AcReal two_pi = (AcReal) 6.28318531;
+
+ // const AcReal xorig = mesh->info.real_params[AC_xorig];
+ // const AcReal zorig = mesh->info.real_params[AC_zorig];
+ // const AcReal trans = mesh->info.real_params[AC_trans];
+
+ // AcReal xx, zz, tanhprof, cosz_wave;
+
+ for (int k = 0; k < mz; k++) {
+ for (int j = 0; j < my; j++) {
+ for (int i = 0; i < mx; i++) {
+ int idx = i + j * mx + k * mx * my;
+ // zz = DZ * AcReal(k) - zorig; // Not used
+ // cosz_wave = ampl_pert*AcReal(cos(k_pert*((zz/zmax)*two_pi))); // Not used
+ // xx = DX * AcReal(i) - xorig + cosz_wave; //ADD WAVE TODO // Not used
+ // tanhprof = AcReal(0.5)*((rho2+rho1) + (rho2-rho1)*AcReal(tanh(xx/trans))); // Not
+ // used Commented out the step function initial codition.
+ // mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(tanhprof);
+ mesh->vertex_buffer[VTXBUF_LNRHO][idx] = lnrho2;
+ }
+ }
+ }
+}
+
+// This is the initial condition type for the infalling vedge in the pseudodisk
+// model.
+void
+inflow_vedge(AcMesh* mesh)
+{
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+ const int mz = mesh->info.int_params[AC_mz];
+
+ // const int nx_min = mesh->info.int_params[AC_nx_min];
+ // const int nx_max = mesh->info.int_params[AC_nx_max];
+ // const int ny_min = mesh->info.int_params[AC_ny_min];
+ // const int ny_max = mesh->info.int_params[AC_ny_max];
+ // const int nz_min = mesh->info.int_params[AC_nz_min];
+ // const int nz_max = mesh->info.int_params[AC_nz_max];
+
+ // const double DX = mesh->info.real_params[AC_dsx];
+ // const double DY = mesh->info.real_params[AC_dsy];
+ const double DZ = mesh->info.real_params[AC_dsz];
+
+ const double AMPL_UU = mesh->info.real_params[AC_ampl_uu];
+ const double ANGL_UU = mesh->info.real_params[AC_angl_uu];
+
+ const double zorig = mesh->info.real_params[AC_zorig];
+ double zz;
+ double trans = mesh->info.real_params[AC_trans];
+
+ // const AcReal range = AcReal(.5);
+
+ // const AcReal zmax = AcReal(DZ * (nz_max - nz_min));
+ // const AcReal gaussr = zmax / AcReal(4.0);
+
+ // for (int k = nz_min; k < nz_max; k++) {
+ // for (int j = ny_min; j < ny_max; j++) {
+ // for (int i = nx_min; i < nx_max; i++) {
+ for (int k = 0; k < mz; k++) {
+ for (int j = 0; j < my; j++) {
+ for (int i = 0; i < mx; i++) {
+ int idx = i + j * mx + k * mx * my;
+ zz = DZ * double(k) - zorig;
+ // mesh->vertex_buffer[VTXBUF_UUX][idx] = -AMPL_UU*cos(ANGL_UU);
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(-AMPL_UU * cos(ANGL_UU) *
+ fabs(tanh(zz / trans)));
+ mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal(0.0);
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(-AMPL_UU * sin(ANGL_UU) *
+ tanh(zz / trans));
+
+ // Variarion to density
+ // AcReal rho = exp(mesh->vertex_buffer[VTXBUF_LNRHO][idx]);
+ // NO GAUSSIAN//rho = rho*exp(-(zz/gaussr)*(zz/gaussr));
+ // mesh->vertex_buffer[VTXBUF_LNRHO][idx] = log(rho + (range*rho) * (randr() -
+ // AcReal(-0.5)));
+ }
+ }
+ }
+}
+
+// This is the initial condition type for the infalling vedge in the pseudodisk
+// model.
+void
+simple_uniform_core(AcMesh* mesh)
+{
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+ const int mz = mesh->info.int_params[AC_mz];
+
+ const double DX = mesh->info.real_params[AC_dsx];
+ const double DY = mesh->info.real_params[AC_dsy];
+ const double DZ = mesh->info.real_params[AC_dsz];
+
+ const double ampl_lnrho = mesh->info.real_params[AC_ampl_lnrho];
+
+ const double xorig = mesh->info.real_params[AC_xorig];
+ const double yorig = mesh->info.real_params[AC_yorig];
+ const double zorig = mesh->info.real_params[AC_zorig];
+
+ const double G_const = mesh->info.real_params[AC_G_const];
+ const double M_sink_init = mesh->info.real_params[AC_M_sink_init];
+ const double cs2_sound = mesh->info.real_params[AC_cs2_sound];
+
+ const double RR_inner_bound = mesh->info.real_params[AC_soft] / AcReal(2.0);
+ const double core_coeff = (exp(ampl_lnrho) * cs2_sound) / (double(4.0) * M_PI * G_const);
+
+ double xx, yy, zz, RR;
+ double core_profile;
+
+ // TEMPORARY TEST INPUT PARAMETERS
+ const double core_radius = DX * 32.0;
+ const double trans = DX * 12.0;
+ // const double epsilon = DX*2.0;
+ const double vel_scale = mesh->info.real_params[AC_ampl_uu];
+ double abso_vel;
+
+ RR = 1.0;
+ printf("%e %e %e \n", RR, trans, core_radius);
+
+ for (int k = 0; k < mz; k++) {
+ for (int j = 0; j < my; j++) {
+ for (int i = 0; i < mx; i++) {
+ int idx = i + j * mx + k * mx * my;
+ xx = DX * double(i) - xorig;
+ yy = DY * double(j) - yorig;
+ zz = DZ * double(k) - zorig;
+
+ RR = sqrt(xx * xx + yy * yy + zz * zz);
+
+ if (RR >= RR_inner_bound) {
+ abso_vel = vel_scale * sqrt(2.0 * G_const * M_sink_init / RR);
+ core_profile = pow(RR, -2.0); // double(1.0);
+ }
+ else {
+ abso_vel = vel_scale * sqrt(2.0 * AC_G_const * AC_M_sink_init / RR_inner_bound);
+ core_profile = pow(RR_inner_bound, -2.0); // double(1.0);
+ }
+
+ if (RR <= sqrt(DX * DX + DY * DY + DZ * DZ)) {
+ abso_vel = 0.0;
+ RR = 1.0;
+ }
+
+ mesh->vertex_buffer[VTXBUF_LNRHO][idx] = AcReal(log(core_coeff * core_profile));
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(-abso_vel * (yy / RR));
+ mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal(abso_vel * (xx / RR));
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(0.0);
+ }
+ }
+ }
+}
+
+// This is the initial condition type for the infalling vedge in the pseudodisk
+// model.
+void
+inflow_vedge_freefall(AcMesh* mesh)
+{
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+ const int mz = mesh->info.int_params[AC_mz];
+
+ // const int nx_min = mesh->info.int_params[AC_nx_min];
+ // const int nx_max = mesh->info.int_params[AC_nx_max];
+ // const int ny_min = mesh->info.int_params[AC_ny_min];
+ // const int ny_max = mesh->info.int_params[AC_ny_max];
+ // const int nz_min = mesh->info.int_params[AC_nz_min];
+ // const int nz_max = mesh->info.int_params[AC_nz_max];
+
+ const double DX = mesh->info.real_params[AC_dsx];
+ // const double DY = mesh->info.real_params[AC_dsy];
+ const double DZ = mesh->info.real_params[AC_dsz];
+
+ // const double AMPL_UU = mesh->info.real_params[AC_ampl_uu];
+ const double ANGL_UU = mesh->info.real_params[AC_angl_uu];
+ const double SQ2GM = mesh->info.real_params[AC_sq2GM_star];
+ // const double GM = mesh->info.real_params[AC_GM_star];
+ // const double M_star = mesh->info.real_params[AC_M_star];
+ // const double G_CONST = mesh->info.real_params[AC_G_CONST];
+
+ // const double unit_length = mesh->info.real_params[AC_unit_length];
+ // const double unit_density = mesh->info.real_params[AC_unit_density];
+ // const double unit_velocity = mesh->info.real_params[AC_unit_velocity];
+
+ const double xorig = mesh->info.real_params[AC_xorig];
+ // const double yorig = mesh->info.real_params[AC_yorig];
+ const double zorig = mesh->info.real_params[AC_zorig];
+ // const double trans = mesh->info.real_params[AC_trans];
+ // double xx, yy, zz, RR;
+ double xx, zz, RR;
+ // double delx, dely, delz;
+ double delx, delz;
+ // double u_x, u_y, u_z, veltot, tanhz;
+ double u_x, u_z, veltot, tanhz;
+
+ const double star_pos_x = mesh->info.real_params[AC_star_pos_x];
+ const double star_pos_z = mesh->info.real_params[AC_star_pos_z];
+
+ for (int k = 0; k < mz; k++) {
+ for (int j = 0; j < my; j++) {
+ for (int i = 0; i < mx; i++) {
+ int idx = i + j * mx + k * mx * my;
+ xx = DX * double(i) - xorig;
+ zz = DZ * double(k) - zorig;
+
+ delx = xx - star_pos_x;
+ delz = zz - star_pos_z;
+ // TODO: Figure out isthis needed. Now a placeholder.
+ // tanhz = fabs(tanh(zz/trans));
+ tanhz = 1.0;
+
+ RR = sqrt(delx * delx + delz * delz);
+ veltot = SQ2GM / sqrt(RR); // Free fall velocity
+
+ // Normal velocity components
+ u_x = -veltot * (delx / RR);
+ u_z = -veltot * (delz / RR);
+
+ // printf("star_pos_z %e, zz %e, delz %e, RR %e\n", star_pos_z, zz, delz, RR);
+
+ // printf("unit_length = %e, unit_density = %e, unit_velocity = %e,\n M_star = %e,
+ // G_CONST = %e, GM = %e, SQ2GM = %e, \n RR = %e, u_x = %e, u_z %e\n",
+ // unit_length, unit_density,
+ // unit_velocity, M_star, G_CONST, GM, SQ2GM, RR, u_x, u_z);
+ // printf("%e\n", unit_length*unit_length*unit_length);
+
+ // Here including an angel tilt due to pseudodisk
+ if (delz >= 0.0) {
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(
+ (u_x * cos(ANGL_UU) - u_z * sin(ANGL_UU)) * tanhz);
+ mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal(0.0);
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(
+ (u_x * sin(ANGL_UU) + u_z * cos(ANGL_UU)) * tanhz);
+ }
+ else {
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = AcReal(
+ (u_x * cos(ANGL_UU) + u_z * sin(ANGL_UU)) * tanhz);
+ mesh->vertex_buffer[VTXBUF_UUY][idx] = AcReal(0.0);
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = AcReal(
+ (-u_x * sin(ANGL_UU) + u_z * cos(ANGL_UU)) * tanhz);
+ }
+ }
+ }
+ }
+}
+
+// Only x-direction free fall
+void
+inflow_freefall_x(AcMesh* mesh)
+{
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+ const int mz = mesh->info.int_params[AC_mz];
+
+ const double DX = mesh->info.real_params[AC_dsx];
+
+ const double SQ2GM = mesh->info.real_params[AC_sq2GM_star];
+ // const double G_CONST = mesh->info.real_params[AC_G_CONST];
+
+ const double xorig = mesh->info.real_params[AC_xorig];
+ double xx, RR;
+ double delx;
+ double /*u_x,*/ veltot;
+
+ const double star_pos_x = mesh->info.real_params[AC_star_pos_x];
+
+ const double ampl_lnrho = mesh->info.real_params[AC_ampl_lnrho];
+
+ for (int k = 0; k < mz; k++) {
+ for (int j = 0; j < my; j++) {
+ for (int i = 0; i < mx; i++) {
+ int idx = i + j * mx + k * mx * my;
+ xx = DX * double(i) - xorig;
+
+ delx = xx - star_pos_x;
+
+ RR = fabs(delx);
+
+ veltot = SQ2GM / sqrt(RR); // Free fall velocity
+
+ if (isinf(veltot) == 1)
+ printf("xx %e star_pos_x %e delz %e RR %e veltot %e\n", xx, star_pos_x, delx,
+ RR, veltot);
+
+ // Normal velocity components
+ // u_x = - veltot; // Not used
+
+ // Freefall condition
+ // mesh->vertex_buffer[VTXBUF_UUX][idx] = u_x;
+ // mesh->vertex_buffer[VTXBUF_UUY][idx] = 0.0;
+ // mesh->vertex_buffer[VTXBUF_UUZ][idx] = 0.0;
+
+ // Starting with steady state
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = 0.0;
+ mesh->vertex_buffer[VTXBUF_UUY][idx] = 0.0;
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = 0.0;
+
+ mesh->vertex_buffer[VTXBUF_LNRHO][idx] = AcReal(ampl_lnrho);
+ }
+ }
+ }
+}
+
+void
+gaussian_radial_explosion(AcMesh* mesh)
+{
+ AcReal* uu_x = mesh->vertex_buffer[VTXBUF_UUX];
+ AcReal* uu_y = mesh->vertex_buffer[VTXBUF_UUY];
+ AcReal* uu_z = mesh->vertex_buffer[VTXBUF_UUZ];
+
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+
+ const int nx_min = mesh->info.int_params[AC_nx_min];
+ const int nx_max = mesh->info.int_params[AC_nx_max];
+ const int ny_min = mesh->info.int_params[AC_ny_min];
+ const int ny_max = mesh->info.int_params[AC_ny_max];
+ const int nz_min = mesh->info.int_params[AC_nz_min];
+ const int nz_max = mesh->info.int_params[AC_nz_max];
+
+ const double DX = mesh->info.real_params[AC_dsx];
+ const double DY = mesh->info.real_params[AC_dsy];
+ const double DZ = mesh->info.real_params[AC_dsz];
+
+ const double xorig = double(XORIG) - 0.000001;
+ const double yorig = double(YORIG) - 0.000001;
+ const double zorig = double(ZORIG) - 0.000001;
+
+ const double INIT_LOC_UU_X = 0.0;
+ const double INIT_LOC_UU_Y = 0.0;
+ const double INIT_LOC_UU_Z = 0.0;
+
+ const double AMPL_UU = mesh->info.real_params[AC_ampl_uu];
+ const double UU_SHELL_R = 0.8;
+ const double WIDTH_UU = 0.2;
+
+ // Outward explosion with gaussian initial velocity profile.
+ int idx;
+ double xx, yy, zz, rr2, rr, theta = 0.0, phi = 0.0;
+ double uu_radial;
+
+ // double theta_old = 0.0;
+
+ for (int k = nz_min; k < nz_max; k++) {
+ for (int j = ny_min; j < ny_max; j++) {
+ for (int i = nx_min; i < nx_max; i++) {
+ // Calculate the value of velocity in a particular radius.
+ idx = i + j * mx + k * mx * my;
+ // Determine the coordinates
+ xx = DX * (i - nx_min) - xorig;
+ xx = xx - INIT_LOC_UU_X;
+
+ yy = DY * (j - ny_min) - yorig;
+ yy = yy - INIT_LOC_UU_Y;
+
+ zz = DZ * (k - nz_min) - zorig;
+ zz = zz - INIT_LOC_UU_Z;
+
+ rr2 = pow(xx, 2.0) + pow(yy, 2.0) + pow(zz, 2.0);
+ rr = sqrt(rr2);
+
+ // Origin is different!
+ double xx_abs, yy_abs, zz_abs;
+ if (rr > 0.0) {
+ // theta range [0, PI]
+ if (zz >= 0.0) {
+ theta = acos(zz / rr);
+ if (theta > M_PI / 2.0 || theta < 0.0) {
+ printf("Explosion THETA WRONG: zz = %.3f, rr = "
+ "%.3f, theta = %.3e/PI, M_PI = %.3e\n",
+ zz, rr, theta / M_PI, M_PI);
+ }
+ }
+ else {
+ zz_abs = -zz; // Needs a posite value for acos
+ theta = M_PI - acos(zz_abs / rr);
+ if (theta < M_PI / 2.0 || theta > 2 * M_PI) {
+ printf("Explosion THETA WRONG: zz = %.3f, rr = "
+ "%.3f, theta = %.3e/PI, M_PI = %.3e\n",
+ zz, rr, theta / M_PI, M_PI);
+ }
+ }
+
+ // phi range [0, 2*PI]i
+ if (xx != 0.0) {
+ if (xx < 0.0 && yy >= 0.0) {
+ //-+
+ xx_abs = -xx; // Needs a posite value for atan
+ phi = M_PI - atan(yy / xx_abs);
+ if (phi < (M_PI / 2.0) || phi > M_PI) {
+ printf("Explosion PHI WRONG -+: xx = %.3f, yy "
+ "= %.3f, phi = %.3e/PI, M_PI = %.3e\n",
+ xx, yy, phi / M_PI, M_PI);
+ }
+ }
+ else if (xx > 0.0 && yy < 0.0) {
+ //+-
+ yy_abs = -yy;
+ phi = 2.0 * M_PI - atan(yy_abs / xx);
+ if (phi < (3.0 * M_PI) / 2.0 || phi > (2.0 * M_PI + 1e-6)) {
+ printf("Explosion PHI WRONG +-: xx = %.3f, yy "
+ "= %.3f, phi = %.3e/PI, M_PI = %.3e\n",
+ xx, yy, phi / M_PI, M_PI);
+ }
+ }
+ else if (xx < 0.0 && yy < 0.0) {
+ //--
+ yy_abs = -yy;
+ xx_abs = -xx;
+ phi = M_PI + atan(yy_abs / xx_abs);
+ if (phi < M_PI || phi > ((3.0 * M_PI) / 2.0 + 1e-6)) {
+ printf("Explosion PHI WRONG --: xx = %.3f, yy "
+ "= %.3f, xx_abs = %.3f, yy_abs = %.3f, "
+ "phi = %.3e, (3.0*M_PI)/2.0 = %.3e\n",
+ xx, yy, xx_abs, yy_abs, phi, (3.0 * M_PI) / 2.0);
+ }
+ }
+ else {
+ //++
+ phi = atan(yy / xx);
+ if (phi < 0 || phi > M_PI / 2.0) {
+ printf("Explosion PHI WRONG --: xx = %.3f, yy = "
+ "%.3f, phi = %.3e, (3.0*M_PI)/2.0 = %.3e\n",
+ xx, yy, phi, (3.0 * M_PI) / 2.0);
+ }
+ }
+ }
+ else { // To avoid div by zero with atan
+ if (yy > 0.0) {
+ phi = M_PI / 2.0;
+ }
+ else if (yy < 0.0) {
+ phi = (3.0 * M_PI) / 2.0;
+ }
+ else {
+ phi = 0.0;
+ }
+ }
+
+ // Set zero for explicit safekeeping
+ if (xx == 0.0 && yy == 0.0) {
+ phi = 0.0;
+ }
+
+ // Gaussian velocity
+ // uu_radial = AMPL_UU*exp( -rr2 / (2.0*pow(WIDTH_UU, 2.0))
+ // ); New distribution, where that gaussion wave is not in
+ // the exact centre coordinates uu_radial = AMPL_UU*exp(
+ // -pow((rr - 4.0*WIDTH_UU),2.0) / (2.0*pow(WIDTH_UU, 2.0))
+ // ); //TODO: Parametrize the peak location.
+ uu_radial = AMPL_UU *
+ exp(-pow((rr - UU_SHELL_R), 2.0) / (2.0 * pow(WIDTH_UU, 2.0)));
+ }
+ else {
+ uu_radial = 0.0; // TODO: There will be a discontinuity in
+ // the origin... Should the shape of the
+ // distribution be different?
+ }
+
+ // Determine the carthesian velocity components and lnrho
+ uu_x[idx] = AcReal(uu_radial * sin(theta) * cos(phi));
+ uu_y[idx] = AcReal(uu_radial * sin(theta) * sin(phi));
+ uu_z[idx] = AcReal(uu_radial * cos(theta));
+
+ // Temporary diagnosticv output (TODO: Remove after not needed)
+ // if (theta > theta_old) {
+ // if (theta > M_PI || theta < 0.0 || phi < 0.0 || phi > 2*M_PI)
+ // {
+ /* printf("Explosion: xx = %.3f, yy = %.3f, zz = %.3f, rr =
+ %.3f, phi = %.3e/PI, theta = %.3e/PI\n, M_PI = %.3e", xx, yy,
+ zz, rr, phi/M_PI, theta/M_PI, M_PI); printf(" uu_radial =
+ %.3e, uu_x[%i] = %.3e, uu_y[%i] = %.3e, uu_z[%i] = %.3e \n",
+ uu_radial, idx, uu_x[idx], idx, uu_y[idx], idx,
+ uu_z[idx]); theta_old = theta;
+ */
+ }
+ }
+ }
+}
+
+void
+acmesh_init_to(const InitType& init_type, AcMesh* mesh)
+{
+ srand(123456789);
+
+ const int n = acVertexBufferSize(mesh->info);
+
+ const int mx = mesh->info.int_params[AC_mx];
+ const int my = mesh->info.int_params[AC_my];
+ const int mz = mesh->info.int_params[AC_mz];
+
+ const int nx_min = mesh->info.int_params[AC_nx_min];
+ const int nx_max = mesh->info.int_params[AC_nx_max];
+ const int ny_min = mesh->info.int_params[AC_ny_min];
+ const int ny_max = mesh->info.int_params[AC_ny_max];
+ const int nz_min = mesh->info.int_params[AC_nz_min];
+ const int nz_max = mesh->info.int_params[AC_nz_max];
+
+ switch (init_type) {
+ case INIT_TYPE_RANDOM: {
+ acMeshClear(mesh);
+ const AcReal range = AcReal(0.01);
+ for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w)
+ for (int i = 0; i < n; ++i)
+ mesh->vertex_buffer[w][i] = 2 * range * randr() - range;
+
+ break;
+ }
+ case INIT_TYPE_GAUSSIAN_RADIAL_EXPL:
+ acMeshClear(mesh);
+ acVertexBufferSet(VTXBUF_LNRHO, mesh->info.real_params[AC_ampl_lnrho], mesh);
+ // acmesh_init_to(INIT_TYPE_RANDOM, mesh);
+ gaussian_radial_explosion(mesh);
+
+ break;
+ case INIT_TYPE_XWAVE:
+ acMeshClear(mesh);
+ acmesh_init_to(INIT_TYPE_RANDOM, mesh);
+ for (int k = 0; k < mz; k++) {
+ for (int j = 0; j < my; j++) {
+ for (int i = 0; i < mx; i++) {
+ int idx = i + j * mx + k * mx * my;
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = 2 * AcReal(sin(j * AcReal(M_PI) / mx)) -
+ 1;
+ }
+ }
+ }
+ break;
+ case INIT_TYPE_SIMPLE_CORE:
+ acMeshClear(mesh);
+ simple_uniform_core(mesh);
+ break;
+ case INIT_TYPE_VEDGE:
+ acMeshClear(mesh);
+ inflow_vedge_freefall(mesh);
+ break;
+ case INIT_TYPE_VEDGEX:
+ acMeshClear(mesh);
+ inflow_freefall_x(mesh);
+ break;
+ case INIT_TYPE_RAYLEIGH_TAYLOR:
+ acMeshClear(mesh);
+ inflow_freefall_x(mesh);
+ lnrho_step(mesh);
+ break;
+ case INIT_TYPE_ABC_FLOW: {
+ acMeshClear(mesh);
+ acmesh_init_to(INIT_TYPE_RANDOM, mesh);
+ for (int k = nz_min; k < nz_max; k++) {
+ for (int j = ny_min; j < ny_max; j++) {
+ for (int i = nx_min; i < nx_max; i++) {
+ const int idx = i + j * mx + k * mx * my;
+
+ /*
+ const double xx = double(
+ mesh->info.real_params[AC_dsx] *
+ (i - mesh->info.int_params[AC_nx_min]) -
+ XORIG + AcReal(.5) * mesh->info.real_params[AC_dsx]);
+ const double yy = double(
+ mesh->info.real_params[AC_dsy] *
+ (j - mesh->info.int_params[AC_ny_min]) -
+ YORIG + AcReal(.5) * mesh->info.real_params[AC_dsy]);
+ const double zz = double(
+ mesh->info.real_params[AC_dsz] *
+ (k - mesh->info.int_params[AC_nz_min]) -
+ ZORIG + AcReal(.5) * mesh->info.real_params[AC_dsz]);
+ */
+
+ const AcReal xx = (i - nx_min) * mesh->info.real_params[AC_dsx] - XORIG;
+ const AcReal yy = (j - ny_min) * mesh->info.real_params[AC_dsy] - YORIG;
+ const AcReal zz = (k - nz_min) * mesh->info.real_params[AC_dsz] - ZORIG;
+
+ const AcReal ampl_uu = 0.5;
+ const AcReal ABC_A = 1.;
+ const AcReal ABC_B = 1.;
+ const AcReal ABC_C = 1.;
+ const AcReal kx_uu = 8.;
+ const AcReal ky_uu = 8.;
+ const AcReal kz_uu = 8.;
+
+ mesh->vertex_buffer[VTXBUF_UUX][idx] = ampl_uu *
+ (ABC_A * (AcReal)sin(kz_uu * zz) +
+ ABC_C * (AcReal)cos(ky_uu * yy));
+ mesh->vertex_buffer[VTXBUF_UUY][idx] = ampl_uu *
+ (ABC_B * (AcReal)sin(kx_uu * xx) +
+ ABC_A * (AcReal)cos(kz_uu * zz));
+ mesh->vertex_buffer[VTXBUF_UUZ][idx] = ampl_uu *
+ (ABC_C * (AcReal)sin(ky_uu * yy) +
+ ABC_B * (AcReal)cos(kx_uu * xx));
+ }
+ }
+ }
+ break;
+ }
+ case INIT_TYPE_RAYLEIGH_BENARD: {
+ acmesh_init_to(INIT_TYPE_RANDOM, mesh);
+#if LTEMPERATURE
+ acVertexBufferSet(VTXBUF_LNRHO, 1, mesh);
+ const AcReal range = AcReal(0.9);
+ for (int k = nz_min; k < nz_max; k++) {
+ for (int j = ny_min; j < ny_max; j++) {
+ for (int i = nx_min; i < nx_max; i++) {
+ const int idx = i + j * mx + k * mx * my;
+ mesh->vertex_buffer[VTXBUF_TEMPERATURE][idx] = (range * (k - nz_min)) /
+ mesh->info
+ .int_params[AC_nz] +
+ 0.1;
+ }
+ }
+ }
+#else
+ WARNING("INIT_TYPE_RAYLEIGH_BERNARD called even though VTXBUF_TEMPERATURE is not used");
+#endif
+ break;
+ }
+ default:
+ ERROR("Unknown init_type");
+ }
+
+ AcReal max_val = AcReal(-1e-32);
+ AcReal min_val = AcReal(1e32);
+ // Normalize the grid
+ for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
+ for (int i = 0; i < n; ++i) {
+ if (mesh->vertex_buffer[w][i] < min_val)
+ min_val = mesh->vertex_buffer[w][i];
+ if (mesh->vertex_buffer[w][i] > max_val)
+ max_val = mesh->vertex_buffer[w][i];
+ }
+ }
+ printf("MAX: %f MIN %f\n", double(max_val), double(min_val));
+ /*
+ const AcReal inv_range = AcReal(1.) / fabs(max_val - min_val);
+ for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
+ for (int i = 0; i < n; ++i) {
+ mesh->vertex_buffer[w][i] = 2*inv_range*(mesh->vertex_buffer[w][i] - min_val) - 1;
+ }
+ }
+ */
+}
diff --git a/samples/standalone_mpi/host_memory.h b/samples/standalone_mpi/host_memory.h
new file mode 100644
index 0000000..f0e3181
--- /dev/null
+++ b/samples/standalone_mpi/host_memory.h
@@ -0,0 +1,49 @@
+/*
+ Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+
+/**
+ * @file
+ * \brief Brief info.
+ *
+ * Detailed info.
+ *
+ */
+#pragma once
+#include "astaroth.h"
+
+// clang-format off
+#define AC_FOR_INIT_TYPES(FUNC)\
+ FUNC(INIT_TYPE_RANDOM), \
+ FUNC(INIT_TYPE_XWAVE), \
+ FUNC(INIT_TYPE_GAUSSIAN_RADIAL_EXPL), \
+ FUNC(INIT_TYPE_ABC_FLOW) , \
+ FUNC(INIT_TYPE_SIMPLE_CORE), \
+ FUNC(INIT_TYPE_VEDGE), \
+ FUNC(INIT_TYPE_VEDGEX), \
+ FUNC(INIT_TYPE_RAYLEIGH_TAYLOR), \
+ FUNC(INIT_TYPE_RAYLEIGH_BENARD)
+// clang-format on
+
+#define AC_GEN_ID(X) X
+typedef enum { AC_FOR_INIT_TYPES(AC_GEN_ID), NUM_INIT_TYPES } InitType;
+#undef AC_GEN_ID
+
+extern const char* init_type_names[]; // Defined in host_memory.cc
+
+void acmesh_init_to(const InitType& type, AcMesh* mesh);
diff --git a/samples/standalone_mpi/main.cc b/samples/standalone_mpi/main.cc
new file mode 100644
index 0000000..c5bb51a
--- /dev/null
+++ b/samples/standalone_mpi/main.cc
@@ -0,0 +1,586 @@
+/*
+ Copyright (C) 2014-2020, Johannes Pekkila, Miikka Vaisala.
+
+ This file is part of Astaroth.
+
+ Astaroth is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ Astaroth is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Astaroth. If not, see .
+*/
+/**
+ Running: mpirun -np
+*/
+#if AC_MPI_ENABLED
+#include "astaroth.h"
+#include "astaroth_utils.h"
+
+#include
+#include
+
+#include "config_loader.h"
+#include "errchk.h"
+#include "host_forcing.h"
+#include "host_memory.h"
+
+#define fprintf(...) \
+ { \
+ int tmppid; \
+ MPI_Comm_rank(MPI_COMM_WORLD, &tmppid); \
+ if (!tmppid) { \
+ fprintf(__VA_ARGS__); \
+ } \
+ }
+
+#define printf(...) \
+ { \
+ int tmppid; \
+ MPI_Comm_rank(MPI_COMM_WORLD, &tmppid); \
+ if (!tmppid) { \
+ printf(__VA_ARGS__); \
+ } \
+ }
+
+#define ARRAY_SIZE(arr) (sizeof(arr) / sizeof(*arr))
+
+// NEED TO BE DEFINED HERE. IS NOT NOTICED BY compile_acc call.
+#define LFORCING (0)
+#define LSINK (0)
+
+// Write all setting info into a separate ascii file. This is done to guarantee
+// that we have the data specifi information in the thing, even though in
+// principle these things are in the astaroth.conf.
+static inline void
+write_info(const AcMeshInfo* config)
+{
+
+ FILE* infotxt;
+
+ infotxt = fopen("purge.sh", "w");
+ fprintf(infotxt, "#!/bin/bash\n");
+ fprintf(infotxt, "rm *.list *.mesh *.ts purge.sh\n");
+ fclose(infotxt);
+
+ infotxt = fopen("info.list", "w");
+
+ // Determine endianness
+ unsigned int EE = 1;
+ char* CC = (char*)&EE;
+ const int endianness = (int)*CC;
+ // endianness = 0 -> big endian
+ // endianness = 1 -> little endian
+
+ fprintf(infotxt, "size_t %s %lu \n", "AcRealSize", sizeof(AcReal));
+
+ fprintf(infotxt, "int %s %i \n", "endian", endianness);
+
+ // JP: this could be done shorter and with smaller chance for errors with the following
+ // (modified from acPrintMeshInfo() in astaroth.cu)
+ // MV: Now adapted into working condition. E.g. removed useless / harmful formatting.
+
+ for (int i = 0; i < NUM_INT_PARAMS; ++i)
+ fprintf(infotxt, "int %s %d\n", intparam_names[i], config->int_params[i]);
+
+ for (int i = 0; i < NUM_INT3_PARAMS; ++i)
+ fprintf(infotxt, "int3 %s %d %d %d\n", int3param_names[i], config->int3_params[i].x,
+ config->int3_params[i].y, config->int3_params[i].z);
+
+ for (int i = 0; i < NUM_REAL_PARAMS; ++i)
+ fprintf(infotxt, "real %s %g\n", realparam_names[i], double(config->real_params[i]));
+
+ for (int i = 0; i < NUM_REAL3_PARAMS; ++i)
+ fprintf(infotxt, "real3 %s %g %g %g\n", real3param_names[i],
+ double(config->real3_params[i].x), double(config->real3_params[i].y),
+ double(config->real3_params[i].z));
+
+ fclose(infotxt);
+}
+
+// This funtion writes a run state into a set of C binaries.
+static inline void
+save_mesh(const AcMesh& save_mesh, const int step, const AcReal t_step)
+{
+ FILE* save_ptr;
+
+ for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
+ const size_t n = acVertexBufferSize(save_mesh.info);
+
+ const char* buffername = vtxbuf_names[w];
+ char cstep[11];
+ char bin_filename[80] = "\0";
+
+ // sprintf(bin_filename, "");
+
+ sprintf(cstep, "%d", step);
+
+ strcat(bin_filename, buffername);
+ strcat(bin_filename, "_");
+ strcat(bin_filename, cstep);
+ strcat(bin_filename, ".mesh");
+
+ printf("Savefile %s \n", bin_filename);
+
+ save_ptr = fopen(bin_filename, "wb");
+
+ // Start file with time stamp
+ AcReal write_long_buf = (AcReal)t_step;
+ fwrite(&write_long_buf, sizeof(AcReal), 1, save_ptr);
+ // Grid data
+ for (size_t i = 0; i < n; ++i) {
+ const AcReal point_val = save_mesh.vertex_buffer[VertexBufferHandle(w)][i];
+ AcReal write_long_buf2 = (AcReal)point_val;
+ fwrite(&write_long_buf2, sizeof(AcReal), 1, save_ptr);
+ }
+ fclose(save_ptr);
+ }
+}
+
+// This funtion reads a run state from a set of C binaries.
+static inline void
+read_mesh(AcMesh& read_mesh, const int step, AcReal* t_step)
+{
+ FILE* read_ptr;
+
+ for (int w = 0; w < NUM_VTXBUF_HANDLES; ++w) {
+ const size_t n = acVertexBufferSize(read_mesh.info);
+
+ const char* buffername = vtxbuf_names[w];
+ char cstep[11];
+ char bin_filename[80] = "\0";
+
+ // sprintf(bin_filename, "");
+
+ sprintf(cstep, "%d", step);
+
+ strcat(bin_filename, buffername);
+ strcat(bin_filename, "_");
+ strcat(bin_filename, cstep);
+ strcat(bin_filename, ".mesh");
+
+ printf("Reading savefile %s \n", bin_filename);
+
+ read_ptr = fopen(bin_filename, "rb");
+
+ // Start file with time stamp
+ size_t result;
+ result = fread(t_step, sizeof(AcReal), 1, read_ptr);
+ // Read grid data
+ AcReal read_buf;
+ for (size_t i = 0; i < n; ++i) {
+ result = fread(&read_buf, sizeof(AcReal), 1, read_ptr);
+ read_mesh.vertex_buffer[VertexBufferHandle(w)][i] = read_buf;
+ if (int(result) != 1) {
+ fprintf(stderr, "Reading error in %s, element %i\n", vtxbuf_names[w], int(i));
+ fprintf(stderr, "Result = %i, \n", int(result));
+ }
+ }
+ fclose(read_ptr);
+ }
+}
+
+// This function prints out the diagnostic values to std.out and also saves and
+// appends an ascii file to contain all the result.
+// JP: MUST BE CALLED FROM PROC 0. Must be rewritten for multiple processes (this implementation
+// has write race condition)
+static inline void
+print_diagnostics_host(const AcMesh mesh, const int step, const AcReal dt, const AcReal t_step,
+ FILE* diag_file, const AcReal sink_mass, const AcReal accreted_mass)
+{
+
+ AcReal buf_rms, buf_max, buf_min;
+ const int max_name_width = 16;
+
+ // Calculate rms, min and max from the velocity vector field
+ buf_max = acModelReduceVec(mesh, RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
+ buf_min = acModelReduceVec(mesh, RTYPE_MIN, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
+ buf_rms = acModelReduceVec(mesh, RTYPE_RMS, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ);
+
+ // MV: The ordering in the earlier version was wrong in terms of variable
+ // MV: name and its diagnostics.
+ printf("Step %d, t_step %.3e, dt %e s\n", step, double(t_step), double(dt));
+ printf(" %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, "uu total", double(buf_min),
+ double(buf_rms), double(buf_max));
+ fprintf(diag_file, "%d %e %e %e %e %e ", step, double(t_step), double(dt), double(buf_min),
+ double(buf_rms), double(buf_max));
+
+ // Calculate rms, min and max from the variables as scalars
+ for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+ buf_max = acModelReduceScal(mesh, RTYPE_MAX, VertexBufferHandle(i));
+ buf_min = acModelReduceScal(mesh, RTYPE_MIN, VertexBufferHandle(i));
+ buf_rms = acModelReduceScal(mesh, RTYPE_RMS, VertexBufferHandle(i));
+
+ printf(" %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, vtxbuf_names[i],
+ double(buf_min), double(buf_rms), double(buf_max));
+ fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max));
+ }
+
+ if ((sink_mass >= AcReal(0.0)) || (accreted_mass >= AcReal(0.0))) {
+ fprintf(diag_file, "%e %e ", double(sink_mass), double(accreted_mass));
+ }
+
+ fprintf(diag_file, "\n");
+}
+
+// This function prints out the diagnostic values to std.out and also saves and
+// appends an ascii file to contain all the result.
+// JP: EXECUTES ON MULTIPLE GPUS, MUST BE CALLED FROM ALL PROCS
+static inline void
+print_diagnostics(const int step, const AcReal dt, const AcReal t_step, FILE* diag_file,
+ const AcReal sink_mass, const AcReal accreted_mass)
+{
+
+ AcReal buf_rms, buf_max, buf_min;
+ const int max_name_width = 16;
+
+ // Calculate rms, min and max from the velocity vector field
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &buf_max);
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_MIN, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &buf_min);
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_RMS, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &buf_rms);
+
+ // MV: The ordering in the earlier version was wrong in terms of variable
+ // MV: name and its diagnostics.
+ printf("Step %d, t_step %.3e, dt %e s\n", step, double(t_step), double(dt));
+ printf(" %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, "uu total", double(buf_min),
+ double(buf_rms), double(buf_max));
+ fprintf(diag_file, "%d %e %e %e %e %e ", step, double(t_step), double(dt), double(buf_min),
+ double(buf_rms), double(buf_max));
+
+ // Calculate rms, min and max from the variables as scalars
+ for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+ acGridReduceScal(STREAM_DEFAULT, RTYPE_MAX, VertexBufferHandle(i), &buf_max);
+ acGridReduceScal(STREAM_DEFAULT, RTYPE_MIN, VertexBufferHandle(i), &buf_min);
+ acGridReduceScal(STREAM_DEFAULT, RTYPE_RMS, VertexBufferHandle(i), &buf_rms);
+
+ printf(" %*s: min %.3e,\trms %.3e,\tmax %.3e\n", max_name_width, vtxbuf_names[i],
+ double(buf_min), double(buf_rms), double(buf_max));
+ fprintf(diag_file, "%e %e %e ", double(buf_min), double(buf_rms), double(buf_max));
+ }
+
+ if ((sink_mass >= AcReal(0.0)) || (accreted_mass >= AcReal(0.0))) {
+ fprintf(diag_file, "%e %e ", double(sink_mass), double(accreted_mass));
+ }
+
+ fprintf(diag_file, "\n");
+}
+
+/*
+ MV NOTE: At the moment I have no clear idea how to calculate magnetic
+ diagnostic variables from grid. Vector potential measures have a limited
+ value. TODO: Smart way to get brms, bmin and bmax.
+*/
+
+#include "math_utils.h"
+AcReal
+calc_timestep(const AcMeshInfo info)
+{
+ AcReal uumax;
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &uumax);
+
+ const long double cdt = info.real_params[AC_cdt];
+ const long double cdtv = info.real_params[AC_cdtv];
+ // const long double cdts = info.real_params[AC_cdts];
+ const long double cs2_sound = info.real_params[AC_cs2_sound];
+ const long double nu_visc = info.real_params[AC_nu_visc];
+ const long double eta = info.real_params[AC_eta];
+ const long double chi = 0; // info.real_params[AC_chi]; // TODO not calculated
+ const long double gamma = info.real_params[AC_gamma];
+ const long double dsmin = info.real_params[AC_dsmin];
+
+ // Old ones from legacy Astaroth
+ // const long double uu_dt = cdt * (dsmin / (uumax + cs_sound));
+ // const long double visc_dt = cdtv * dsmin * dsmin / nu_visc;
+
+ // New, closer to the actual Courant timestep
+ // See Pencil Code user manual p. 38 (timestep section)
+ const long double uu_dt = cdt * dsmin / (fabsl(uumax) + sqrtl(cs2_sound + 0.0l));
+ const long double visc_dt = cdtv * dsmin * dsmin / max(max(nu_visc, eta), max(gamma, chi));
+
+ const long double dt = min(uu_dt, visc_dt);
+ ERRCHK_ALWAYS(is_valid((AcReal)dt));
+ return AcReal(dt);
+}
+
+int
+main(int argc, char** argv)
+{
+ MPI_Init(NULL, NULL);
+ int nprocs, pid;
+ MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+ MPI_Comm_rank(MPI_COMM_WORLD, &pid);
+
+ /////////// Simple example START
+ (void)argc; // Unused
+ (void)argv; // Unused
+
+ // Set random seed for reproducibility
+ srand(321654987);
+
+ AcMeshInfo info;
+ acLoadConfig(AC_DEFAULT_CONFIG, &info);
+ load_config(AC_DEFAULT_CONFIG, &info);
+ update_config(&info);
+
+ AcMesh mesh;
+ if (pid == 0) {
+ acMeshCreate(info, &mesh);
+ acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, &mesh);
+ }
+ acGridInit(info);
+ acGridLoadMesh(STREAM_DEFAULT, mesh);
+
+ for (int t_step = 0; t_step < 100; ++t_step) {
+ const AcReal dt = calc_timestep(info);
+ acGridIntegrate(STREAM_DEFAULT, dt);
+
+ if (1) { // Diag step
+ AcReal uumin, uumax, uurms;
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_MIN, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &uumin);
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &uumax);
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_RMS, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &uurms);
+
+ printf("Step %d, dt: %g\n", t_step, (double)dt);
+ printf("%*s: %.3e, %.3e, %.3e\n", 16, "UU", (double)uumin, (double)uumax,
+ (double)uurms);
+
+ for (size_t vtxbuf = 0; vtxbuf < NUM_VTXBUF_HANDLES; ++vtxbuf) {
+ AcReal scalmin, scalmax, scalrms;
+ acGridReduceScal(STREAM_DEFAULT, RTYPE_MIN, (VertexBufferHandle)vtxbuf, &scalmin);
+ acGridReduceScal(STREAM_DEFAULT, RTYPE_MAX, (VertexBufferHandle)vtxbuf, &scalmax);
+ acGridReduceScal(STREAM_DEFAULT, RTYPE_RMS, (VertexBufferHandle)vtxbuf, &scalrms);
+
+ printf("%*s: %.3e, %.3e, %.3e\n", 16, vtxbuf_names[vtxbuf], (double)scalmin,
+ (double)scalmax, (double)scalrms);
+ }
+ }
+ }
+ if (pid == 0)
+ acMeshDestroy(&mesh);
+
+ acGridQuit();
+ /////////////// Simple example END
+
+// JP: The following is directly from standalone/simulation.cc and modified to work with MPI
+// However, not extensively tested
+#if 0
+ // Load config to AcMeshInfo
+ AcMeshInfo info;
+ if (argc > 1) {
+ if (argc == 3 && (!strcmp(argv[1], "-c") || !strcmp(argv[1], "--config"))) {
+ acLoadConfig(argv[2], &info);
+ load_config(argv[2], &info);
+ acUpdateBuiltinParams(&info);
+ }
+ else {
+ printf("Usage: ./ac_run\n");
+ printf("Usage: ./ac_run -c \n");
+ printf("Usage: ./ac_run --config \n");
+ return EXIT_FAILURE;
+ }
+ }
+ else {
+ acLoadConfig(AC_DEFAULT_CONFIG, &info);
+ load_config(AC_DEFAULT_CONFIG, &info);
+ acUpdateBuiltinParams(&info);
+ }
+
+ const int start_step = info.int_params[AC_start_step];
+ const int max_steps = info.int_params[AC_max_steps];
+ const int save_steps = info.int_params[AC_save_steps];
+ const int bin_save_steps = info.int_params[AC_bin_steps];
+
+ const AcReal max_time = info.real_params[AC_max_time];
+ const AcReal bin_save_t = info.real_params[AC_bin_save_t];
+ AcReal bin_crit_t = bin_save_t;
+ AcReal t_step = 0.0;
+ FILE* diag_file = fopen("timeseries.ts", "a");
+ ERRCHK_ALWAYS(diag_file);
+
+ AcMesh mesh;
+ ///////////////////////////////// PROC 0 BLOCK START ///////////////////////////////////////////
+ if (pid == 0) {
+ acMeshCreate(info, &mesh);
+ // TODO: This need to be possible to define in astaroth.conf
+ acmesh_init_to(INIT_TYPE_GAUSSIAN_RADIAL_EXPL, &mesh);
+ // acmesh_init_to(INIT_TYPE_SIMPLE_CORE, mesh); //Initial condition for a collapse test
+
+#if LSINK
+ acVertexBufferSet(VTXBUF_ACCRETION, 0.0, &mesh);
+#endif
+ // Read old binary if we want to continue from an existing snapshot
+ // WARNING: Explicit specification of step needed!
+ if (start_step > 0) {
+ read_mesh(mesh, start_step, &t_step);
+ }
+
+ // Generate the title row.
+ if (start_step == 0) {
+ fprintf(diag_file, "step t_step dt uu_total_min uu_total_rms uu_total_max ");
+ for (int i = 0; i < NUM_VTXBUF_HANDLES; ++i) {
+ fprintf(diag_file, "%s_min %s_rms %s_max ", vtxbuf_names[i], vtxbuf_names[i],
+ vtxbuf_names[i]);
+ }
+ }
+
+#if LSINK
+ fprintf(diag_file, "sink_mass accreted_mass ");
+#endif
+ fprintf(diag_file, "\n");
+
+ write_info(&info);
+
+ if (start_step == 0) {
+#if LSINK
+ print_diagnostics_host(mesh, 0, AcReal(.0), t_step, diag_file,
+ info.real_params[AC_M_sink_init], 0.0);
+#else
+ print_diagnostics_host(mesh, 0, AcReal(.0), t_step, diag_file, -1.0, -1.0);
+#endif
+ }
+
+ acMeshApplyPeriodicBounds(&mesh);
+ if (start_step == 0) {
+ save_mesh(mesh, 0, t_step);
+ }
+ }
+ ////////////////////////////////// PROC 0 BLOCK END ////////////////////////////////////////////
+
+ // Init GPU
+ acGridInit(info);
+ acGridLoadMesh(STREAM_DEFAULT, mesh);
+
+ /* initialize random seed: */
+ srand(312256655);
+
+ /* Step the simulation */
+ AcReal accreted_mass = 0.0;
+ AcReal sink_mass = 0.0;
+ for (int i = start_step + 1; i < max_steps; ++i) {
+ AcReal umax;
+ acGridReduceVec(STREAM_DEFAULT, RTYPE_MAX, VTXBUF_UUX, VTXBUF_UUY, VTXBUF_UUZ, &umax);
+ const AcReal dt = host_timestep(umax, info);
+
+#if LSINK
+ AcReal sum_mass;
+ acGridReduceScal(STREAM_DEFAULT, RTYPE_SUM, VTXBUF_ACCRETION, &sum_mass);
+ accreted_mass = accreted_mass + sum_mass;
+ sink_mass = 0.0;
+ sink_mass = info.real_params[AC_M_sink_init] + accreted_mass;
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_M_sink, sink_mass);
+
+ // JP: !!! WARNING !!! acVertexBufferSet operates in host memory. The mesh is
+ // never loaded to device memory. Is this intended?
+ if (pid == 0)
+ acVertexBufferSet(VTXBUF_ACCRETION, 0.0, mesh);
+
+ int on_off_switch;
+ if (i < 1) {
+ on_off_switch = 0; // accretion is off till certain amount of steps.
+ }
+ else {
+ on_off_switch = 1;
+ }
+ acGridLoadScalarUniform(STREAM_DEFAULT, AC_switch_accretion, on_off_switch);
+#else
+ accreted_mass = -1.0;
+ sink_mass = -1.0;
+#endif
+
+#if LFORCING
+ const ForcingParams forcing_params = generateForcingParams(info);
+ loadForcingParamsToGrid(forcing_params);
+#endif
+
+ acGridIntegrate(STREAM_DEFAULT, dt);
+
+ t_step += dt;
+
+ /* Save the simulation state and print diagnostics */
+ if ((i % save_steps) == 0) {
+
+ /*
+ print_diagnostics() writes out both std.out printout from the
+ results and saves the diagnostics into a table for ascii file
+ timeseries.ts.
+ */
+
+ print_diagnostics(i, dt, t_step, diag_file, sink_mass, accreted_mass);
+#if LSINK
+ printf("sink mass is: %.15e \n", double(sink_mass));
+ printf("accreted mass is: %.15e \n", double(accreted_mass));
+#endif
+ /*
+ We would also might want an XY-average calculating funtion,
+ which can be very useful when observing behaviour of turbulent
+ simulations. (TODO)
+ */
+ }
+
+ /* Save the simulation state and print diagnostics */
+ if ((i % bin_save_steps) == 0 || t_step >= bin_crit_t) {
+
+ /*
+ This loop saves the data into simple C binaries which can be
+ used for analysing the data snapshots closely.
+
+ The updated mesh will be located on the GPU. Also all calls
+ to the astaroth interface (functions beginning with ac*) are
+ assumed to be asynchronous, so the meshes must be also synchronized
+ before transferring the data to the CPU. Like so:
+
+ acBoundcondStep();
+ acStore(mesh);
+ */
+ acGridPeriodicBoundconds(STREAM_DEFAULT);
+ acGridStoreMesh(STREAM_DEFAULT, &mesh);
+
+ if (pid == 0)
+ save_mesh(mesh, i, t_step);
+
+ bin_crit_t += bin_save_t;
+ }
+
+ // End loop if max time reached.
+ if (max_time > AcReal(0.0)) {
+ if (t_step >= max_time) {
+ printf("Time limit reached! at t = %e \n", double(t_step));
+ break;
+ }
+ }
+ }
+
+ //////Save the final snapshot
+ ////acSynchronize();
+ ////acStore(mesh);
+
+ ////save_mesh(*mesh, , t_step);
+
+ acGridQuit();
+ if (pid == 0)
+ acMeshDestroy(&mesh);
+ fclose(diag_file);
+#endif
+
+ MPI_Finalize();
+ return EXIT_SUCCESS;
+}
+
+#else
+#include
+#include
+
+int
+main(void)
+{
+ printf("The library was built without MPI support, cannot run mpitest. Rebuild Astaroth with "
+ "cmake -DMPI_ENABLED=ON .. to enable.\n");
+ return EXIT_FAILURE;
+}
+#endif // AC_MPI_ENABLES