commit 5d95f6fc19e0c30bd68f0142558ec2dce788655f
parent e9dcf9319aa389468cbf24ceeeb396b11b1b778a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Tue, 23 Mar 2021 14:44:24 +0100
Add the atrstm_fetch_rdgfa function
Diffstat:
15 files changed, 441 insertions(+), 83 deletions(-)
diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt
@@ -69,6 +69,7 @@ set(ATRSTM_FILES_SRC
atrstm_log.c
atrstm_partition.c
atrstm_radcoefs.c
+ atrstm_rdgfa.c
atrstm_setup_octrees.c
atrstm_setup_uvm.c
atrstm_svx.c)
@@ -90,7 +91,9 @@ set(ATRSTM_FILES_DOC COPYING README.md)
# SIMD files
if(USE_SIMD)
set(ATRSTM_FILES_SRC ${ATRSTM_FILES_SRC} atrstm_radcoefs_simd4.c)
- set(ATRSTM_FILES_INC ${ATRSTM_FILES_INC} atrstm_radcoefs_simd4.h)
+ set(ATRSTM_FILES_INC ${ATRSTM_FILES_INC}
+ atrstm_radcoefs_simd4.h
+ atrstm_rdgfa_simd4.h)
endif()
# Prepend each file in the `ATRSTM_FILES_<SRC|INC>' list by `ATRSTM_SOURCE_DIR'
diff --git a/src/atrstm.c b/src/atrstm.c
@@ -44,7 +44,7 @@ check_args(const struct atrstm_args* args)
|| !args->name
|| (unsigned)args->spectral_type >= ATRSTM_SPECTRAL_TYPES_COUNT__
|| args->wlen_range[0] > args->wlen_range[1]
- || args->gyration_radius_prefactor <= 0
+ || args->fractal_prefactor <= 0
|| args->fractal_dimension <= 0
|| args->optical_thickness < 0
|| !args->nthreads) {
@@ -180,7 +180,7 @@ atrstm_create
atrstm->allocator = allocator;
atrstm->verbose = args->verbose;
atrstm->nthreads = MMIN(args->nthreads, (unsigned)nthreads_max);
- atrstm->gyration_radius_prefactor = args->gyration_radius_prefactor;
+ atrstm->fractal_prefactor = args->fractal_prefactor;
atrstm->fractal_dimension = args->fractal_dimension;
atrstm->spectral_type = args->spectral_type;
atrstm->wlen_range[0] = args->wlen_range[0];
@@ -317,6 +317,43 @@ atrstm_get_name(const struct atrstm* atrstm)
return str_cget(&atrstm->name);
}
+res_T
+atrstm_at
+ (const struct atrstm* atrstm,
+ const double pos[3],
+ struct suvm_primitive* prim, /* Volumetric primitive */
+ double bcoords[4]) /* `pos' in `prim' */
+{
+ res_T res = RES_OK;
+
+ if(!atrstm || !pos || !prim || !bcoords) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Find the primitive into which pos lies */
+ res = suvm_volume_at(atrstm->volume, pos, prim, bcoords);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Error accessing the volumetric mesh at {%g, %g, %g} -- %s.\n",
+ SPLIT3(pos), res_to_cstr(res));
+ goto error;
+ }
+
+ /* No primitive found */
+ if(SUVM_PRIMITIVE_NONE(prim)) {
+ log_warn(atrstm,
+ "%s: the position {%g, %g, %g} does not belong to the semi-transparent "
+ "medium.\n", FUNC_NAME, SPLIT3(pos));
+ goto exit;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
/*******************************************************************************
* Local functions
******************************************************************************/
diff --git a/src/atrstm.h b/src/atrstm.h
@@ -16,6 +16,7 @@
#ifndef ATRSTM_H
#define ATRSTM_H
+#include <star/suvm.h>
#include <star/svx.h>
#include <rsys/rsys.h>
@@ -101,7 +102,7 @@ struct atrstm_args {
enum atrstm_spectral_type spectral_type; /* Longwave/shortwave */
double wlen_range[2]; /* Spectral range to handle In nm */
- double gyration_radius_prefactor;
+ double fractal_prefactor;
double fractal_dimension;
unsigned grid_max_definition[3]; /* Fixed grid definition along the 3 axes */
@@ -127,7 +128,7 @@ struct atrstm_args {
ATRSTM_SPECTRAL_SW, /* Spectral type */ \
{500,500}, /* Spectral integration range */ \
\
- 1.70, /* Gyration radius prefactor */ \
+ 1.70, /* Fractal prefactor */ \
1.75, /* Fractal dimension */ \
\
{256, 256, 256}, /* Acceleration grid max definition */ \
@@ -144,7 +145,8 @@ struct atrstm_args {
static const struct atrstm_args ATRSTM_ARGS_DEFAULT = ATRSTM_ARGS_DEFAULT__;
struct atrstm_fetch_radcoefs_args {
- double pos[3]; /* Position to query */
+ struct suvm_primitive prim; /* Volumetric primitive to query */
+ double bcoords[4]; /* Barycentric coordinates of the pos in prim to query */
size_t iband; /* Spectral band index. Not use in shortwave */
size_t iquad; /* Quadrature point index. Not use in shortwave */
double wavelength; /* In nm */
@@ -158,7 +160,8 @@ struct atrstm_fetch_radcoefs_args {
};
#define ATRSTM_FETCH_RADCOEFS_ARGS_DEFAULT__ { \
- {0, 0, 0}, /* Position */ \
+ SUVM_PRIMITIVE_NULL__, \
+ {0, 0, 0, 1}, /* Barycentric coordinates */ \
SIZE_MAX, /* Index of the spectral band */ \
SIZE_MAX, /* Index of the quadrature point */ \
DBL_MAX, /* Wavelength */ \
@@ -256,6 +259,30 @@ struct atrstm_trace_ray_args {
static const struct atrstm_trace_ray_args ATRSTM_TRACE_RAY_ARGS_DEFAULT =
ATRSTM_TRACE_RAY_ARGS_DEFAULT__;
+struct atrstm_fetch_rdgfa_args {
+ struct suvm_primitive prim; /* Volumetric primitive to query */
+ double bcoords[4]; /* Barycentric coordinates of the pos in prim to query */
+ double wavelength; /* In nm */
+};
+
+#define ATRSTM_FETCH_RDGFA_ARGS_DEFAULT__ { \
+ SUVM_PRIMITIVE_NULL__, \
+ {0, 0, 0, 1}, /* Barycentric coordinates */ \
+ DBL_MAX, /* Wavelength */ \
+}
+static const struct atrstm_fetch_rdgfa_args
+ATRSTM_FETCH_RDGFA_ARGS_DEFAULT = ATRSTM_FETCH_RDGFA_ARGS_DEFAULT__;
+
+struct atrstm_rdgfa {
+ double wavelength; /* In nm */
+ double fractal_dimension;
+ double gyration_radius;
+};
+
+#define ATRSTM_RDGFA_NULL__ {0,0,0}
+static const struct atrstm_rdgfa ATRSTM_RDGFA_NULL =
+ ATRSTM_RDGFA_NULL__;
+
/* Types used as syntactic sugar that store the radiative coefficients */
typedef double atrstm_radcoefs_T[ATRSTM_RADCOEFS_COUNT__];
typedef double atrstm_radcoefs_svx_T[ATRSTM_RADCOEFS_COUNT__][ATRSTM_SVX_OPS_COUNT__];
@@ -292,6 +319,13 @@ atrstm_get_name
(const struct atrstm* atrstm);
ATRSTM_API res_T
+atrstm_at
+ (const struct atrstm* atrstm,
+ const double pos[3],
+ struct suvm_primitive* prim, /* Volumetric primitive */
+ double barycentric_coords[4]); /* `pos' in `prim' */
+
+ATRSTM_API res_T
atrstm_fetch_radcoefs
(const struct atrstm* atrstm,
const struct atrstm_fetch_radcoefs_args* args,
@@ -322,6 +356,12 @@ atrstm_dump_svx_octree
const struct atrstm_dump_svx_octree_args* args,
FILE* stream);
+ATRSTM_API res_T
+atrstm_fetch_rdgfa
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_rdgfa_args* args,
+ struct atrstm_rdgfa* rdgfa);
+
END_DECLS
#endif /* ATRSTM_H */
diff --git a/src/atrstm_c.h b/src/atrstm_c.h
@@ -46,7 +46,7 @@ struct atrstm {
struct str name; /* Name of the semi-transparent medium */
- double gyration_radius_prefactor;
+ double fractal_prefactor;
double fractal_dimension;
enum atrstm_spectral_type spectral_type;
diff --git a/src/atrstm_cache.c b/src/atrstm_cache.c
@@ -161,7 +161,7 @@ write_cache_header(struct cache* cache)
} \
} (void)0
WRITE(&CACHE_VERSION, 1);
- WRITE(&cache->atrstm->gyration_radius_prefactor, 1);
+ WRITE(&cache->atrstm->fractal_prefactor, 1);
WRITE(&cache->atrstm->fractal_dimension, 1);
WRITE(&cache->atrstm->spectral_type, 1);
WRITE(cache->atrstm->wlen_range, 2);
@@ -185,7 +185,7 @@ read_cache_header(struct cache* cache)
{
struct hash cached_hash;
struct hash current_hash;
- double gyration_radius_prefactor = 0;
+ double fractal_prefactor = 0;
double fractal_dimension = 0;
enum atrstm_spectral_type spectral_type = ATRSTM_SPECTRAL_TYPES_COUNT__;
double wlen_range[2] = {0,0};
@@ -230,9 +230,9 @@ read_cache_header(struct cache* cache)
} \
} (void)0
- READ(&gyration_radius_prefactor, 1);
- CHK_VAR(gyration_radius_prefactor, cache->atrstm->gyration_radius_prefactor,
- "gyration radius prefactor", "%g");
+ READ(&fractal_prefactor, 1);
+ CHK_VAR(fractal_prefactor, cache->atrstm->fractal_prefactor,
+ "fractal prefactor", "%g");
READ(&fractal_dimension, 1);
CHK_VAR(fractal_dimension, cache->atrstm->fractal_dimension,
diff --git a/src/atrstm_radcoefs.c b/src/atrstm_radcoefs.c
@@ -17,6 +17,10 @@
#include "atrstm_log.h"
#include "atrstm_radcoefs.h"
+#ifdef ATRSTM_USE_SIMD
+ #include "atrstm_radcoefs_simd4.h"
+#endif
+
#include <astoria/atrri.h>
#include <astoria/atrtp.h>
@@ -27,9 +31,51 @@
#include <string.h>
/*******************************************************************************
- * Helper function
+ * Exported functions
+ ******************************************************************************/
+#ifndef ATRSTM_USE_SIMD
+res_T
+atrstm_fetch_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ atrstm_radcoefs_T radcoefs)
+{
+ return fetch_radcoefs(atrstm, args, radcoefs);
+}
+
+#else /* ATRSTM_USE_SIMD */
+res_T
+atrstm_fetch_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ atrstm_radcoefs_T radcoefs)
+{
+ const res_T res = fetch_radcoefs_simd4(atrstm, args, radcoefs);
+ if(res != RES_OK) goto error;
+
+ #ifndef NDEBUG
+ {
+ atrstm_radcoefs_T radcoefs_scalar;
+ int i;
+
+ ASSERT(fetch_radcoefs(atrstm, args, radcoefs_scalar) == RES_OK);
+ FOR_EACH(i, 0, ATRSTM_RADCOEFS_COUNT__) {
+ ASSERT(eq_eps(radcoefs[i], radcoefs_scalar[i], radcoefs[i]*1e-6));
+ }
+ }
+ #endif
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+#endif /* !ATRSTM_USE_SIMD */
+
+/*******************************************************************************
+ * Local functions functions
******************************************************************************/
-static INLINE int
+int
check_fetch_radcoefs_args
(const struct atrstm* atrstm,
const struct atrstm_fetch_radcoefs_args* args,
@@ -39,6 +85,7 @@ check_fetch_radcoefs_args
ASSERT(atrstm && args && func_name);
if(!args
+ || SUVM_PRIMITIVE_NONE(&args->prim)
|| args->wavelength < atrstm->wlen_range[0]
|| args->wavelength > atrstm->wlen_range[1]
|| !(args->radcoefs_mask & ATRSTM_RADCOEFS_MASK_ALL)
@@ -60,18 +107,13 @@ check_fetch_radcoefs_args
return 1;
}
-/*******************************************************************************
- * Exported functions
- ******************************************************************************/
res_T
-atrstm_fetch_radcoefs
+fetch_radcoefs
(const struct atrstm* atrstm,
const struct atrstm_fetch_radcoefs_args* args,
- double radcoefs[ATRSTM_RADCOEFS_COUNT__])
+ atrstm_radcoefs_T radcoefs)
{
struct radcoefs prim_radcoefs[4];
- double bcoords[4];
- struct suvm_primitive prim = SUVM_PRIMITIVE_NULL;
res_T res = RES_OK;
if(!atrstm
@@ -82,30 +124,13 @@ atrstm_fetch_radcoefs
}
memset(radcoefs, 0, sizeof(double[ATRSTM_RADCOEFS_COUNT__]));
- /* Find the primitive into which pos lies */
- res = suvm_volume_at(atrstm->volume, args->pos, &prim, bcoords);
- if(res != RES_OK) {
- log_err(atrstm,
- "Error accessing the volumetric mesh at {%g, %g, %g} -- %s.\n",
- SPLIT3(args->pos), res_to_cstr(res));
- goto error;
- }
-
- /* No primitive found */
- if(SUVM_PRIMITIVE_NONE(&prim)) {
- log_warn(atrstm,
- "%s: the position {%g, %g, %g} does not belong to the semi-transparent "
- "medium.\n", FUNC_NAME, SPLIT3(args->pos));
- goto exit;
- }
-
/* Compute the radiative properties of the primitive nodes */
- res = primitive_compute_radcoefs
- (atrstm, &atrstm->refract_id, &prim, args->radcoefs_mask, prim_radcoefs);
+ res = primitive_compute_radcoefs(atrstm, &atrstm->refract_id, &args->prim,
+ args->radcoefs_mask, prim_radcoefs);
if(res != RES_OK) {
log_err(atrstm,
"Error computing the radiative properties for the primitive %lu -- %s.\n",
- (unsigned long)prim.iprim,
+ (unsigned long)args->prim.iprim,
res_to_cstr(res));
goto error;
}
@@ -113,28 +138,28 @@ atrstm_fetch_radcoefs
/* Linearly interpolate the node radiative properties */
if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) {
radcoefs[ATRSTM_RADCOEF_Ka] =
- prim_radcoefs[0].ka * bcoords[0]
- + prim_radcoefs[1].ka * bcoords[1]
- + prim_radcoefs[2].ka * bcoords[2]
- + prim_radcoefs[3].ka * bcoords[3];
+ prim_radcoefs[0].ka * args->bcoords[0]
+ + prim_radcoefs[1].ka * args->bcoords[1]
+ + prim_radcoefs[2].ka * args->bcoords[2]
+ + prim_radcoefs[3].ka * args->bcoords[3];
ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] >= args->k_min[ATRSTM_RADCOEF_Ka]);
ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] <= args->k_max[ATRSTM_RADCOEF_Ka]);
}
if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) {
radcoefs[ATRSTM_RADCOEF_Ks] =
- prim_radcoefs[0].ks * bcoords[0]
- + prim_radcoefs[1].ks * bcoords[1]
- + prim_radcoefs[2].ks * bcoords[2]
- + prim_radcoefs[3].ks * bcoords[3];
+ prim_radcoefs[0].ks * args->bcoords[0]
+ + prim_radcoefs[1].ks * args->bcoords[1]
+ + prim_radcoefs[2].ks * args->bcoords[2]
+ + prim_radcoefs[3].ks * args->bcoords[3];
ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] >= args->k_min[ATRSTM_RADCOEF_Ks]);
ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] <= args->k_max[ATRSTM_RADCOEF_Ks]);
}
if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) {
radcoefs[ATRSTM_RADCOEF_Kext] =
- prim_radcoefs[0].kext * bcoords[0]
- + prim_radcoefs[1].kext * bcoords[1]
- + prim_radcoefs[2].kext * bcoords[2]
- + prim_radcoefs[3].kext * bcoords[3];
+ prim_radcoefs[0].kext * args->bcoords[0]
+ + prim_radcoefs[1].kext * args->bcoords[1]
+ + prim_radcoefs[2].kext * args->bcoords[2]
+ + prim_radcoefs[3].kext * args->bcoords[3];
ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] >= args->k_min[ATRSTM_RADCOEF_Kext]);
ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] <= args->k_max[ATRSTM_RADCOEF_Kext]);
}
@@ -145,9 +170,6 @@ error:
goto exit;
}
-/*******************************************************************************
- * Local functions
- ******************************************************************************/
res_T
primitive_compute_radcoefs
(const struct atrstm* atrstm,
@@ -170,7 +192,7 @@ primitive_compute_radcoefs
args.lambda = refract_id->wavelength;
args.n = refract_id->n;
args.kappa = refract_id->kappa;
- args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
+ args.fractal_prefactor = atrstm->fractal_prefactor;
args.fractal_dimension = atrstm->fractal_dimension;
args.radcoefs_mask = radcoefs_mask;
@@ -221,7 +243,7 @@ dump_radcoefs
args.lambda = refract_id->wavelength;
args.n = refract_id->n;
args.kappa = refract_id->kappa;
- args.gyration_radius_prefactor = atrstm->gyration_radius_prefactor;
+ args.fractal_prefactor = atrstm->fractal_prefactor;
args.fractal_dimension = atrstm->fractal_dimension;
args.radcoefs_mask = ATRSTM_RADCOEFS_MASK_ALL;
diff --git a/src/atrstm_radcoefs.h b/src/atrstm_radcoefs.h
@@ -17,6 +17,7 @@
#define ATRSTM_RADCOEFS_H
#include "atrstm.h"
+#include "atrstm_rdgfa.h"
#include <rsys/math.h>
#include <rsys/rsys.h>
@@ -34,7 +35,7 @@ struct radcoefs_compute_args {
double lambda; /* wavelength [nm] */
double n; /* Refractive index (real component) */
double kappa; /* Refractive index (imaginary component) */
- double gyration_radius_prefactor;
+ double fractal_prefactor;
double fractal_dimension;
double soot_volumic_fraction; /* In m^3(soot) / m^3 */
double soot_primary_particles_count;
@@ -50,13 +51,17 @@ struct atrri_refractive_index;
struct atrstm;
struct suvm_primitive;
+extern LOCAL_SYM int
+check_fetch_radcoefs_args
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ const char* func_name);
+
extern LOCAL_SYM res_T
-primitive_compute_radcoefs
+fetch_radcoefs
(const struct atrstm* atrstm,
- const struct atrri_refractive_index* refract_id,
- const struct suvm_primitive* prim,
- const int radcoefs_mask,
- struct radcoefs radcoefs[]); /* Per primitive node radiative coefficients */
+ const struct atrstm_fetch_radcoefs_args* args,
+ atrstm_radcoefs_T radcoefs);
/* Write per node radiative coefficients regarding the submitted refracted
* index */
@@ -66,6 +71,14 @@ dump_radcoefs
const struct atrri_refractive_index* refract_id,
FILE* stream);
+extern LOCAL_SYM res_T
+primitive_compute_radcoefs
+ (const struct atrstm* atrstm,
+ const struct atrri_refractive_index* refract_id,
+ const struct suvm_primitive* prim,
+ const int radcoefs_mask,
+ struct radcoefs radcoefs[]); /* Per primitive node radiative coefficients */
+
static INLINE res_T
primitive_compute_radcoefs_range
(const struct atrstm* atrstm,
@@ -121,7 +134,7 @@ radcoefs_compute
const double Fv = args->soot_volumic_fraction; /* [nm^3(soot)/nm] */
const double Np = args->soot_primary_particles_count;
const double Dp = args->soot_primary_particles_diameter; /* [nm] */
- const double kf = args->gyration_radius_prefactor;
+ const double kf = args->fractal_prefactor;
const double Df = args->fractal_dimension;
const int ka = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) != 0;
const int ks = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) != 0;
@@ -165,7 +178,7 @@ radcoefs_compute
const double Fm = tmp2*tmp2 + Em*Em;
/* Gyration radius */
- const double Rg = 0.5 * Dp * pow(Np/kf, 1.0/Df); /* [nm] */
+ const double Rg = compute_gyration_radius(kf, Df, Np, Dp); /* [nm] */
const double Rg2 = Rg*Rg;
const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5);
diff --git a/src/atrstm_radcoefs_simd4.c b/src/atrstm_radcoefs_simd4.c
@@ -14,15 +14,77 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include "atrstm_c.h"
+#include "atrstm_log.h"
+#include "atrstm_radcoefs.h"
#include "atrstm_radcoefs_simd4.h"
#include <astoria/atrtp.h>
#include <star/suvm.h>
+#include <rsys/cstr.h>
+
/*******************************************************************************
* Local functions
******************************************************************************/
res_T
+fetch_radcoefs_simd4
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ atrstm_radcoefs_T radcoefs)
+{
+ struct radcoefs_simd4 prim_radcoefs;
+ v4f_T bcoords;
+ res_T res = RES_OK;
+
+ if(!atrstm
+ || !check_fetch_radcoefs_args(atrstm, args, FUNC_NAME)
+ || !radcoefs) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+ memset(radcoefs, 0, sizeof(double[ATRSTM_RADCOEFS_COUNT__]));
+
+ /* Compute the radiative properties of the primitive nodes */
+ res = primitive_compute_radcoefs_simd4(atrstm, &atrstm->refract_id,
+ &args->prim, args->radcoefs_mask, &prim_radcoefs);
+ if(res != RES_OK) {
+ log_err(atrstm,
+ "Error computing the radiative properties for the primitive %lu -- %s.\n",
+ (unsigned long)args->prim.iprim,
+ res_to_cstr(res));
+ goto error;
+ }
+
+ bcoords = v4f_set
+ ((float)args->bcoords[0],
+ (float)args->bcoords[1],
+ (float)args->bcoords[2],
+ (float)args->bcoords[3]);
+
+ /* Linearly interpolate the node radiative properties */
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) {
+ radcoefs[ATRSTM_RADCOEF_Ka] = v4f_x(v4f_dot(prim_radcoefs.ka, bcoords));
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] >= args->k_min[ATRSTM_RADCOEF_Ka]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ka] <= args->k_max[ATRSTM_RADCOEF_Ka]);
+ }
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) {
+ radcoefs[ATRSTM_RADCOEF_Ks] = v4f_x(v4f_dot(prim_radcoefs.ks, bcoords));
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] >= args->k_min[ATRSTM_RADCOEF_Ks]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Ks] <= args->k_max[ATRSTM_RADCOEF_Ks]);
+ }
+ if(args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Kext) {
+ radcoefs[ATRSTM_RADCOEF_Kext] = v4f_x(v4f_dot(prim_radcoefs.kext, bcoords));
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] >= args->k_min[ATRSTM_RADCOEF_Kext]);
+ ASSERT(radcoefs[ATRSTM_RADCOEF_Kext] <= args->k_max[ATRSTM_RADCOEF_Kext]);
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+res_T
primitive_compute_radcoefs_simd4
(const struct atrstm* atrstm,
const struct atrri_refractive_index* refract_id,
@@ -44,7 +106,7 @@ primitive_compute_radcoefs_simd4
args.lambda = v4f_set1((float)refract_id->wavelength);
args.n = v4f_set1((float)refract_id->n);
args.kappa = v4f_set1((float)refract_id->kappa);
- args.gyration_radius_prefactor = v4f_set1((float)atrstm->gyration_radius_prefactor);
+ args.fractal_prefactor = v4f_set1((float)atrstm->fractal_prefactor);
args.fractal_dimension = v4f_set1((float)atrstm->fractal_dimension);
args.radcoefs_mask = radcoefs_mask;
diff --git a/src/atrstm_radcoefs_simd4.h b/src/atrstm_radcoefs_simd4.h
@@ -18,6 +18,7 @@
#include "atrstm.h"
#include "atrstm_radcoefs.h"
+#include "atrstm_rdgfa_simd4.h"
#include <rsimd/math.h>
#include <rsimd/rsimd.h>
@@ -37,7 +38,7 @@ struct radcoefs_compute_simd4_args {
v4f_T lambda; /* wavelength [nm] */
v4f_T n; /* Refractive index (real component) */
v4f_T kappa; /* Refractive index (imaginary component) */
- v4f_T gyration_radius_prefactor;
+ v4f_T fractal_prefactor;
v4f_T fractal_dimension;
v4f_T soot_volumic_fraction; /* In m^3(soot) / m^3 */
v4f_T soot_primary_particles_count;
@@ -52,6 +53,12 @@ struct atrri_refractive_index;
struct suvm_primitive;
extern LOCAL_SYM res_T
+fetch_radcoefs_simd4
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_radcoefs_args* args,
+ atrstm_radcoefs_T radcoefs);
+
+extern LOCAL_SYM res_T
primitive_compute_radcoefs_simd4
(const struct atrstm* atrstm,
const struct atrri_refractive_index* refract_id,
@@ -112,7 +119,7 @@ radcoefs_compute_simd4
const v4f_T Fv = args->soot_volumic_fraction; /* [nm^3(soot)/nm] */
const v4f_T Np = args->soot_primary_particles_count;
const v4f_T Dp = args->soot_primary_particles_diameter; /* [nm] */
- const v4f_T kf = args->gyration_radius_prefactor;
+ const v4f_T kf = args->fractal_prefactor;
const v4f_T Df = args->fractal_dimension;
const int ka = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ka) != 0;
const int ks = (args->radcoefs_mask & ATRSTM_RADCOEF_FLAG_Ks) != 0;
@@ -172,13 +179,12 @@ radcoefs_compute_simd4
const v4f_T Fm = v4f_madd(tmp3, tmp3, v4f_mul(Em, Em));
/* Gyration radius */
- const v4f_T tmp4 = v4f_pow(v4f_div(Np, kf), v4f_rcp(Df));
- const v4f_T Rg = v4f_mul(v4f_mul(v4f_set1(0.5f), Dp), tmp4); /* [nm] */
+ const v4f_T Rg = compute_gyration_radius_simd4(kf, Df, Np, Dp); /* [nm] */
const v4f_T Rg2 = v4f_mul(Rg, Rg);
const v4f_T four_k2_Rg2 = v4f_mul(v4f_mul(v4f_set1(4), k2), Rg2);
- const v4f_T tmp5 = v4f_div(four_k2_Rg2, v4f_mul(v4f_set1(3), Df));
- const v4f_T tmp6 = v4f_add(v4f_set1(1), tmp5);
- const v4f_T g = v4f_pow(tmp6, v4f_mul(Df, v4f_set1(-0.5f)));
+ const v4f_T tmp4 = v4f_div(four_k2_Rg2, v4f_mul(v4f_set1(3), Df));
+ const v4f_T tmp5 = v4f_add(v4f_set1(1), tmp4);
+ const v4f_T g = v4f_pow(tmp5, v4f_mul(Df, v4f_set1(-0.5f)));
/* Scattering cross section, [nm^3/aggrefate] */
const v4f_T Np2 = v4f_mul(Np, Np);
diff --git a/src/atrstm_rdgfa.c b/src/atrstm_rdgfa.c
@@ -0,0 +1,96 @@
+/* Copyright (C) 2020, 2021 CNRS
+ *
+ * This program 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.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#include "atrstm_c.h"
+#include "atrstm_rdgfa.h"
+
+#include <astoria/atrtp.h>
+
+/*******************************************************************************
+ * Helper function
+ ******************************************************************************/
+static INLINE int
+check_fetch_rdgfa_args
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_rdgfa_args* args)
+{
+ return args
+ && !SUVM_PRIMITIVE_NONE(&args->prim)
+ && args->wavelength >= atrstm->wlen_range[0]
+ && args->wavelength <= atrstm->wlen_range[1];
+}
+
+/*******************************************************************************
+ * Exported functions
+ ******************************************************************************/
+res_T
+atrstm_fetch_rdgfa
+ (const struct atrstm* atrstm,
+ const struct atrstm_fetch_rdgfa_args* args,
+ struct atrstm_rdgfa* rdgfa)
+{
+ struct atrtp_desc atrtp_desc = ATRTP_DESC_NULL;
+ double soot_primary_particles_count[4];
+ double soot_primary_particles_diameter[4];
+ double Np;
+ double Dp;
+ int inode;
+ res_T res = RES_OK;
+
+ if(!atrstm || !check_fetch_rdgfa_args(atrstm, args) || !rdgfa) {
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ res = atrtp_get_desc(atrstm->atrtp, &atrtp_desc);
+ if(res != RES_OK) goto error;
+
+ FOR_EACH(inode, 0, 4) {
+ const double* node;
+
+ /* Fetch the thermodynamic properties of the node */
+ node = atrtp_desc_get_node_properties(&atrtp_desc, args->prim.indices[inode]);
+
+ /* Fetch the required node attributes */
+ soot_primary_particles_count[inode] =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT];
+ soot_primary_particles_diameter[inode] =
+ node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER];
+ }
+
+ /* Linearly interpolate the node attributes */
+ Np =
+ soot_primary_particles_count[0] * args->bcoords[0]
+ + soot_primary_particles_count[1] * args->bcoords[1]
+ + soot_primary_particles_count[2] * args->bcoords[2]
+ + soot_primary_particles_count[3] * args->bcoords[3];
+ Dp =
+ soot_primary_particles_diameter[0] * args->bcoords[0]
+ + soot_primary_particles_diameter[1] * args->bcoords[1]
+ + soot_primary_particles_diameter[2] * args->bcoords[2]
+ + soot_primary_particles_diameter[3] * args->bcoords[3];
+
+ /* Setup output data */
+ rdgfa->wavelength = args->wavelength;
+ rdgfa->fractal_dimension = atrstm->fractal_dimension;
+ rdgfa->gyration_radius = compute_gyration_radius
+ (atrstm->fractal_prefactor, atrstm->fractal_dimension, Np, Dp);
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
diff --git a/src/atrstm_rdgfa.h b/src/atrstm_rdgfa.h
@@ -0,0 +1,38 @@
+/* Copyright (C) 2020, 2021 CNRS
+ *
+ * This program 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.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef ATRSTM_RDGFA_H
+#define ATRSTM_RDGFA_H
+
+#include <rsys/rsys.h>
+#include <math.h>
+
+static FINLINE double /* In nanometer */
+compute_gyration_radius
+ (const double fractal_prefactor,
+ const double fractal_dimension,
+ const double soot_primary_particles_count,
+ const double soot_primary_particles_diameter)
+{
+ const double kf = fractal_prefactor;
+ const double Df = fractal_dimension;
+ const double Np = soot_primary_particles_count;
+ const double Dp = soot_primary_particles_diameter;
+ const double Rg = 0.5 * Dp * pow(Np/kf, 1.0/Df); /* [nm] */
+ return Rg;
+}
+
+#endif /* ATRSTM_RDGFA_H */
+
diff --git a/src/atrstm_rdgfa_simd4.h b/src/atrstm_rdgfa_simd4.h
@@ -0,0 +1,41 @@
+/* Copyright (C) 2020, 2021 CNRS
+ *
+ * This program 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.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>. */
+
+#ifndef ATRSTM_RDGFA_SIMD4_H
+#define ATRSTM_RDGFA_SIMD4_H
+
+#include <rsys/rsys.h>
+
+#include <rsimd/math.h>
+#include <rsimd/rsimd.h>
+
+static FINLINE v4f_T /* In nanometer */
+compute_gyration_radius_simd4
+ (const v4f_T fractal_prefactor,
+ const v4f_T fractal_dimension,
+ const v4f_T soot_primary_particles_count,
+ const v4f_T soot_primary_particles_diameter)
+{
+ const v4f_T kf = fractal_prefactor;
+ const v4f_T Df = fractal_dimension;
+ const v4f_T Np = soot_primary_particles_count;
+ const v4f_T Dp = soot_primary_particles_diameter;
+ const v4f_T tmp = v4f_pow(v4f_div(Np, kf), v4f_rcp(Df));
+ const v4f_T Rg = v4f_mul(v4f_mul(v4f_set1(0.5f), Dp), tmp); /* [nm] */
+ return Rg;
+}
+
+#endif /* ATRSTM_RDGFA_SIMD4_H */
+
diff --git a/src/test_atrstm.c b/src/test_atrstm.c
@@ -50,8 +50,8 @@ print_help(const char* cmd)
" -f FRACTAL_DIM fractal dimension. Its default value is %g.\n",
ATRSTM_ARGS_DEFAULT.fractal_dimension);
printf(
-" -g PREFACTOR gyration radius prefactor. Its default value is %g.\n",
- ATRSTM_ARGS_DEFAULT.gyration_radius_prefactor);
+" -g PREFACTOR fractal prefactor. Its default value is %g.\n",
+ ATRSTM_ARGS_DEFAULT.fractal_prefactor);
printf(
" -h display this help and exit.\n");
printf(
@@ -163,8 +163,8 @@ args_init(struct args* args, int argc, char** argv)
res = RES_BAD_ARG;
break;
case 'g':
- res = cstr_to_double(optarg, &args->atrstm.gyration_radius_prefactor);
- if(res == RES_OK && args->atrstm.gyration_radius_prefactor <= 0)
+ res = cstr_to_double(optarg, &args->atrstm.fractal_prefactor);
+ if(res == RES_OK && args->atrstm.fractal_prefactor <= 0)
res = RES_BAD_ARG;
break;
case 'h':
diff --git a/src/test_atrstm_radcoefs.c b/src/test_atrstm_radcoefs.c
@@ -28,7 +28,7 @@ main(int argc, char** argv)
args.lambda = 633;
args.n = 1.90;
args.kappa = 0.55;
- args.gyration_radius_prefactor = 1.70;
+ args.fractal_prefactor = 1.70;
args.fractal_dimension = 1.75;
args.soot_volumic_fraction = 1.e-7;
args.soot_primary_particles_count = 100;
diff --git a/src/test_atrstm_radcoefs_simd.c b/src/test_atrstm_radcoefs_simd.c
@@ -13,7 +13,7 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>. */
-#include "atrstm_radcoefs_simd4.c"
+#include "atrstm_radcoefs_simd4.h"
int
main(int argc, char** argv)
@@ -41,7 +41,7 @@ main(int argc, char** argv)
args.lambda = v4f_set1(633.f);
args.n = v4f_set1(1.90f);
args.kappa = v4f_set1(0.55f);
- args.gyration_radius_prefactor = v4f_set1(1.70f);
+ args.fractal_prefactor = v4f_set1(1.70f);
args.fractal_dimension = v4f_set1(1.75f);
args.soot_volumic_fraction = v4f_set(1e-7f, 0.f, 1e-7f, 1e-7f);
args.soot_primary_particles_count = v4f_set(100.f, 100.f, 0.f, 100.f);
@@ -109,7 +109,7 @@ main(int argc, char** argv)
args.lambda = v4f_set1(633.f);
args.n = v4f_set1(1.75f);
args.kappa = v4f_set1(0.435f);
- args.gyration_radius_prefactor = v4f_set1(1.70f);
+ args.fractal_prefactor = v4f_set1(1.70f);
args.fractal_dimension = v4f_set1(1.75f);
args.soot_volumic_fraction = v4f_set
(9.9999999999999995e-008f,
