cwpearson/astaroth
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added a stripped down MPI version of standalone: standalone_mpi. In fact, it's more like a pure simulation module since I've dropped real-time rendering and other old parts that do not work with MPI without heavy modifications. The most important functionalities in addition to simulation have already been adapted to work with MPI (samples/benchmark and samples/mpi) so there's no need to re-create them in standalone_mpi. The current version of standalone_mpi is able to run a basic simulation and I get an agreement with non-mpi and mpi versions after 100 timesteps. There's also draft that's a direct adaptation of what's currently in standalone/simulation.cc (it should be 100% equivalent), but it's currently commented out as I haven't done extensive tests with it.

Browse Source
This commit is contained in:
jpekkila
2020-08-24 19:03:03 +03:00
parent f21c6a8c0b
commit cec9a23dc0
9 changed files with 1928 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@@ -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)

2
samples/standalone_mpi/CMakeLists.txt Normal file
View File

@@ -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)

184
samples/standalone_mpi/config_loader.cc Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @file
* \brief Brief info.
*
* Detailed info.
*
*/
#include "config_loader.h"
#include <limits.h> // UINT_MAX
#include <stdint.h> // uint8_t, uint32_t
#include <stdio.h> // print
#include <string.h> // 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;
}

34
samples/standalone_mpi/config_loader.h Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @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);

294
samples/standalone_mpi/host_forcing.cc Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @file
* \brief Brief info.
*
* Detailed info.
*
*/
#include "host_forcing.h"
// #include "math_utils.h"
#include <cmath>
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(&params.ff_hel_re, &params.ff_hel_im, params.k_force, e_force,
relhel);
return params;
}

60
samples/standalone_mpi/host_forcing.h Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @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);

718
samples/standalone_mpi/host_memory.cc Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @file
* \brief Brief info.
*
* Detailed info.
*
*/
#include "host_memory.h"
#include <math.h>
#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;
}
}
*/
}

49
samples/standalone_mpi/host_memory.h Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @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);

586
samples/standalone_mpi/main.cc Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
Running: mpirun -np <num processes> <executable>
*/
#if AC_MPI_ENABLED
#include "astaroth.h"
#include "astaroth_utils.h"
#include <mpi.h>
#include <string.h>
#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 <config_path>\n");
printf("Usage: ./ac_run --config <config_path>\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 <stdio.h>
#include <stdlib.h>
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