atrstm

Load and structure a combustion gas mixture
git clone git://git.meso-star.fr/atrstm.git
Log | Files | Refs | README | LICENSE

commit f05b326a6f3b6f24d9083740094a7efcc440c8df
parent 16f9712ee6acaa269a691992046fdf3c13a5e41a
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Wed,  3 Feb 2021 13:27:51 +0100

Remove unused files

Diffstat:
Dsrc/atrstm_radprops.c | 271-------------------------------------------------------------------------------
Dsrc/atrstm_radprops.h | 166-------------------------------------------------------------------------------
2 files changed, 0 insertions(+), 437 deletions(-)

diff --git a/src/atrstm_radprops.c b/src/atrstm_radprops.c @@ -1,271 +0,0 @@ -/* 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_log.h" -#include "atrstm_radcoefs.h" - -#include <astoria/atrri.h> -#include <astoria/atrtp.h> - -#include <star/suvm.h> - -#include <rsys/cstr.h> - -#include <string.h> - -/******************************************************************************* - * Helper function - ******************************************************************************/ -static INLINE int -check_fetch_radcoefs_args - (const struct atrstm* atrstm, - const struct atrstm_fetch_radcoefs_args* args) -{ - int i; - ASSERT(atrstm); - - if(!args - || args->wavelength < atrstm->wlen_range[0] - || args->wavelength > atrstm->wlen_range[1] - || !(args->radcoefs_mask & ATRSTM_RADCOEFS_MASK_ALL) - || !(args->components_mask & ATRSTM_CPNTS_MASK_ALL)) - return 0; - - FOR_EACH(i, 0, ATRSTM_RADCOEFS_COUNT__) { - if(args->k_min[i] > args->k_max[i]) return 0; - } - - return 1; -} - -/******************************************************************************* - * Exported functions - ******************************************************************************/ -res_T -atrstm_fetch_radcoefs - (const struct atrstm* atrstm, - const struct atrstm_fetch_radcoefs_args* args, - double radcoefs[ATRSTM_RADCOEFS_COUNT__]) -{ - struct radcoefs prim_radcoefs[4]; - double bcoords[4]; - struct suvm_primitive prim = SUVM_PRIMITIVE_NULL; - res_T res = RES_OK; - - if(!atrstm || !check_fetch_radcoefs_args(atrstm, args) || !radcoefs) { - res = RES_BAD_ARG; - goto error; - } - 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 belongs 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, prim_radcoefs); - if(res != RES_OK) { - log_err(atrstm, - "Error computing the radiative properties for the primitive %lu -- %s.\n", - (unsigned long)prim.iprim, - res_to_cstr(res)); - goto error; - } - - /* 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]; - 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]; - 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]; - 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; -} - -/******************************************************************************* - * Local functions - ******************************************************************************/ -res_T -primitive_compute_radcoefs - (const struct atrstm* atrstm, - const struct atrri_refractive_index* refract_id, - const struct suvm_primitive* prim, - struct radcoefs radcoefs[]) -{ - struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL; - struct atrtp_desc desc = ATRTP_DESC_NULL; - size_t inode; - res_T res = RES_OK; - ASSERT(atrstm && prim && prim->nvertices == 4 && refract_id && radcoefs); - - res = atrtp_get_desc(atrstm->atrtp, &desc); - if(res != RES_OK) goto error; - - /* Setup the constant input arguments of the optical properties computation - * routine */ - 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_dimension = atrstm->fractal_dimension; - - FOR_EACH(inode, 0, 4) { - const double* node; - - /* Fetch the thermodynamic properties of the node */ - node = atrtp_desc_get_node_properties(&desc, prim->indices[inode]); - - /* Setup the per node input arguments of the optical properties computation - * routine */ - args.soot_volumic_fraction = - node[ATRTP_SOOT_VOLFRAC]; - args.soot_primary_particles_count = - node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT]; - args.soot_primary_particles_diameter = - node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER]; - - /* Compute the node optical properties */ - radcoefs_compute(&radcoefs[inode], &args); - } - -exit: - return res; -error: - goto exit; -} - -res_T -dump_radcoefs - (const struct atrstm* atrstm, - const struct atrri_refractive_index* refract_id, - FILE* stream) -{ - struct radcoefs_compute_args args = RADCOEFS_COMPUTE_ARGS_NULL; - struct atrtp_desc desc = ATRTP_DESC_NULL; - struct radcoefs props_min; - struct radcoefs props_max; - size_t inode; - - res_T res = RES_OK; - ASSERT(atrstm && stream); - - res = atrtp_get_desc(atrstm->atrtp, &desc); - if(res != RES_OK) goto error; - - /* Setup the constant input params of the optical properties computation */ - 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_dimension = atrstm->fractal_dimension; - - props_min.ka = DBL_MAX; - props_max.ka =-DBL_MAX; - props_min.ks = DBL_MAX; - props_max.ks =-DBL_MAX; - props_min.kext = DBL_MAX; - props_max.kext =-DBL_MAX; - - #define FPRINTF(Text, Args) { \ - const int n = fprintf(stream, Text COMMA_##Args LIST_##Args); \ - if(n < 0) { \ - log_err(atrstm, "Error writing optical properties.\n"); \ - res = RES_IO_ERR; \ - goto error; \ - } \ - } (void)0 - - FPRINTF("# Ka Ks Kext\n", ARG0()); - - FOR_EACH(inode, 0, desc.nnodes) { - struct radcoefs props; - const double* node; - - /* Fetch the thermodynamic properties of the node */ - node = atrtp_desc_get_node_properties(&desc, inode); - - /* Setup the per node input args of the optical properties computation */ - args.soot_volumic_fraction = - node[ATRTP_SOOT_VOLFRAC]; - args.soot_primary_particles_count = - node[ATRTP_SOOT_PRIMARY_PARTICLES_COUNT]; - args.soot_primary_particles_diameter = - node[ATRTP_SOOT_PRIMARY_PARTICLES_DIAMETER]; - - /* Compute the node optical properties */ - radcoefs_compute(&props, &args); - - FPRINTF("%g %g %g\n", ARG3(props.ka, props.ks, props.kext)); - - props_min.ka = MMIN(props_min.ka, props.ka); - props_max.ka = MMAX(props_max.ka, props.ka); - props_min.ks = MMIN(props_min.ks, props.ks); - props_max.ks = MMAX(props_max.ks, props.ks); - props_min.kext = MMIN(props_min.kext, props.kext); - props_max.kext = MMAX(props_max.kext, props.kext); - - } - - FPRINTF("# Bounds = [%g %g %g], [%g, %g, %g]\n", ARG6 - (props_min.ka, props_min.ks, props_min.kext, - props_max.ka, props_max.ks, props_max.kext)); - - #undef FPRINTF - -exit: - return res; -error: - goto exit; -} diff --git a/src/atrstm_radprops.h b/src/atrstm_radprops.h @@ -1,166 +0,0 @@ -/* 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_RADCOEFS_H -#define ATRSTM_RADCOEFS_H - -#include <rsys/math.h> -#include <rsys/rsys.h> - -/* Radiative coefficients */ -struct radcoefs { - double ka; /* Absorption coefficient */ - double ks; /* Scattering coefficient */ - double kext; /* Extinction coefficient */ -}; -#define RADCOEFS_NULL__ {0} -static const struct radcoefs RADCOEFS_NULL = RADCOEFS_NULL__; - -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_dimension; - double soot_volumic_fraction; /* In m^3(soot) / m^3 */ - double soot_primary_particles_count; - double soot_primary_particles_diameter; /* In nm */ -}; -#define RADCOEFS_COMPUTE_ARGS_NULL__ {0} -static const struct radcoefs_compute_args RADCOEFS_COMPUTE_ARGS_NULL = - RADCOEFS_COMPUTE_ARGS_NULL__; - -/* Forward declarations */ -struct atrri_refractive_index; -struct atrstm; -struct suvm_primitive; - -extern LOCAL_SYM res_T -primitive_compute_radcoefs - (const struct atrstm* atrstm, - const struct atrri_refractive_index* refract_id, - const struct suvm_primitive* prim, - struct radcoefs radcoefs[]); /* Per primitive node radiative coefficients */ - -/* Write per node radiative coefficients regarding the submitted refracted - * index */ -extern LOCAL_SYM res_T -dump_radcoefs - (const struct atrstm* atrstm, - const struct atrri_refractive_index* refract_id, - FILE* stream); - -static INLINE res_T -primitive_compute_radcoefs_range - (const struct atrstm* atrstm, - const struct atrri_refractive_index* refract_id, - const struct suvm_primitive* prim, - struct radcoefs* radcoefs_min, - struct radcoefs* radcoefs_max) -{ - struct radcoefs props[4]; - res_T res = RES_OK; - ASSERT(radcoefs_min && radcoefs_max); - - res = primitive_compute_radcoefs(atrstm, refract_id, prim, props); - if(res != RES_OK) return res; - - radcoefs_min->ka = MMIN - (MMIN(props[0].ka, props[1].ka), - MMIN(props[2].ka, props[3].ka)); - radcoefs_min->ks = MMIN - (MMIN(props[0].ks, props[1].ks), - MMIN(props[2].ks, props[3].ks)); - radcoefs_min->kext = MMIN - (MMIN(props[0].kext, props[1].kext), - MMIN(props[2].kext, props[3].kext)); - - radcoefs_max->ka = MMAX - (MMAX(props[0].ka, props[1].ka), - MMAX(props[2].ka, props[3].ka)); - radcoefs_max->ks = MMAX - (MMAX(props[0].ks, props[1].ks), - MMAX(props[2].ks, props[3].ks)); - radcoefs_max->kext = MMAX - (MMAX(props[0].kext, props[1].kext), - MMAX(props[2].kext, props[3].kext)); - - return RES_OK; -} - -/* Compute the radiative coefficients of the semi transparent medium as described - * in "Effects of multiple scattering on radiative properties of soot fractal - * aggregates" J. Yon et al., Journal of Quantitative Spectroscopy and - * Radiative Transfer Vol 133, pp. 374-381, 2014 */ -static INLINE void -radcoefs_compute - (struct radcoefs* radcoefs, - const struct radcoefs_compute_args* args) -{ - /* Input variables */ - const double lambda = args->lambda; /* [nm] */ - const double n = args->n; - const double kappa = args->kappa; - 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 Df = args->fractal_dimension; - - if(Np == 0 || Fv == 0) { - radcoefs->ka = 0; - radcoefs->ks = 0; - radcoefs->kext = 0; - } else { - /* Precomputed values */ - const double n2 = n*n; - const double kappa2 = kappa*kappa; - const double four_n2_kappa2 = 4*n2*kappa2; - const double xp = (PI*Dp) / lambda; - const double xp3 = xp*xp*xp; - const double k = (2*PI) / lambda; - const double k2 = k*k; - - /* E(m), m = n + i*kappa */ - const double tmp0 = (n2 - kappa2 + 2); - const double denom = tmp0*tmp0 + four_n2_kappa2; - const double Em = (6*n*kappa) / denom; - - /* F(m), m = n + i*kappa */ - const double tmp1 = ((n2-kappa2 - 1)*(n2-kappa2+2) + four_n2_kappa2)/denom; - const double Fm = tmp1*tmp1 + Em*Em; - - /* Gyration radius */ - const double Rg = 0.5 * Dp * pow(Np/kf, 1.0/Df); /* [nm] */ - const double Rg2 = Rg*Rg; - const double g = pow(1 + 4*k2*Rg2/(3*Df), -Df*0.5); - - /* Cross sections, [nm^3/aggrefate] */ - const double Cabs = Np * (4*PI*xp3)/(k*k) * Em; - const double Csca = Np*Np * (8*PI*xp3*xp3) / (3*k2) * Fm * g; - - /* Radiative coefficients */ - const double Va = (Np*PI*Dp*Dp*Dp) / 6; /* [nm^3] */ - const double rho = Fv / Va; /* [#aggregates / nm^3] */ - const double tmp2 = 1.e-9 * rho; - - ASSERT(denom != 0 && Df != 0 && Dp != 0); - radcoefs->ka = tmp2 * Cabs; - radcoefs->ks = tmp2 * Csca; - radcoefs->kext = radcoefs->ka + radcoefs->ks; - } -} - -#endif /* ATRSTM_RADCOEFS_H */