commit 126da4ebe36becfe0669bcf7b03393d591c465f8 parent 22f2c11ea95188399a306999b699fc448577115f Author: Vincent Forest <vincent.forest@meso-star.com> Date: Thu, 13 Nov 2025 16:53:25 +0100 Merge branch 'release_0.16.2' Diffstat:
127 files changed, 1285 insertions(+), 860 deletions(-)
diff --git a/Makefile b/Makefile @@ -30,6 +30,7 @@ MPI_DEF = -DSDIS_ENABLE_MPI MPI_SRC = src/sdis_mpi.c SRC =\ src/sdis.c\ + src/sdis_brdf.c\ src/sdis_camera.c\ src/sdis_data.c\ src/sdis_device.c\ diff --git a/README.md b/README.md @@ -158,6 +158,13 @@ Edit config.mk as needed, then run: ## Release notes +### Version 0.16.2 + +- Continue to improve numerical robustness when sampling the radiative + path, both in long waves and when processing an external source. +- Update the attenuation of numerical inaccuracies when the candidate + path for reinjection is located approximately on one of the vertices + of a triangle. ### Version 0.16.1 @@ -697,9 +704,8 @@ First version and implementation of the Stardis-Solver API. ## License -Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) Stardis-Solver is free software released under the GPLv3+ license: GNU GPL version 3 or later. You are welcome to redistribute it under certain conditions; refer to the COPYING files for details. - diff --git a/config.mk b/config.mk @@ -1,6 +1,6 @@ VERSION_MAJOR = 0 VERSION_MINOR = 16 -VERSION_PATCH = 1 +VERSION_PATCH = 2 VERSION = $(VERSION_MAJOR).$(VERSION_MINOR).$(VERSION_PATCH) PREFIX = /usr/local diff --git a/src/sdis.c b/src/sdis.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis.h b/src/sdis.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -377,7 +377,7 @@ struct sdis_interface_side_shader { /* Fixed temperature/flux. May be NULL if the temperature/flux is unknown * onto the whole interface */ sdis_interface_getter_T temperature; /* [K]. SDIS_TEMPERATURE_NONE = Unknown */ - sdis_interface_getter_T flux; /* [W.m^-2]. SDIS_FLUX_NONE = no flux */ + sdis_interface_getter_T flux; /* Toward solid. [W.m^-2]. SDIS_FLUX_NONE = no flux */ /* Control the emissivity of the interface. May be NULL for solid/solid * interface or if the emissivity is 0 onto the whole interface. */ @@ -862,7 +862,7 @@ struct sdis_solve_boundary_flux_args { 10000, /* #realisations */ \ NULL, /* List or primitive ids */ \ 0, /* #primitives */ \ - {DBL_MAX,DBL_MAX}, /* Time range */ \ + {DBL_MAX, DBL_MAX}, /* Time range */ \ 1, /* Picard order */ \ NULL, /* RNG state */ \ SSP_RNG_THREEFRY, /* RNG type */ \ diff --git a/src/sdis_Xd_begin.h b/src/sdis_Xd_begin.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_Xd_end.h b/src/sdis_Xd_end.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_brdf.c b/src/sdis_brdf.c @@ -0,0 +1,127 @@ +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) + * + * 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 "sdis.h" +#include "sdis_brdf.h" +#include "sdis_interface_c.h" + +#include <star/ssp.h> +#include <rsys/double3.h> + +/******************************************************************************* + * Helper functions + ******************************************************************************/ +static INLINE res_T +check_brdf_setup_args(const struct brdf_setup_args* args) +{ + const struct sdis_medium* mdm = NULL; + + if(!args) return RES_BAD_ARG; + + if(args->interf == NULL || args->frag == NULL) + return RES_BAD_ARG; + + switch(args->frag->side) { + case SDIS_FRONT: mdm = args->interf->medium_front; break; + case SDIS_BACK: mdm = args->interf->medium_back; break; + default: FATAL("Unreachable code\n"); break; + } + + if(sdis_medium_get_type(mdm) != SDIS_FLUID) + return RES_BAD_ARG; + + return RES_OK; +} + +/* Reflect the V wrt the normal N. By convention V points outward the surface. + * In fact, this function is a double-precision version of the reflect_3d + * function. TODO Clean this "repeat" */ +static FINLINE double* +reflect(double res[3], const double V[3], const double N[3]) +{ + double tmp[3]; + double cos_V_N; + ASSERT(res && V && N); + ASSERT(d3_is_normalized(V) && d3_is_normalized(N)); + cos_V_N = d3_dot(V, N); + d3_muld(tmp, N, 2*cos_V_N); + d3_sub(res, tmp, V); + return res; +} + +/******************************************************************************* + * Local functions + ******************************************************************************/ +res_T +brdf_setup + (struct sdis_device* dev, + const struct brdf_setup_args* args, + struct brdf* brdf) +{ + res_T res = RES_OK; + + ASSERT(check_brdf_setup_args(args) == RES_OK); + + #define GET(Attr) { \ + brdf->Attr = interface_side_get_##Attr \ + (args->interf, args->source_id, args->frag); \ + \ + res = interface_side_check_##Attr \ + (dev, brdf->Attr, args->frag->P, args->frag->time); \ + if(res != RES_OK) goto error; \ + } (void)0 + + GET(emissivity); + GET(specular_fraction); + + #undef GET + +exit: + return res; +error: + *brdf = BRDF_NULL; + goto exit; +} + +void +brdf_sample + (const struct brdf* brdf, + struct ssp_rng* rng, + const double wi[3], /* Incident direction. Point away from the surface */ + const double N[3], /* Surface normal */ + struct brdf_sample* sample) +{ + double r = 0; /* Random number */ + + /* Preconditions */ + ASSERT(brdf && rng && wi && N && sample); + ASSERT(d3_is_normalized(wi) && d3_is_normalized(N)); + ASSERT(d3_dot(wi, N) > 0); + + r = ssp_rng_canonical(rng); + + /* Sample the specular part */ + if(r < brdf->specular_fraction) { + reflect(sample->dir, wi, N); + sample->pdf = 1; + sample->cpnt = BRDF_SPECULAR; + + /* Sample the diffuse part */ + } else { + ssp_ran_hemisphere_cos(rng, N, sample->dir, NULL); + sample->pdf = 1.0/PI; + sample->cpnt = BRDF_DIFFUSE; + } +} diff --git a/src/sdis_brdf.h b/src/sdis_brdf.h @@ -0,0 +1,71 @@ +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) + * + * 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 SDIS_BRDF_H +#define SDIS_BRDF_H + +#include <rsys/rsys.h> + +enum brdf_component { + BRDF_SPECULAR, + BRDF_DIFFUSE, + BRDF_NONE +}; + +struct brdf_sample { + double dir[3]; + double pdf; + enum brdf_component cpnt; +}; +#define BRDF_SAMPLE_NULL__ {{0}, 0, BRDF_NONE} +static const struct brdf_sample BRDF_SAMPLE_NULL = BRDF_SAMPLE_NULL__; + +struct brdf { + double emissivity; + double specular_fraction; +}; +#define BRDF_NULL__ {0, 0} +static const struct brdf BRDF_NULL = BRDF_NULL__; + +/* forward declarations */ +struct sdis_device; +struct sdis_interface; +struct sdis_interface_fragment; +struct ssp_rng; + +struct brdf_setup_args { + struct sdis_interface* interf; + struct sdis_interface_fragment* frag; + unsigned source_id; +}; +#define BRDF_SETUP_ARGS_NULL__ {NULL, NULL, SDIS_INTERN_SOURCE_ID} +static const struct brdf_setup_args BRDF_SETUP_ARGS_NULL = + BRDF_SETUP_ARGS_NULL__; + +extern LOCAL_SYM res_T +brdf_setup + (struct sdis_device* dev, + const struct brdf_setup_args* args, + struct brdf* brdf); + +extern LOCAL_SYM void +brdf_sample + (const struct brdf* brdf, + struct ssp_rng* rng, + const double wi[3], /* Incident direction. Point away from the surface */ + const double N[3], /* Surface normal */ + struct brdf_sample* sample); + +#endif /* SDIS_BRDF_H */ diff --git a/src/sdis_c.h b/src/sdis_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_camera.c b/src/sdis_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_camera.h b/src/sdis_camera.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_data.c b/src/sdis_data.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_device.c b/src/sdis_device.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_estimator.c b/src/sdis_estimator.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_estimator_buffer.c b/src/sdis_estimator_buffer.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_estimator_buffer_X_obs.h b/src/sdis_estimator_buffer_X_obs.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_estimator_buffer_c.h b/src/sdis_estimator_buffer_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_estimator_c.h b/src/sdis_estimator_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_green.c b/src/sdis_green.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_green.h b/src/sdis_green.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path.c b/src/sdis_heat_path.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path.h b/src/sdis_heat_path.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -327,6 +327,51 @@ radiative_path_3d struct ssp_rng* rng, struct temperature* temperature); +extern LOCAL_SYM void +trace_ray_2d + (struct sdis_scene* scn, + const double pos[2], + const double dir[3], /* Always in 3D */ + const double distance, + const unsigned enc_id, + const struct s2d_hit* hit_from, + struct s2d_hit* hit); + +extern LOCAL_SYM void +trace_ray_3d + (struct sdis_scene* scn, + const double pos[3], + const double dir[3], /* Always in 3D */ + const double distance, + const unsigned enc_id, + const struct s3d_hit* hit_from, + struct s3d_hit* hit); + +/* Trace a ray and setup the fragment at the intersection found, if any. */ +extern LOCAL_SYM res_T +find_next_fragment_2d + (struct sdis_scene* scn, + const double in_pos[2], + const double in_dir[3], /* Always in 3D */ + const struct s2d_hit* in_hit, + const double time, + const unsigned enc_id, + struct s2d_hit* out_hit, + struct sdis_interface** out_interf, + struct sdis_interface_fragment* out_frag); + +extern LOCAL_SYM res_T +find_next_fragment_3d + (struct sdis_scene* scn, + const double in_pos[3], + const double in_dir[3], /* Always in 3D */ + const struct s3d_hit* in_hit, + const double time, + const unsigned enc_id, + struct s3d_hit* out_hit, + struct sdis_interface** out_interf, + struct sdis_interface_fragment* out_frag); + /******************************************************************************* * Convective path ******************************************************************************/ diff --git a/src/sdis_heat_path_boundary.c b/src/sdis_heat_path_boundary.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_boundary_Xd.h b/src/sdis_heat_path_boundary_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_boundary_Xd_c.h b/src/sdis_heat_path_boundary_Xd_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -246,130 +246,6 @@ XD(sample_reinjection_dir) #endif } - -#if DIM == 2 -static void -XD(move_away_primitive_boundaries) - (const struct rwalk* rwalk, - const double delta, - double position[DIM]) /* Position to move */ -{ - struct sXd(attrib) attr; - float pos[DIM]; - float dir[DIM]; - float len; - const float st = 0.5f; - ASSERT(rwalk && !SXD_HIT_NONE(&rwalk->XD(hit)) && delta > 0); - - SXD(primitive_get_attrib(&rwalk->XD(hit).prim, SXD_POSITION, st, &attr)); - - fX_set_dX(pos, position); - fX(sub)(dir, attr.value, pos); - len = fX(normalize)(dir, dir); - len = MMIN(len, (float)(delta*0.1)); - - XD(move_pos)(position, dir, len); -} -#else -/* Move the submitted position away from the primitive boundaries to avoid - * numerical issues leading to inconsistent random walks. */ -static void -XD(move_away_primitive_boundaries) - (const struct rwalk* rwalk, - const double delta, - double position[DIM]) -{ - struct s3d_attrib v0, v1, v2; /* Triangle vertices */ - float E[3][4]; /* 3D edge equations */ - float dst[3]; /* Distance from current position to edge equation */ - float N[3]; /* Triangle normal */ - float P[3]; /* Random walk position */ - float tmp[3]; - float min_dst, max_dst; - float cos_a1, cos_a2; - float len; - int imax = 0; - int imin = 0; - int imid = 0; - int i; - ASSERT(rwalk && delta > 0 && !S3D_HIT_NONE(&rwalk->XD(hit))); - - fX_set_dX(P, position); - - /* Fetch triangle vertices */ - S3D(triangle_get_vertex_attrib(&rwalk->XD(hit).prim, 0, S3D_POSITION, &v0)); - S3D(triangle_get_vertex_attrib(&rwalk->XD(hit).prim, 1, S3D_POSITION, &v1)); - S3D(triangle_get_vertex_attrib(&rwalk->XD(hit).prim, 2, S3D_POSITION, &v2)); - - /* Compute the edge vector */ - f3_sub(E[0], v1.value, v0.value); - f3_sub(E[1], v2.value, v1.value); - f3_sub(E[2], v0.value, v2.value); - - /* Compute the triangle normal */ - f3_cross(N, E[1], E[0]); - - /* Compute the 3D edge equation */ - f3_normalize(E[0], f3_cross(E[0], E[0], N)); - f3_normalize(E[1], f3_cross(E[1], E[1], N)); - f3_normalize(E[2], f3_cross(E[2], E[2], N)); - E[0][3] = -f3_dot(E[0], v0.value); - E[1][3] = -f3_dot(E[1], v1.value); - E[2][3] = -f3_dot(E[2], v2.value); - - /* Compute the distance from current position to the edges */ - dst[0] = f3_dot(E[0], P) + E[0][3]; - dst[1] = f3_dot(E[1], P) + E[1][3]; - dst[2] = f3_dot(E[2], P) + E[2][3]; - - /* Retrieve the min and max distance from random walk position to triangle - * edges */ - min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]); - max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]); - - /* Sort the edges with respect to their distance to the random walk position */ - FOR_EACH(i, 0, 3) { - if(dst[i] == min_dst) { - imin = i; - } else if(dst[i] == max_dst) { - imax = i; - } else { - imid = i; - } - } - (void)imax; - - /* TODO if the current position is near a vertex, one should move toward the - * farthest edge along its normal to avoid too small displacement */ - - /* Compute the distance `dst' from the current position to the edges to move - * to, along the normal of the edge from which the random walk is the nearest - * - * +. cos(a) = d / dst => dst = d / cos_a - * / `*. - * / | `*. - * / dst| a /`*. - * / | / `*. - * / | / d `*. - * / |/ `*. - * +---------o----------------+ */ - cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid])); - cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax])); - dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX; - dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX; - - /* Compute the maximum displacement distance into the triangle along the - * normal of the edge from which the random walk is the nearest */ - len = MMIN(dst[imid], dst[imax]); - ASSERT(len != FLT_MAX); - - /* Define the displacement distance as the minimum between 10 percent of - * delta and len / 2. */ - len = MMIN(len*0.5f, (float)(delta*0.1)); - XD(move_pos)(position, E[imin], len); -} -#endif - static res_T XD(find_reinjection_ray) (struct sdis_scene* scn, @@ -411,7 +287,7 @@ XD(find_reinjection_ray) ASSERT(XD(check_find_reinjection_ray_args)(scn, args) == RES_OK); *ray = XD(REINJECTION_RAY_NULL); - MAX_ATTEMPTS = args->can_move ? 2 : 1; + MAX_ATTEMPTS = args->can_move ? 20 : 1; dst_adjusted = args->distance * RAY_RANGE_MAX_SCALE; reinject_threshold = (float)args->distance * REINJECT_DST_MIN_SCALE; @@ -421,7 +297,13 @@ XD(find_reinjection_ray) do { fX_set_dX(org, ray->org); filter_data.XD(hit) = args->rwalk->XD(hit); + + /* Limit the epsilon to 1.e-6, as Star-3D's single-precision floating-point + * representation will inevitably present numerical accuracy problems below + * this threshold. There's no point in going any lower */ + /*filter_data.epsilon = MMAX(args->distance * 0.01, 1e-6);*/ filter_data.epsilon = args->distance * 0.01; + SXD(scene_view_trace_ray (scn->sXd(view), org, args->dir0, range, &filter_data, &hit0)); SXD(scene_view_trace_ray @@ -479,19 +361,15 @@ XD(find_reinjection_ray) * and retry to find a valid reinjection. */ if(dst0 == -1 && dst1 == -1 && iattempt < MAX_ATTEMPTS - 1) { /* Is there still a trial to be done? */ - XD(move_away_primitive_boundaries)(args->rwalk, args->distance, ray->org); + XD(move_away_primitive_boundaries) + (&args->rwalk->XD(hit), args->distance, ray->org); ray->position_was_moved = 1; } } while(dst0 == -1 && dst1 == -1 && ++iattempt < MAX_ATTEMPTS); if(dst0 == -1 && dst1 == -1) { /* No valid reinjection */ -#if DIM == 2 - log_err(scn->dev, "%s: no valid reinjection direction at {%g, %g}.\n", - FUNC_NAME, SPLIT2(ray->org)); -#else - log_err(scn->dev, "%s: no valid reinjection direction at {%g, %g, %g}.\n", - FUNC_NAME, SPLIT3(ray->org)); -#endif + log_err(scn->dev, "%s: no valid reinjection direction at {"FORMAT_VECX"}.\n", + FUNC_NAME, SPLITX(ray->org)); res = RES_BAD_OP_IRRECOVERABLE; goto error; } @@ -1187,4 +1065,130 @@ error: goto exit; } +#if DIM == 2 +void +XD(move_away_primitive_boundaries) + (const struct sXd(hit)* hit, + const double delta, + double position[DIM]) /* Position to move */ +{ + struct sXd(attrib) attr; + float pos[DIM]; + float dir[DIM]; + float len; + const float st = 0.5f; + ASSERT(!SXD_HIT_NONE(hit) && delta > 0); + + SXD(primitive_get_attrib(&hit->prim, SXD_POSITION, st, &attr)); + + fX_set_dX(pos, position); + fX(sub)(dir, attr.value, pos); + len = fX(normalize)(dir, dir); + len = MMIN(len, (float)(delta*0.1)); + + XD(move_pos)(position, dir, len); +} +#else +/* Move the submitted position away from the primitive boundaries to avoid + * numerical issues leading to inconsistent random walks. */ +void +XD(move_away_primitive_boundaries) + (const struct sXd(hit)* hit, + const double delta, + double position[DIM]) +{ + struct s3d_attrib v0, v1, v2; /* Triangle vertices */ + float E[3][4]; /* 3D edge equations */ + float dst[3]; /* Distance from current position to edge equation */ + float N[3]; /* Triangle normal */ + float P[3]; /* Random walk position */ + float tmp[3]; + float min_dst, max_dst; + float len; + int imax = 0; + int imin = 0; + int imid = 0; + int i; + ASSERT(delta > 0 && !S3D_HIT_NONE(hit)); + + fX_set_dX(P, position); + + /* Fetch triangle vertices */ + S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0)); + S3D(triangle_get_vertex_attrib(&hit->prim, 1, S3D_POSITION, &v1)); + S3D(triangle_get_vertex_attrib(&hit->prim, 2, S3D_POSITION, &v2)); + + /* Compute the edge vector */ + f3_sub(E[0], v1.value, v0.value); + f3_sub(E[1], v2.value, v1.value); + f3_sub(E[2], v0.value, v2.value); + + /* Compute the triangle normal */ + f3_cross(N, E[1], E[0]); + + /* Compute the 3D edge equation */ + f3_normalize(E[0], f3_cross(E[0], E[0], N)); + f3_normalize(E[1], f3_cross(E[1], E[1], N)); + f3_normalize(E[2], f3_cross(E[2], E[2], N)); + E[0][3] = -f3_dot(E[0], v0.value); + E[1][3] = -f3_dot(E[1], v1.value); + E[2][3] = -f3_dot(E[2], v2.value); + + /* Compute the distance from current position to the edges */ + dst[0] = f3_dot(E[0], P) + E[0][3]; + dst[1] = f3_dot(E[1], P) + E[1][3]; + dst[2] = f3_dot(E[2], P) + E[2][3]; + + /* Retrieve the min and max distance from random walk position to triangle + * edges */ + min_dst = MMIN(MMIN(dst[0], dst[1]), dst[2]); + max_dst = MMAX(MMAX(dst[0], dst[1]), dst[2]); + + /* Sort the edges with respect to their distance to the random walk position */ + FOR_EACH(i, 0, 3) { + if(dst[i] == min_dst) { + imin = i; + } else if(dst[i] == max_dst) { + imax = i; + } else { + imid = i; + } + } + (void)imax; + + if(eq_eps(dst[imin], 0, delta*1e-3) && eq_eps(dst[imid], 0, delta*1e-3)) { + /* The random position is in a corner, meaning that its distance to the two + * nearest edges is approximately equal to 0. Move it toward the farthest + * edge along its normal to avoid moving too little. */ + len = MMIN(dst[imax]*0.5f, (float)delta*0.1f); + XD(move_pos)(position, f3_minus(tmp, E[imax]), len); + + } else { + /* Compute the distance `dst' from the current position to the edges to move + * to, along the normal of the edge from which the random walk is the nearest + * + * +. cos(a) = d / dst => dst = d / cos_a + * / `*. + * / | `*. + * / dst| a /`*. + * / | / `*. + * / | / d `*. + * / |/ `*. + * +---------o----------------+ */ + const float cos_a1 = f3_dot(E[imin], f3_minus(tmp, E[imid])); + const float cos_a2 = f3_dot(E[imin], f3_minus(tmp, E[imax])); + dst[imid] = cos_a1 > 0 ? dst[imid] / cos_a1 : FLT_MAX; + dst[imax] = cos_a2 > 0 ? dst[imax] / cos_a2 : FLT_MAX; + len = MMIN(dst[imid], dst[imax]); + ASSERT(len != FLT_MAX); + + /* Define the displacement distance as the minimum between 10 percent of + * delta and len / 2. */ + len = MMIN(len*0.5f, (float)(delta*0.1)); + XD(move_pos)(position, E[imin], len); + } + +} +#endif + #include "sdis_Xd_end.h" diff --git a/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h b/src/sdis_heat_path_boundary_Xd_handle_external_net_flux.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -13,6 +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 "sdis_brdf.h" #include "sdis_heat_path_boundary_c.h" #include "sdis_interface_c.h" #include "sdis_log.h" @@ -29,27 +30,6 @@ #ifndef SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H #define SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H -enum brdf_component { - BRDF_SPECULAR, - BRDF_DIFFUSE, - BRDF_NONE -}; - -struct brdf_sample { - double dir[3]; - double pdf; - enum brdf_component cpnt; -}; -#define BRDF_SAMPLE_NULL__ {{0}, 0, BRDF_NONE} -static const struct brdf_sample BRDF_SAMPLE_NULL = BRDF_SAMPLE_NULL__; - -struct brdf { - double emissivity; - double specular_fraction; -}; -#define BRDF_NULL__ {0, 0} -static const struct brdf BRDF_NULL = BRDF_NULL__; - /* Incident diffuse flux is made up of two components. One corresponds to the * diffuse flux due to the reflection of the source on surfaces. The other is * the diffuse flux due to the source's radiation scattering at least once in @@ -63,107 +43,6 @@ struct incident_diffuse_flux { static const struct incident_diffuse_flux INCIDENT_DIFFUSE_FLUX_NULL = INCIDENT_DIFFUSE_FLUX_NULL__; -/* Reflect the V wrt the normal N. By convention V points outward the surface. - * In fact, this function is a double-precision version of the reflect_3d - * function. TODO Clean this "repeat" */ -static FINLINE double* -reflect(double res[3], const double V[3], const double N[3]) -{ - double tmp[3]; - double cos_V_N; - ASSERT(res && V && N); - ASSERT(d3_is_normalized(V) && d3_is_normalized(N)); - cos_V_N = d3_dot(V, N); - d3_muld(tmp, N, 2*cos_V_N); - d3_sub(res, tmp, V); - return res; -} - -static void -sample_brdf - (const struct brdf* brdf, - struct ssp_rng* rng, - const double wi[3], /* Incident direction. Point away from the surface */ - const double N[3], /* Surface normal */ - struct brdf_sample* sample) -{ - double r = 0; /* Random number */ - - /* Preconditions */ - ASSERT(brdf && rng && wi && N && sample); - ASSERT(d3_is_normalized(wi) && d3_is_normalized(N)); - ASSERT(d3_dot(wi, N) > 0); - - r = ssp_rng_canonical(rng); - - /* Sample the specular part */ - if(r < brdf->specular_fraction) { - reflect(sample->dir, wi, N); - sample->pdf = 1; - sample->cpnt = BRDF_SPECULAR; - - /* Sample the diffuse part */ - } else { - ssp_ran_hemisphere_cos(rng, N, sample->dir, NULL); - sample->pdf = 1.0/PI; - sample->cpnt = BRDF_DIFFUSE; - } -} - -/* Check that the trajectory reaches a valid interface, i.e. that it is on a - * fluid/solid interface and has reached it from the fluid */ -static res_T -check_interface - (const struct sdis_interface* interf, - const struct sdis_interface_fragment* frag) -{ - enum sdis_medium_type mdm_frt_type = SDIS_MEDIUM_TYPES_COUNT__; - enum sdis_medium_type mdm_bck_type = SDIS_MEDIUM_TYPES_COUNT__; - enum sdis_side fluid_side = SDIS_SIDE_NULL__; - res_T res = RES_OK; - - mdm_frt_type = sdis_medium_get_type(interf->medium_front); - mdm_bck_type = sdis_medium_get_type(interf->medium_back); - - /* Semi-transparent materials are not supported. This means that a solid/solid - * interface must not be intersected when tracing radiative paths */ - if(mdm_frt_type == SDIS_SOLID && mdm_bck_type == SDIS_SOLID) { - log_err(interf->dev, - "Error when sampling the trajectory to calculate the incident diffuse " - "flux. The trajectory reaches a solid/solid interface, whereas this is " - "supposed to be impossible (path position: %g, %g, %g).\n", - SPLIT3(frag->P)); - res = RES_BAD_OP; - goto error; - } - - /* Find out which side of the interface the fluid is on */ - if(mdm_frt_type == SDIS_FLUID) { - fluid_side = SDIS_FRONT; - } else if(mdm_bck_type == SDIS_FLUID) { - fluid_side = SDIS_BACK; - } else { - FATAL("Unreachable code\n"); - } - - /* Check that the current position is on the correct side of the interface */ - if(frag->side != fluid_side) { - log_err(interf->dev, - "Inconsistent intersection when sampling the trajectory to calculate the " - "incident diffuse flux. The radiative path reaches an interface on " - "its solid side, whereas this is supposed to be impossible " - "(path position: %g, %g, %g).\n", - SPLIT3(frag->P)); - res = RES_BAD_OP; - goto error; - } - -exit: - return res; -error: - goto exit; -} - #endif /* SDIS_HEAT_PATH_BOUNDARY_XD_HANDLE_EXTERNAL_NET_FLUX_H */ /******************************************************************************* @@ -209,48 +88,19 @@ XD(check_handle_external_net_flux_args) return RES_OK; } -static INLINE void -XD(trace_ray) - (const struct sdis_scene* scn, - const double pos[DIM], - const double dir[3], - const double distance, - const struct sXd(hit)* hit_from, - struct sXd(hit)* hit) -{ - struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL; - float ray_org[DIM] = {0}; - float ray_dir[3] = {0}; - float ray_range[2] = {0}; - ASSERT(scn && pos && dir && distance >= 0 && hit_from && hit); - - fX_set_dX(ray_org, pos); - f3_set_d3(ray_dir, dir); - ray_range[0] = 0; - ray_range[1] = (float)distance; - filter_data.XD(hit) = *hit_from; - filter_data.epsilon = 1.e-4; -#if DIM == 2 - SXD(scene_view_trace_ray_3d - (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit)); -#else - SXD(scene_view_trace_ray - (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit)); -#endif -} - static INLINE double /* [W/m^2/sr] */ XD(direct_contribution) - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct source_sample* sample, const double pos[DIM], + const unsigned enc_id, /* Current enclosure */ const struct sXd(hit)* hit_from) { struct sXd(hit) hit = SXD_HIT_NULL; ASSERT(scn && sample && pos && hit_from); /* Is the source hidden */ - XD(trace_ray)(scn, pos, sample->dir, sample->dst, hit_from, &hit); + XD(trace_ray)(scn, pos, sample->dir, sample->dst, enc_id, hit_from, &hit); if(!SXD_HIT_NONE(&hit)) return 0; /* [W/m^2/sr] */ /* Note that the value returned is not the source's actual radiance, but the @@ -261,71 +111,15 @@ XD(direct_contribution) return sample->radiance_term; /* [W/m^2/sr] */ } -static INLINE void -XD(setup_fragment) - (struct sdis_interface_fragment* frag, - const double pos[DIM], - const double dir[DIM], /* Direction _toward_ the hit position */ - const double time, /* Current time */ - const double N[DIM],/* Surface normal */ - const struct sXd(hit)* hit) -{ - struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; - enum sdis_side side = SDIS_SIDE_NULL__; - ASSERT(frag && pos && dir && N); - ASSERT(dX(is_normalized)(N)); - - /* Setup the interface fragment at the intersection position */ - dX(set)(vtx.P, pos); - vtx.time = time; - side = dX(dot)(dir, N) < 0 ? SDIS_FRONT : SDIS_BACK; - XD(setup_interface_fragment)(frag, &vtx, hit, side); -} - -static INLINE res_T -XD(setup_brdf) - (struct sdis_device* dev, - const struct sdis_source* src, - struct brdf* brdf, - const struct sdis_interface* interf, - const struct sdis_interface_fragment* frag) -{ - double epsilon = 0; - double alpha = 0; - unsigned src_id = 0; - res_T res = RES_OK; - ASSERT(brdf && frag); - ASSERT((frag->side == SDIS_FRONT - && sdis_medium_get_type(interf->medium_front) == SDIS_FLUID) - || sdis_medium_get_type(interf->medium_back) == SDIS_FLUID); - - src_id = sdis_source_get_id(src); - - epsilon = interface_side_get_emissivity(interf, src_id, frag); - res = interface_side_check_emissivity(dev, epsilon, frag->P, frag->time); - if(res != RES_OK) goto error; - - alpha = interface_side_get_specular_fraction(interf, src_id, frag); - res = interface_side_check_specular_fraction(dev, alpha, frag->P, frag->time); - if(res != RES_OK) goto error; - - brdf->emissivity = epsilon; - brdf->specular_fraction = alpha; - -exit: - return res; -error: - *brdf = BRDF_NULL; - goto exit; -} - static res_T XD(compute_incident_diffuse_flux) - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct ssp_rng* rng, + const struct source_props* props, const double in_pos[DIM], /* position */ const double in_N[DIM], /* Surface normal. (Away from the surface) */ const double time, + const unsigned enc_id, /* Current enclosure */ const struct sXd(hit)* in_hit, /* Current intersection */ struct incident_diffuse_flux* diffuse_flux) /* [W/m^2] */ { @@ -333,8 +127,9 @@ XD(compute_incident_diffuse_flux) double pos[3] = {0}; /* In 3D for ray tracing ray to the source */ double dir[3] = {0}; /* Incident direction (toward the surface). Always 3D.*/ double N[3] = {0}; /* Surface normal. Always 3D */ + size_t nbounces = 0; /* For debug */ res_T res = RES_OK; - ASSERT(in_pos && in_N && in_hit && diffuse_flux); + ASSERT(props && in_pos && in_N && in_hit && diffuse_flux); /* Local copy of input argument */ dX(set)(pos, in_pos); @@ -348,7 +143,7 @@ XD(compute_incident_diffuse_flux) for(;;) { /* External sources */ - struct source_sample src_sample = SOURCE_SAMPLE_NULL; + struct source_sample samp = SOURCE_SAMPLE_NULL; /* Interface */ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; @@ -356,17 +151,19 @@ XD(compute_incident_diffuse_flux) /* BRDF */ struct brdf brdf = BRDF_NULL; - struct brdf_sample brdf_sample = BRDF_SAMPLE_NULL; + struct brdf_sample bounce = BRDF_SAMPLE_NULL; + struct brdf_setup_args brdf_setup_args = BRDF_SETUP_ARGS_NULL; /* Miscellaneous */ double L = 0; /* incident flux to bounce position */ double wi[3] = {0}; /* Incident direction (outward the surface). Always 3D */ - double vec[DIM] = {0}; /* Temporary variable */ d3_minus(wi, dir); /* Always in 3D */ - /* Find the following surface along the direction of propagation */ - XD(trace_ray)(scn, pos, dir, INF, &hit, &hit); + res = XD(find_next_fragment) + (scn, pos, dir, &hit, time, enc_id, &hit, &interf, &frag); + if(res != RES_OK) goto error; + if(SXD_HIT_NONE(&hit)) { /* No surface. Handle the radiance emitted by the source and scattered at * least once in the environment. Note that the value returned is not the @@ -382,49 +179,47 @@ XD(compute_incident_diffuse_flux) break; } - /* Retrieve the current position and normal */ - dX(add)(pos, pos, dX(muld)(vec, dir, hit.distance)); - dX_set_fX(N, hit.normal); - dX(normalize(N, N)); - - /* Retrieve the current interface properties */ - interf = scene_get_interface(scn, hit.prim.prim_id); - XD(setup_fragment)(&frag, pos, dir, time, N, &hit); + d3_set(pos, frag.P); - /* Check that the path reaches a valid interface */ - res = check_interface(interf, &frag); + /* Retrieve BRDF at current interface position */ + brdf_setup_args.interf = interf; + brdf_setup_args.frag = &frag; + brdf_setup_args.source_id = sdis_source_get_id(scn->source); + res = brdf_setup(scn->dev, &brdf_setup_args, &brdf); if(res != RES_OK) goto error; - XD(setup_brdf)(scn->dev, scn->source, &brdf, interf, &frag); - /* Check if path is absorbed */ if(ssp_rng_canonical(rng) < brdf.emissivity) break; /* Sample rebound direction */ - if(frag.side == SDIS_BACK) dX(minus)(N, N); /* Revert normal if necessary */ - sample_brdf(&brdf, rng, wi, N, &brdf_sample); - d3_set(dir, brdf_sample.dir); /* Always in 3D */ + switch(frag.side) { + case SDIS_FRONT: dX(set)(N, frag.Ng); break; + case SDIS_BACK: dX(minus)(N, frag.Ng); break; + default: FATAL("Unreachable code\n"); + } + brdf_sample(&brdf, rng, wi, N, &bounce); + d3_set(dir, bounce.dir); /* Always in 3D */ /* Calculate the direct contribution if the rebound is specular */ - if(brdf_sample.cpnt == BRDF_SPECULAR) { - res = source_trace_to(scn->source, pos, brdf_sample.dir, time, &src_sample); + if(bounce.cpnt == BRDF_SPECULAR) { + res = source_trace_to(scn->source, props, pos, bounce.dir, &samp); if(res != RES_OK) goto error; - if(!SOURCE_SAMPLE_NONE(&src_sample)) { - const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit); + if(!SOURCE_SAMPLE_NONE(&samp)) { + double Ld = XD(direct_contribution)(scn, &samp, pos, enc_id, &hit); L = Ld; /* [W/m^2/sr] */ } /* Calculate the direct contribution of the rebound is diffuse */ } else { double cos_theta = 0; - ASSERT(brdf_sample.cpnt == BRDF_DIFFUSE); + ASSERT(bounce.cpnt == BRDF_DIFFUSE); /* Sample an external source to handle its direct contribution at the * bounce position */ - res = source_sample(scn->source, rng, pos, time, &src_sample); + res = source_sample(scn->source, props, rng, pos, &samp); CHK(res == RES_OK); - cos_theta = d3_dot(src_sample.dir, N); + cos_theta = d3_dot(samp.dir, N); /* The source is behind the surface */ if(cos_theta <= 0) { @@ -432,11 +227,12 @@ XD(compute_incident_diffuse_flux) /* The source is above the surface */ } else { - const double Ld = XD(direct_contribution)(scn, &src_sample, pos, &hit); - L = Ld * cos_theta / (PI * src_sample.pdf); /* [W/m^2/sr] */ + double Ld = XD(direct_contribution)(scn, &samp, pos, enc_id, &hit); + L = Ld * cos_theta / (PI * samp.pdf); /* [W/m^2/sr] */ } } diffuse_flux->reflected += L; /* [W/m^2/sr] */ + ++nbounces; } diffuse_flux->reflected *= PI; /* [W/m^2] */ @@ -451,7 +247,7 @@ error: ******************************************************************************/ res_T XD(handle_external_net_flux) - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct ssp_rng* rng, const struct handle_external_net_flux_args* args, struct temperature* T) @@ -460,7 +256,8 @@ XD(handle_external_net_flux) struct sdis_green_external_flux_terms green = SDIS_GREEN_EXTERNAL_FLUX_TERMS_NULL; - /* Sampling external sources */ + /* External source */ + struct source_props src_props = SOURCE_PROPS_NULL; struct source_sample src_sample = SOURCE_SAMPLE_NULL; /* External flux */ @@ -477,6 +274,7 @@ XD(handle_external_net_flux) struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; /* Miscellaneous */ + unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; double sum_h = 0; double emissivity = 0; /* Emissivity */ double Ld = 0; /* Incident radiance [W/m^2/sr] */ @@ -489,7 +287,7 @@ XD(handle_external_net_flux) res = XD(check_handle_external_net_flux_args)(scn, FUNC_NAME, args); if(res != RES_OK) goto error; - /* Setup the interface fragment on flud side */ + /* Setup the interface fragment on fluid side */ frag = *args->frag; if(sdis_medium_get_type(args->interf->medium_front) == SDIS_FLUID) { frag.side = SDIS_FRONT; @@ -498,6 +296,9 @@ XD(handle_external_net_flux) frag.side = SDIS_BACK; } + /* Retrieve the enclosures */ + scene_get_enclosure_ids(scn, args->XD(hit)->prim.prim_id, enc_ids); + /* No external sources <=> no external fluxes. Nothing to do */ handle_flux = interface_side_is_external_flux_handled(args->interf, &frag); handle_flux = handle_flux && (scn->source != NULL); @@ -508,11 +309,14 @@ XD(handle_external_net_flux) emissivity = interface_side_get_emissivity(args->interf, src_id, &frag); res = interface_side_check_emissivity(scn->dev, emissivity, frag.P, frag.time); if(res != RES_OK) goto error; + if(emissivity == 0) goto exit; - /* Sample the external source */ - res = source_sample - (scn->source, rng, frag.P, frag.time, &src_sample); + res = source_get_props(scn->source, frag.time, &src_props); + if(res != RES_OK) goto error; + + /* Sample a direction toward the source to add its direct contribution */ + res = source_sample(scn->source, &src_props, rng, frag.P, &src_sample); if(res != RES_OK) goto error; /* Setup the normal to ensure that it points toward the fluid medium */ @@ -523,13 +327,14 @@ XD(handle_external_net_flux) * interface side */ cos_theta = d3_dot(N, src_sample.dir); if(cos_theta > 0) { - Ld = XD(direct_contribution)(scn, &src_sample, frag.P, args->XD(hit)); + Ld = XD(direct_contribution) + (scn, &src_sample, frag.P, enc_ids[frag.side], args->XD(hit)); incident_flux_direct = cos_theta * Ld / src_sample.pdf; /* [W/m^2] */ } /* Calculate the incident diffuse flux [W/m^2] */ - res = XD(compute_incident_diffuse_flux) - (scn, rng, frag.P, N, frag.time, args->XD(hit), &incident_flux_diffuse); + res = XD(compute_incident_diffuse_flux)(scn, rng, &src_props, frag.P, N, + frag.time, enc_ids[frag.side], args->XD(hit), &incident_flux_diffuse); if(res != RES_OK) goto error; /* Calculate the incident flux without the part scattered by the environment. @@ -562,7 +367,7 @@ XD(handle_external_net_flux) green.dir[1] = incident_flux_diffuse.dir[1]; green.dir[2] = incident_flux_diffuse.dir[2]; - T->value += green.term_wrt_power * source_get_power(scn->source, green.time); + T->value += green.term_wrt_power * src_props.power; if(green.term_wrt_diffuse_radiance) { T->value += green.term_wrt_diffuse_radiance diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picard1.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -48,6 +48,7 @@ XD(rwalk_get_Tref) ray.dir[0] = rwalk->dir[0]; ray.dir[1] = rwalk->dir[1]; ray.dir[2] = rwalk->dir[2]; + ray.time = rwalk->vtx.time; Tref = radiative_env_get_reference_temperature(scn->radenv, &ray); } else { struct sdis_interface_fragment frag; diff --git a/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h b/src/sdis_heat_path_boundary_Xd_solid_fluid_picardN.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_boundary_Xd_solid_solid.h b/src/sdis_heat_path_boundary_Xd_solid_solid.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_boundary_c.h b/src/sdis_heat_path_boundary_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -175,14 +175,14 @@ HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL = HANDLE_EXTERNAL_NET_FLUX_ARGS_NULL__; extern LOCAL_SYM res_T handle_external_net_flux_2d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct ssp_rng* rng, const struct handle_external_net_flux_args* args, struct temperature* T); extern LOCAL_SYM res_T handle_external_net_flux_3d - (const struct sdis_scene* scn, + (struct sdis_scene* scn, struct ssp_rng* rng, const struct handle_external_net_flux_args* args, struct temperature* T); @@ -222,6 +222,20 @@ query_medium_temperature_from_boundary_3d struct rwalk* rwalk, struct temperature* T); +/* Move the submitted position away from the primitive boundaries to avoid + * numerical issues leading to inconsistent random walks. */ +extern LOCAL_SYM void +move_away_primitive_boundaries_2d + (const struct s2d_hit* hit, + const double delta, + double position[2]); /* Position to move */ + +extern LOCAL_SYM void +move_away_primitive_boundaries_3d + (const struct s3d_hit* hit, + const double delta, + double position[3]); /* Position to move */ + /******************************************************************************* * Boundary sub-paths ******************************************************************************/ diff --git a/src/sdis_heat_path_conductive.c b/src/sdis_heat_path_conductive.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_conductive_Xd.h b/src/sdis_heat_path_conductive_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_conductive_c.h b/src/sdis_heat_path_conductive_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_conductive_custom_Xd.h b/src/sdis_heat_path_conductive_custom_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_conductive_delta_sphere_Xd.h b/src/sdis_heat_path_conductive_delta_sphere_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_conductive_wos_Xd.h b/src/sdis_heat_path_conductive_wos_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -22,6 +22,9 @@ #include "sdis_Xd_begin.h" +/* Define epsilon shell from delta */ +#define EPSILON_SHELL(Delta) ((Delta)*1e-2) + /******************************************************************************* * Non generic helper functions ******************************************************************************/ @@ -218,6 +221,7 @@ exit: static INLINE enum sdis_side compute_hit_side_2d (const struct s2d_hit* hit, + const double delta, const double pos[2]) /* Position from which intersection occurs */ { struct s2d_attrib p0, p1; /* Segment positions */ @@ -226,7 +230,14 @@ compute_hit_side_2d double z = 0; /* Check pre-conditions */ - ASSERT(hit && pos && !S2D_HIT_NONE(hit)); + ASSERT(hit && delta > 0 && pos && !S2D_HIT_NONE(hit)); + + /* Delta is not yet used. It could be used to check confidence in the impact + * side calculation. A small value of Z (in relation to epsilon) means that the + * position of the input is close to the boundary and that the calculation of + * the side must be carried out with care. So far, however, no such problems + * have been observed in 2D. */ + (void)delta; /* Retrieve the positions of the intersected segment */ S2D(segment_get_vertex_attrib(&hit->prim, 0, S2D_POSITION, &p0)); @@ -245,9 +256,13 @@ compute_hit_side_2d #endif #if DIM == 3 +/* May return SDIS_SIDE_NULL__ if the side cannot be calculated with confidence, + * i.e. if the input position is too close to the boundary and the calculation + * may therefore present numerical problems */ static INLINE enum sdis_side compute_hit_side_3d (const struct s3d_hit* hit, + const double delta, const double pos[3]) /* Position from which intersection occurs */ { struct s3d_attrib v0; /* Position of the 1st triangle vertex */ @@ -257,7 +272,13 @@ compute_hit_side_3d double dst = 0; /* Distance of pos to the plane */ /* Check pre-conditions */ - ASSERT(hit && pos && !S3D_HIT_NONE(hit)); + ASSERT(hit && delta > 0 && pos && !S3D_HIT_NONE(hit)); + + /* The distance is close to the border and its calculation can suffer from + * numerical problems. No side can therefore be estimated with confidence */ + if(hit->distance < 1e-4 && eq_eps(hit->distance, 0, EPSILON_SHELL(delta))) { + return SDIS_SIDE_NULL__; + } /* Retrieve the positions of the intersected triangle */ S3D(triangle_get_vertex_attrib(&hit->prim, 0, S3D_POSITION, &v0)); @@ -287,6 +308,7 @@ XD(check_diffusion_position) enum sdis_side side = SDIS_SIDE_NULL__; struct sXd(hit) hit = SXD_HIT_NULL; + struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL; float wos_pos[DIM] = {0}; float wos_radius = 0; res_T res = RES_OK; @@ -295,6 +317,11 @@ XD(check_diffusion_position) ASSERT(scn && pos); ASSERT(expected_enc_id != ENCLOSURE_ID_NULL); + /* Filter positions on/near a primitive boundary that don't look towards the + * query position */ + filter_data.scn = scn; + filter_data.enc_id = expected_enc_id; + /* Look for the nearest surface of the position to be checked. By limiting the * search radius to delta we speed up the closest point query. If no surface * is found, we assume that the position is in the intended medium. We rely @@ -309,17 +336,27 @@ XD(check_diffusion_position) * problem of this kind could have arisen. */ wos_radius = (float)delta; fX_set_dX(wos_pos, pos); - SXD(scene_view_closest_point(scn->sXd(view), wos_pos, wos_radius, NULL, &hit)); + SXD(scene_view_closest_point + (scn->sXd(view), wos_pos, wos_radius, &filter_data, &hit)); if(SXD_HIT_NONE(&hit)) goto exit; /* Check path consistency */ scene_get_enclosure_ids(scn, hit.prim.prim_id, enc_ids); - side = XD(compute_hit_side)(&hit, pos); - if(enc_ids[side] != expected_enc_id) { + side = XD(compute_hit_side)(&hit, delta, pos); + if(side != SDIS_SIDE_NULL__ && enc_ids[side] != expected_enc_id) { res = RES_BAD_ARG; goto error; } + /* Position is close of the border: check both sides to handle numerical + * problems in calculating hit side */ + if(side == SDIS_SIDE_NULL__ + && enc_ids[SDIS_FRONT] != expected_enc_id + && enc_ids[SDIS_BACK] != expected_enc_id) { + res = RES_BAD_ARG; + goto error; + } + exit: return res; error: @@ -330,6 +367,7 @@ static res_T XD(setup_hit_wos) (struct sdis_scene* scn, const struct sXd(hit)* hit, + const double delta, struct rwalk* rwalk) { /* Geometry */ @@ -355,18 +393,30 @@ XD(setup_hit_wos) SXD(primitive_get_attrib(&prim, SXD_POSITION, hit->uv, &attr)); #endif + scene_get_enclosure_ids(scn, hit->prim.prim_id, enc_ids); + /* Calculate on which side the intersection occurs */ dX_set_fX(tgt, attr.value); - side = XD(compute_hit_side)(hit, rwalk->vtx.P); + side = XD(compute_hit_side)(hit, delta, rwalk->vtx.P); + + /* Calculating the side of the intersection can suffer from numerical problems + * when the position of the path is close to the intersected surface (i.e. + * side == SDIS_SIDE_NULL). It is therefore reasonable to assume that there is + * no cause for concern as long as the enclosure identifier on the other side + * of the triangle is the expected one. */ + if(side == SDIS_SIDE_NULL__) { + if(rwalk->enc_id == enc_ids[SDIS_FRONT]) side = SDIS_FRONT; + if(rwalk->enc_id == enc_ids[SDIS_BACK]) side = SDIS_BACK; + } /* Check path consistency */ - scene_get_enclosure_ids(scn, hit->prim.prim_id, enc_ids); - if(enc_ids[side] != rwalk->enc_id) { + if(side == SDIS_SIDE_NULL__ || enc_ids[side] != rwalk->enc_id) { + res = RES_BAD_OP_IRRECOVERABLE; log_err(scn->dev, "%s:%s: the conductive path has reached an invalid interface. " "Unexpected enclosure -- pos=("FORMAT_VECX"), side=%s\n", - __FILE__, FUNC_NAME, SPLITX(tgt), side == SDIS_FRONT ? "front" : "back"); - res = RES_BAD_OP_IRRECOVERABLE; + __FILE__, FUNC_NAME, SPLITX(tgt), + side == SDIS_FRONT ? "front" : "back"); goto error; } @@ -415,11 +465,7 @@ XD(setup_hit_rt) /* Fetch interface properties and check path consistency */ scene_get_enclosure_ids(scn, hit->prim.prim_id, enc_ids); if(enc_ids[side] != rwalk->enc_id) { - log_err(scn->dev, - "%s:%s: the conductive path has reached an invalid interface. " - "Unexpected enclosure -- pos=("FORMAT_VECX"), side=%s\n", - __FILE__, FUNC_NAME, SPLITX(tgt), side == SDIS_FRONT ? "front" : "back"); - res = RES_BAD_OP_IRRECOVERABLE; + res = RES_BAD_OP; goto error; } @@ -465,11 +511,11 @@ XD(sample_next_position) CHK(!SXD_HIT_NONE(&hit)); wos_distance = hit.distance; - /* The current position is in the epsilon shell, i.e. 1% of delta: + /* The current position is in the epsilon shell, * move it to the nearest interface position */ - wos_epsilon = delta*1.e-2; + wos_epsilon = EPSILON_SHELL(delta); if(wos_distance <= wos_epsilon) { - res = XD(setup_hit_wos)(scn, &hit, rwalk); + res = XD(setup_hit_wos)(scn, &hit, delta, rwalk); if(res != RES_OK) goto error; /* Uniformly sample a new position on the surrounding sphere */ @@ -508,6 +554,7 @@ XD(sample_next_position) * want to move to is actually inside the epsilon shell. In this case, the * trajectory will be moved to this interface in the next step anyway. */ } else { + struct sXd(hit) hit_rt = SXD_HIT_NULL; float rt_pos[DIM] = {0}; float rt_dir[DIM] = {0}; float rt_range[2] = {0, 0}; @@ -516,27 +563,33 @@ XD(sample_next_position) fX_set_dX(rt_dir, dir); rt_range[0] = 0; rt_range[1] = (float)INF; - SXD(scene_view_trace_ray(scn->sXd(view), rt_pos, rt_dir, rt_range, NULL, &hit)); - - /* An intersection should be found. If not, we can do nothing and simply - * reject the path. - * - * TODO: we could take the treatment of numerical problems a step further - * by sampling other directions and trying to move in them again. But at - * present, and until we have proof to the contrary, we assume that the - * rejection of a path should not occur, or that it will be so rare that - * we don't care to save it. */ - if(SXD_HIT_NONE(&hit)) { - log_err(scn->dev, - "%s:%s: unable to find the next diffusion position -- " - "position=("FORMAT_VECX"), direction=("FORMAT_VECX"), distance=%g\n", - __FILE__, FUNC_NAME, SPLITX(pos), SPLITX(dir), wos_distance); - res = RES_BAD_OP_IRRECOVERABLE; - goto error; + SXD(scene_view_trace_ray + (scn->sXd(view), rt_pos, rt_dir, rt_range, NULL, &hit_rt)); + + if(SXD_HIT_NONE(&hit_rt)) { + /* The lack of intersection is probably due to a current position close + * to the boundary. And although it is detected in the solid by WoS, the + * specific numerical errors of the ray-tracing operator may contradict + * the WoS algorithm, which relies on the closest point operator. But + * since the position is close to the boundary, it can be snaped to it*/ + res = XD(setup_hit_wos)(scn, &hit, delta, rwalk); + if(res != RES_OK) goto error; + + } else { + res = XD(setup_hit_rt)(scn, rwalk->vtx.P, dir, &hit_rt, rwalk); + if(res != RES_OK) { + /* An error occurs while handling ray intersection. This means that + * the Ray-Tracing operator find an invalid intersection regarding the + * current enclosure in which the path should be sampled. As in the + * case of the lack of intersection (see above) this means that the + * position is close to the enclosure boundary and that the ray missed + * it. So, As previously, the position can be simply snaped to it + * since one can assumes that the current position is in the right + * enclosure */ + res = XD(setup_hit_wos)(scn, &hit, delta, rwalk); + if(res != RES_OK) goto error; + } } - - res = XD(setup_hit_rt)(scn, rwalk->vtx.P, dir, &hit, rwalk); - if(res != RES_OK) goto error; } } diff --git a/src/sdis_heat_path_convective_Xd.h b/src/sdis_heat_path_convective_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_heat_path_radiative_Xd.h b/src/sdis_heat_path_radiative_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -13,9 +13,11 @@ * 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 "sdis_brdf.h" #include "sdis_device_c.h" #include "sdis_green.h" #include "sdis_heat_path.h" +#include "sdis_heat_path_boundary_c.h" #include "sdis_interface_c.h" #include "sdis_log.h" #include "sdis_medium_c.h" @@ -28,7 +30,7 @@ #include "sdis_Xd_begin.h" /******************************************************************************* - * Non generic local functions + * Non generic helper functions ******************************************************************************/ #ifndef SDIS_HEAT_PATH_RADIATIVE_XD_H #define SDIS_HEAT_PATH_RADIATIVE_XD_H @@ -39,7 +41,7 @@ set_limit_radiative_temperature struct rwalk_context* ctx, struct rwalk* rwalk, /* Direction along which the random walk reached the radiative environment */ - const float dir[3], + const double dir[3], const int branch_id, struct temperature* T) { @@ -52,9 +54,10 @@ set_limit_radiative_temperature ASSERT(SXD_HIT_NONE(&rwalk->XD(hit))); rwalk->hit_side = SDIS_SIDE_NULL__; - d3_set_f3(rwalk->dir, dir); + d3_set(rwalk->dir, dir); d3_normalize(rwalk->dir, rwalk->dir); d3_set(ray.dir, rwalk->dir); + ray.time = rwalk->vtx.time; trad = radiative_env_get_temperature(scn->radenv, &ray); if(SDIS_TEMPERATURE_IS_UNKNOWN(trad)) { @@ -85,7 +88,7 @@ set_limit_radiative_temperature * distance of 0.1 meters from the surface along the direction reaching the * radiative environment */ if(ctx->heat_path) { - const float empirical_dst = 0.1f * (float)scn->fp_to_meter; + const double empirical_dst = 0.1 * (float)scn->fp_to_meter; struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; vtx = rwalk->vtx; @@ -103,9 +106,91 @@ error: goto exit; } +/* Check that the trajectory reaches a valid interface, i.e. that it is on a + * fluid/solid interface and has reached it from the fluid */ +static res_T +check_interface + (const struct sdis_interface* interf, + const struct sdis_interface_fragment* frag, + const int verbose) /* Control the verbosity of the function */ +{ + enum sdis_medium_type mdm_frt_type = SDIS_MEDIUM_TYPES_COUNT__; + enum sdis_medium_type mdm_bck_type = SDIS_MEDIUM_TYPES_COUNT__; + enum sdis_side fluid_side = SDIS_SIDE_NULL__; + res_T res = RES_OK; + + mdm_frt_type = sdis_medium_get_type(interf->medium_front); + mdm_bck_type = sdis_medium_get_type(interf->medium_back); + + /* Semi-transparent materials are not supported. This means that a solid/solid + * interface must not be intersected when tracing radiative paths */ + if(mdm_frt_type == SDIS_SOLID && mdm_bck_type == SDIS_SOLID) { + if(verbose) { + log_err(interf->dev, + "Error when sampling the radiatve path. The trajectory reaches a " + "solid/solid interface, whereas this is supposed to be impossible " + "(path position: %g, %g, %g).\n", + SPLIT3(frag->P)); + } + res = RES_BAD_OP; + goto error; + } + + /* Find out which side of the interface the fluid is on */ + if(mdm_frt_type == SDIS_FLUID) { + fluid_side = SDIS_FRONT; + } else if(mdm_bck_type == SDIS_FLUID) { + fluid_side = SDIS_BACK; + } else { + FATAL("Unreachable code\n"); + } + + /* Check that the current position is on the correct side of the interface */ + if(frag->side != fluid_side) { + if(verbose) { + log_err(interf->dev, + "Inconsistent intersection when sampling the radiative path. " + "The path reaches an interface on its solid side, whereas this is " + "supposed to be impossible (path position: %g, %g, %g).\n", + SPLIT3(frag->P)); + } + res = RES_BAD_OP; + goto error; + } + +exit: + return res; +error: + goto exit; +} + #endif /* SDIS_HEAT_PATH_RADIATIVE_XD_H */ /******************************************************************************* + * Generic helper functions + ******************************************************************************/ +static INLINE void +XD(setup_fragment) + (struct sdis_interface_fragment* frag, + const double pos[DIM], + const double dir[DIM], /* Direction _toward_ the hit position */ + const double time, /* Current time */ + const double N[DIM],/* Surface normal */ + const struct sXd(hit)* hit) +{ + struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; + enum sdis_side side = SDIS_SIDE_NULL__; + ASSERT(frag && pos && dir && N); + ASSERT(dX(is_normalized)(N)); + + /* Setup the interface fragment at the intersection position */ + dX(set)(vtx.P, pos); + vtx.time = time; + side = dX(dot)(dir, N) < 0 ? SDIS_FRONT : SDIS_BACK; + XD(setup_interface_fragment)(frag, &vtx, hit, side); +} + +/******************************************************************************* * Local functions ******************************************************************************/ res_T @@ -119,46 +204,40 @@ XD(trace_radiative_path) { /* The radiative random walk is always performed in 3D. In 2D, the geometry * are assumed to be extruded to the infinity along the Z dimension. */ - float N[3] = {0, 0, 0}; - float dir[3] = {0, 0, 0}; + double N[3] = {0,0,0}; + double dir[3] = {0,0,0}; + double pos[3] = {0,0,0}; int branch_id; + size_t nbounces = 0; /* For debug */ res_T res = RES_OK; ASSERT(scn && ray_dir && ctx && rwalk && rng && T); - f3_set(dir, ray_dir); + d3_set_f3(dir, ray_dir); + d3_normalize(dir, dir); /* (int)ctx->nbranchings < 0 <=> Beginning of the realisation */ branch_id = MMAX((int)ctx->nbranchings, 0); /* Launch the radiative random walk */ for(;;) { - const struct sdis_interface* interf = NULL; - struct sdis_medium* chk_mdm = NULL; - struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL; + /* BRDF */ + struct brdf brdf = BRDF_NULL; + struct brdf_sample bounce = BRDF_SAMPLE_NULL; + struct brdf_setup_args brdf_setup_args = BRDF_SETUP_ARGS_NULL; + + /* Miscellaneous */ struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; - unsigned enc_ids[2] = {ENCLOSURE_ID_NULL, ENCLOSURE_ID_NULL}; - unsigned chk_enc_id = ENCLOSURE_ID_NULL; - double alpha; - double epsilon; - double r; - float pos[DIM]; - const float range[2] = { 0, FLT_MAX }; - - fX_set_dX(pos, rwalk->vtx.P); - - /* Trace the radiative ray */ - filter_data.XD(hit) = rwalk->XD(hit); - filter_data.epsilon = 1.e-6; - filter_data.scn = scn; /* Enable the filtering wrt the enclosure id */ - filter_data.enc_id = rwalk->enc_id; -#if (SDIS_XD_DIMENSION == 2) - SXD(scene_view_trace_ray_3d - (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->XD(hit))); -#else - SXD(scene_view_trace_ray - (scn->sXd(view), pos, dir, range, &filter_data, &rwalk->XD(hit))); -#endif + struct sdis_interface* interf = NULL; + struct sdis_medium* chk_mdm = NULL; + double wi[3] = {0,0,0}; + + d3_set(pos, rwalk->vtx.P); + d3_minus(wi, dir); + + res = XD(find_next_fragment)(scn, pos, dir, &rwalk->XD(hit), + rwalk->vtx.time, rwalk->enc_id, &rwalk->XD(hit), &interf, &frag); + if(res != RES_OK) goto error; /* The path reaches the radiative environment */ if(SXD_HIT_NONE(&rwalk->XD(hit))) { @@ -169,60 +248,18 @@ XD(trace_radiative_path) break; /* Stop the radiative path */ } - /* Define the hit side */ - rwalk->hit_side = fX(dot)(dir, rwalk->XD(hit).normal) < 0 - ? SDIS_FRONT : SDIS_BACK; - - /* Move the random walk to the hit position */ - XD(move_pos)(rwalk->vtx.P, dir, rwalk->XD(hit).distance); - - /* Register the random walk vertex against the heat path */ - res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, - SDIS_HEAT_VERTEX_RADIATIVE, branch_id); - if(res != RES_OK) goto error; - - /* Fetch the new interface and setup the hit fragment */ - interf = scene_get_interface(scn, rwalk->XD(hit).prim.prim_id); - XD(setup_interface_fragment)(&frag, &rwalk->vtx, &rwalk->XD(hit), rwalk->hit_side); - - /* Fetch the interface emissivity */ - epsilon = interface_side_get_emissivity(interf, SDIS_INTERN_SOURCE_ID, &frag); - if(epsilon > 1 || epsilon < 0) { - log_err(scn->dev, - "%s: invalid overall emissivity `%g' at position `%g %g %g'.\n", - FUNC_NAME, epsilon, SPLIT3(rwalk->vtx.P)); - res = RES_BAD_OP; - goto error; - } - - /* Switch in boundary temperature ? */ - r = ssp_rng_canonical(rng); - if(r < epsilon) { - T->func = XD(boundary_path); - rwalk->enc_id = ENCLOSURE_ID_NULL; /* Interface between 2 enclosures */ - break; - } - - /* Normalize the normal of the interface and ensure that it points toward the - * current medium */ - fX(normalize)(N, rwalk->XD(hit).normal); - if(rwalk->hit_side == SDIS_BACK) fX(minus)(N, N); - - /* Check that the radiative path is still within the same enclosure. Note - * that this may not be the case, even if the filtering of intersections - * relative to the current enclosure is enabled. This filtering is only - * performed for intersections on a boundary between primitives. As a - * consequence, a threshold effect on how "intersections on a boundary" are - * detected could lead to this situation */ - scene_get_enclosure_ids(scn, rwalk->XD(hit).prim.prim_id, enc_ids); - chk_enc_id = rwalk->hit_side == SDIS_FRONT ? enc_ids[0] : enc_ids[1]; - if(chk_enc_id != rwalk->enc_id) { - log_warn(scn->dev, - "%s: the radiative path has escaped from its cavity -- pos=(%g, %g, %g)\n", - FUNC_NAME, SPLIT3(rwalk->vtx.P)); - res = RES_BAD_OP; - goto error; - } + /* Move the random walk to the hit position, i.e., the next position on the + * interface returned as a fragment by the find_next_fragment function. Do + * not use the sampled direction and distance to the hit point to + * calculate the new position, as the current position may have been slightly + * shifted on the starting triangle by the find_next_fragment function in + * order to avoid numerical inaccuracy issues, making it impossible to + * reconstruct the position actually returned by the function. The starting + * point and distance returned are not, in any case, those used by the + * function to calculate the new wall position. So simply use the position + * returned by this function. */ + d3_set(rwalk->vtx.P, frag.P); + rwalk->hit_side = frag.side; /* Verify that the intersection, although in the same enclosure, touches the * interface of a fluid. We verify this by interface, since a radiative path @@ -235,7 +272,8 @@ XD(trace_radiative_path) * when semi-transparent solids are not yet supported by Stardis. This error * is therefore fatal for the calculation */ chk_mdm = rwalk->hit_side == SDIS_FRONT - ? interf->medium_front : interf->medium_back; + ? interf->medium_front + : interf->medium_back; if(sdis_medium_get_type(chk_mdm) == SDIS_SOLID) { log_err(scn->dev, "%s: a radiative path cannot evolve in a solid -- pos=(%g, %g, %g)\n", @@ -244,13 +282,36 @@ XD(trace_radiative_path) goto error; } - alpha = interface_side_get_specular_fraction(interf, SDIS_INTERN_SOURCE_ID, &frag); - r = ssp_rng_canonical(rng); - if(r < alpha) { /* Sample specular part */ - reflect_3d(dir, f3_minus(dir, dir), N); - } else { /* Sample diffuse part */ - ssp_ran_hemisphere_cos_float(rng, N, dir, NULL); + /* Register the random walk vertex against the heat path */ + res = register_heat_vertex(ctx->heat_path, &rwalk->vtx, T->value, + SDIS_HEAT_VERTEX_RADIATIVE, branch_id); + if(res != RES_OK) goto error; + + /* Retrieve BRDF at current interface position */ + brdf_setup_args.interf = interf; + brdf_setup_args.frag = &frag; + brdf_setup_args.source_id = SDIS_INTERN_SOURCE_ID; + res = brdf_setup(scn->dev, &brdf_setup_args, &brdf); + if(res != RES_OK) goto error; + + /* Switch in boundary temperature? */ + if(ssp_rng_canonical(rng) < brdf.emissivity) { + T->func = XD(boundary_path); + rwalk->enc_id = ENCLOSURE_ID_NULL; /* Interface between 2 enclosures */ + break; } + + /* Normalize the normal of the interface and ensure that it points toward the + * current medium */ + switch(frag.side) { + case SDIS_FRONT: d3_set(N, frag.Ng); break; + case SDIS_BACK: d3_minus(N, frag.Ng); break; + default: FATAL("Unreachable code\n"); break; + } + brdf_sample(&brdf, rng, wi, N, &bounce); + d3_set(dir, bounce.dir); /* Always in 3D */ + + ++nbounces; } exit: @@ -296,4 +357,133 @@ error: goto exit; } +void +XD(trace_ray) + (struct sdis_scene* scn, + const double pos[DIM], + const double dir[3], + const double distance, + const unsigned enc_id, + const struct sXd(hit)* hit_from, + struct sXd(hit)* hit) +{ + struct hit_filter_data filter_data = HIT_FILTER_DATA_NULL; + float ray_org[DIM] = {0}; + float ray_dir[3] = {0}; + float ray_range[2] = {0}; + ASSERT(scn && pos && dir && distance >= 0 && hit_from && hit); + + fX_set_dX(ray_org, pos); + f3_set_d3(ray_dir, dir); + ray_range[0] = 0; + ray_range[1] = (float)distance; + filter_data.XD(hit) = *hit_from; + filter_data.epsilon = 1.e-6; + filter_data.scn = scn; /* Enable the filtering wrt the enclosure id */ + filter_data.enc_id = enc_id; +#if DIM == 2 + SXD(scene_view_trace_ray_3d + (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit)); +#else + SXD(scene_view_trace_ray + (scn->sXd(view), ray_org, ray_dir, ray_range, &filter_data, hit)); +#endif +} + +res_T +XD(find_next_fragment) + (struct sdis_scene* scn, + const double in_pos[DIM], + const double in_dir[3], /* Always in 3D */ + const struct sXd(hit)* in_hit, + const double time, + const unsigned enc_id, + struct sXd(hit)* out_hit, + struct sdis_interface** out_interf, + struct sdis_interface_fragment* out_frag) +{ + int NATTEMPTS_MAX = 10; + int nattempts = 1; + + /* Stardis */ + struct sdis_interface_fragment frag = SDIS_INTERFACE_FRAGMENT_NULL; + struct sdis_interface* interf = NULL; + + struct sXd(hit) hit = SXD_HIT_NULL; + double rt_pos[DIM] = {0}; + res_T res = RES_OK; + + ASSERT(scn && in_pos && in_dir && in_hit); + ASSERT(out_hit && out_interf && out_frag); + + /* Only one attempt is allowed when the ray does not start from a primitive */ + NATTEMPTS_MAX = S3D_HIT_NONE(in_hit) ? 1 : 10; + + dX(set)(rt_pos, in_pos); + + do { + struct sdis_rwalk_vertex vtx = SDIS_RWALK_VERTEX_NULL; + struct sdis_medium* solid = NULL; + double pos[3] = {0}; + double vec[3] = {0}; + double N[3] = {0}; + double delta = 0; + + /* Reset result code. It may have been modified during a previous attempt */ + res = RES_OK; + + /* Find the following surface along the direction of propagation */ + XD(trace_ray)(scn, rt_pos, in_dir, INF, enc_id, in_hit, &hit); + if(SXD_HIT_NONE(&hit)) break; + + /* Retrieve the current position and normal */ + dX(add)(pos, rt_pos, dX(muld)(vec, in_dir, hit.distance)); + dX_set_fX(N, hit.normal); + dX(normalize(N, N)); + + /* Retrieve the current interface properties */ + interf = scene_get_interface(scn, hit.prim.prim_id); + XD(setup_fragment)(&frag, pos, in_dir, time, N, &hit); + + /* Check that the path reaches a valid interface. + * An invalid fragment may mean that the ray position is in a corner and the + * traced ray has missed the surface of that corner. To correct this, the + * ray position is moved slightly away from the corner before a ray is drawn + * in the same direction. This fallback solution is executed a number of + * times, after which, if the fragment is still invalid, it is considered + * that the numerical error cannot be mitigated. */ + res = check_interface(interf, &frag, nattempts == NATTEMPTS_MAX); + if(res != RES_OK && nattempts == NATTEMPTS_MAX) goto error; + ++nattempts; + + if(res != RES_OK) { /* Mitigate numerical error (see above) */ + if(sdis_medium_get_type(interf->medium_front) == SDIS_SOLID) { + solid = interf->medium_front; + } else { + ASSERT(sdis_medium_get_type(interf->medium_back) == SDIS_SOLID); + solid = interf->medium_back; + } + + /* Retrieves the delta of the solid that surrounds the boundary, as it is + * actually the only numerical parameter that says something about the + * system. */ + vtx.P[0] = pos[0]; + vtx.P[1] = pos[1]; + vtx.P[2] = pos[2]; + vtx.time = time; + delta = solid_get_delta(solid, &vtx); + + XD(move_away_primitive_boundaries)(in_hit, delta, rt_pos); + } + } while(res != RES_OK); + +exit: + *out_hit = hit; + *out_interf = interf; + *out_frag = frag; + return res; +error: + goto exit; +} + #include "sdis_Xd_end.h" diff --git a/src/sdis_interface.c b/src/sdis_interface.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_interface_c.h b/src/sdis_interface_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_log.c b/src/sdis_log.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_log.h b/src/sdis_log.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_medium.c b/src/sdis_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_medium_c.h b/src/sdis_medium_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_misc.c b/src/sdis_misc.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_misc.h b/src/sdis_misc.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_mpi.c b/src/sdis_mpi.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_mpi.h b/src/sdis_mpi.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_primkey.c b/src/sdis_primkey.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_radiative_env.c b/src/sdis_radiative_env.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_radiative_env_c.h b/src/sdis_radiative_env_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_realisation.c b/src/sdis_realisation.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_realisation.h b/src/sdis_realisation.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_realisation_Xd.h b/src/sdis_realisation_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_scene.c b/src/sdis_scene.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_scene_Xd.h b/src/sdis_scene_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_scene_c.h b/src/sdis_scene_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_solve.c b/src/sdis_solve.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_solve_boundary_Xd.h b/src/sdis_solve_boundary_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_solve_camera.c b/src/sdis_solve_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_solve_medium_Xd.h b/src/sdis_solve_medium_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_solve_probe_Xd.h b/src/sdis_solve_probe_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_solve_probe_boundary_Xd.h b/src/sdis_solve_probe_boundary_Xd.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_source.c b/src/sdis_source.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -163,43 +163,25 @@ sdis_source_get_id(const struct sdis_source* source) res_T source_sample (const struct sdis_source* src, + const struct source_props* props, struct ssp_rng* rng, const double pos[3], - const double time, struct source_sample* sample) { - double src_pos[3]; /* [m] */ double main_dir[3]; double half_angle; /* [radians] */ double cos_half_angle; /* [radians] */ double dst; /* [m] */ - double radius; /* Source radius [m] */ - double power; /* Source power [W] */ - double area; /* Source area [m^2] */ res_T res = RES_OK; ASSERT(src && rng && pos && sample); - /* Retrieve current source position, radius and power */ - src->spherical.position(time, src_pos, src->data); - power = src->spherical.power(time, src->data); - radius = src->spherical.radius; - - if(power < 0) { - log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n", - FUNC_NAME, power); - res = RES_BAD_ARG; - goto error; - } - - area = 4*PI*radius*radius; /* [m^2] */ - /* compute the direction of `pos' toward the center of the source */ - d3_sub(main_dir, src_pos, pos); + d3_sub(main_dir, props->pos, pos); /* Normalize the direction and keep the distance from `pos' to the center of * the source */ dst = d3_normalize(main_dir, main_dir); - if(dst <= radius) { + if(dst <= props->radius) { log_err(src->dev, "%s: the position from which the external source is sampled " "is included in the source:\n" @@ -208,34 +190,32 @@ source_sample "\tposition = %g, %g, %g\n" "\ttime = %g\n" "\tdistance from position to source = %g\n", - FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst); + FUNC_NAME, SPLIT3(props->pos), props->radius, SPLIT3(pos), props->time, dst); res = RES_BAD_ARG; goto error; } /* Point source */ - if(area == 0) { + if(props->area == 0) { d3_set(sample->dir, main_dir); sample->pdf = 1; sample->dst = dst; - sample->power = power; /* [W] */ sample->radiance_term = 1.0 / (4*PI*dst*dst); /* [W/m^2/sr] */ - sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */ + sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */ /* Spherical source */ } else { /* Sample the source according to its solid angle, * i.e. 2*PI*(1 - cos(half_angle)) */ - half_angle = asin(radius/dst); + half_angle = asin(props->radius/dst); cos_half_angle = cos(half_angle); ssp_ran_sphere_cap_uniform /* pdf = 1/(2*PI*(1-cos(half_angle))) */ (rng, main_dir, cos_half_angle, sample->dir, &sample->pdf); /* Set other sample variables */ - sample->dst = dst - radius; /* From pos to source boundaries [m] */ - sample->power = power; /* [W] */ - sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */ - sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */ + sample->dst = dst - props->radius; /* From pos to source boundaries [m] */ + sample->radiance_term = 1.0 / (PI*props->area); /* [W/m^2/sr] */ + sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */ } exit: @@ -247,47 +227,31 @@ error: res_T source_trace_to (const struct sdis_source* src, + const struct source_props* props, const double pos[3], /* Ray origin */ const double dir[3], /* Ray direction */ - const double time, /* Time at which ray is traced */ struct source_sample* sample) { - double src_pos[3]; /* [m] */ double main_dir[3]; - double radius; /* [m] */ - double power; /* [W] */ double dst; /* Distance from pos to the source center [m] */ double half_angle; /* [radian] */ res_T res = RES_OK; - ASSERT(src && pos && dir && sample); + ASSERT(src && props && pos && dir && sample); ASSERT(d3_is_normalized(dir)); - radius = src->spherical.radius; - /* Point sources cannot be targeted */ - if(radius == 0) { + if(props->radius == 0) { *sample = SOURCE_SAMPLE_NULL; goto exit; } - /* Retrieve current source position and power */ - src->spherical.position(time, src_pos, src->data); - power = src->spherical.power(time, src->data); - - if(power < 0) { - log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n", - FUNC_NAME, power); - res = RES_BAD_ARG; - goto error; - } - /* compute the direction of `pos' toward the center of the source */ - d3_sub(main_dir, src_pos, pos); + d3_sub(main_dir, props->pos, pos); /* Normalize the direction and keep the distance from `pos' to the center of * the source */ dst = d3_normalize(main_dir, main_dir); - if(dst <= radius) { + if(dst <= props->radius) { log_err(src->dev, "%s: the position from which the external source is targeted " "is included in the source:\n" @@ -296,13 +260,13 @@ source_trace_to "\tposition = %g, %g, %g\n" "\ttime = %g\n" "\tdistance from position to source = %g\n", - FUNC_NAME, SPLIT3(src_pos), radius, SPLIT3(pos), time, dst); + FUNC_NAME, SPLIT3(props->pos), props->radius, SPLIT3(pos), props->time, dst); res = RES_BAD_ARG; goto error; } /* Compute the half angle of the source as seen from pos */ - half_angle = asin(radius/dst); + half_angle = asin(props->radius/dst); /* The source is missed */ if(d3_dot(dir, main_dir) < cos(half_angle)) { @@ -310,14 +274,11 @@ source_trace_to /* The source is intersected */ } else { - const double area = 4*PI*radius*radius; /* [m^2] */ - d3_set(sample->dir, dir); sample->pdf = 1; - sample->dst = dst - radius; /* From pos to source boundaries [m] */ - sample->power = power; /* [W] */ - sample->radiance_term = 1.0 / (PI*area); /* [W/m^2/sr] */ - sample->radiance = sample->power * sample->radiance_term; /* [W/m^2/sr] */ + sample->dst = dst - props->radius; /* From pos to source boundaries [m] */ + sample->radiance_term = 1.0 / (PI*props->area); /* [W/m^2/sr] */ + sample->radiance = props->power * sample->radiance_term; /* [W/m^2/sr] */ } exit: @@ -327,6 +288,36 @@ error: goto exit; } +res_T +source_get_props + (const struct sdis_source* src, + const double time, /* [s] */ + struct source_props* props) +{ + res_T res = RES_OK; + ASSERT(src && props); + + /* Retrieve the source properties */ + src->spherical.position(time, props->pos, src->data); + props->power = src->spherical.power(time, src->data); + props->radius = src->spherical.radius; + + if(props->power < 0) { + log_err(src->dev, "%s: invalid source power '%g' W. It cannot be negative.\n", + FUNC_NAME, props->power); + res = RES_BAD_ARG; + goto error; + } + + props->area = 4*PI*props->radius*props->radius; /* [m^2] */ + props->time = time; /* [s] */ + +exit: + return res; +error: + goto exit; +} + double /* [W] */ source_get_power(const struct sdis_source* src, const double time /* [s] */) { diff --git a/src/sdis_source_c.h b/src/sdis_source_c.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -22,11 +22,21 @@ struct sdis_source; struct ssp_rng; +struct source_props { + double pos[3]; /* [m/fp_to_meter] */ + double radius; /* [m/fp_to_meter] */ + double power; /* [W] */ + double area; /* [m^2/fp_to_meter] */ + double time; /* [s] */ +}; +#define SOURCE_PROPS_NULL__ {0} +static const struct source_props SOURCE_PROPS_NULL = SOURCE_PROPS_NULL__; + struct source_sample { double dir[3]; /* Direction _to_ the source */ double pdf; /* pdf of sampled direction */ double dst; /* Distance to the source [m] */ - double power; /* [W] */ + double radiance; /* [W/m^2/sr] */ /* Radiance relative to power, i.e. the source power is assumed to be equal to @@ -36,18 +46,24 @@ struct source_sample { * green function */ double radiance_term; /* [W/m^2/sr] */ }; -#define SOURCE_SAMPLE_NULL__ {{0,0,0}, 0, 0, 0, 0, 0} +#define SOURCE_SAMPLE_NULL__ {0} static const struct source_sample SOURCE_SAMPLE_NULL = SOURCE_SAMPLE_NULL__; /* Helper macro used to define whether a sample is valid or not */ #define SOURCE_SAMPLE_NONE(Sample) ((Sample)->pdf == 0) extern LOCAL_SYM res_T +source_get_props + (const struct sdis_source* source, + const double time, /* Time at which props are retrieved [s] */ + struct source_props* props); + +extern LOCAL_SYM res_T source_sample (const struct sdis_source* source, + const struct source_props* props, struct ssp_rng* rng, const double pos[3], /* Position from which the source is sampled */ - const double time, /* Time at which the source is sampled */ struct source_sample* sample); /* Trace a ray toward the source. The returned sample has a pdf of 1 or 0 @@ -56,9 +72,9 @@ source_sample extern LOCAL_SYM res_T source_trace_to (const struct sdis_source* source, + const struct source_props* props, const double pos[3], /* Ray origin */ const double dir[3], /* Ray direction */ - const double time, /* Time at which ray is traced */ struct source_sample* sample); /* pdf == 0 if no source is reached */ extern LOCAL_SYM double /* [W] */ diff --git a/src/sdis_tile.c b/src/sdis_tile.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_tile.h b/src/sdis_tile.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/sdis_version.h.in b/src/sdis_version.h.in @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis.c b/src/test_sdis.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_accum_buffer.c b/src/test_sdis_accum_buffer.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_camera.c b/src/test_sdis_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_compute_power.c b/src/test_sdis_compute_power.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_conducto_radiative.c b/src/test_sdis_conducto_radiative.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_conducto_radiative_2d.c b/src/test_sdis_conducto_radiative_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_contact_resistance.c b/src/test_sdis_contact_resistance.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_contact_resistance.h b/src/test_sdis_contact_resistance.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_contact_resistance_2.c b/src/test_sdis_contact_resistance_2.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_convection.c b/src/test_sdis_convection.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_convection_non_uniform.c b/src/test_sdis_convection_non_uniform.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_custom_solid_path_sampling.c b/src/test_sdis_custom_solid_path_sampling.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_custom_solid_path_sampling_2d.c b/src/test_sdis_custom_solid_path_sampling_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_data.c b/src/test_sdis_data.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_device.c b/src/test_sdis_device.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_draw_external_flux.c b/src/test_sdis_draw_external_flux.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_enclosure_limit_conditions.c b/src/test_sdis_enclosure_limit_conditions.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_external_flux.c b/src/test_sdis_external_flux.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_external_flux_with_diffuse_radiance.c b/src/test_sdis_external_flux_with_diffuse_radiance.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_flux.c b/src/test_sdis_flux.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_flux2.c b/src/test_sdis_flux2.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_flux_with_h.c b/src/test_sdis_flux_with_h.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_interface.c b/src/test_sdis_interface.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_medium.c b/src/test_sdis_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_mesh.h b/src/test_sdis_mesh.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_picard.c b/src/test_sdis_picard.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -225,7 +225,7 @@ solid_get_delta (const struct sdis_rwalk_vertex* vtx, struct sdis_data* data) { CHK(vtx && data); - return 0.005; + return 0.0025; } static double @@ -435,6 +435,7 @@ static void test_picard (struct sdis_scene* scn, const size_t picard_order, + const enum sdis_diffusion_algorithm algo, const struct reference_result* ref) { struct sdis_solve_probe_args probe_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; @@ -458,6 +459,7 @@ test_picard probe_args.position[1] = 0; probe_args.position[2] = 0; probe_args.picard_order = picard_order; + probe_args.diff_algo = algo; OK(sdis_solve_probe(scn, &probe_args, &estimator)); OK(sdis_estimator_get_temperature(estimator, &mc)); printf("Temperature at `%g %g %g' = %g ~ %g +/- %g\n", @@ -687,7 +689,7 @@ main(int argc, char** argv) interf_props.temperature = 350; interf_props.h = -1; interf_props.emissivity = 1; - interf_props.specular_fraction = -1; + interf_props.specular_fraction = 0; interf_props.Tref = 350; create_interface(dev, fluid, dummy, &interf_props, interfaces+BOUNDARY_pX); @@ -712,8 +714,8 @@ main(int argc, char** argv) pinterf_props[BOUNDARY_pX]->Tref = 300; radenv_props->temperature = 280; radenv_props->reference = 300; - test_picard(scn_2d, 1/*Picard order*/, &ref); - test_picard(scn_3d, 1/*Picard order*/, &ref); + test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); + test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); printf("\n"); /* Test picard1 using T4 as a reference */ @@ -726,8 +728,8 @@ main(int argc, char** argv) pinterf_props[BOUNDARY_pX]->Tref = 350; radenv_props->temperature = 280; radenv_props->reference = 280; - test_picard(scn_2d, 1/*Picard order*/, &ref); - test_picard(scn_3d, 1/*Picard order*/, &ref); + test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); + test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); printf("\n"); /* Test picard2 */ @@ -740,8 +742,8 @@ main(int argc, char** argv) pinterf_props[BOUNDARY_pX]->Tref = 300; radenv_props->temperature = 280; radenv_props->reference = 300; - test_picard(scn_2d, 2/*Picard order*/, &ref); - test_picard(scn_3d, 2/*Picard order*/, &ref); + test_picard(scn_2d, 2/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); + test_picard(scn_3d, 2/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); printf("\n"); t_range[0] = 200; @@ -761,8 +763,8 @@ main(int argc, char** argv) pinterf_props[BOUNDARY_pX]->Tref = pinterf_props[BOUNDARY_pX]->temperature; radenv_props->temperature = t_range[0]; radenv_props->reference = t_range[0]; - test_picard(scn_2d, 3/*Picard order*/, &ref); - test_picard(scn_3d, 3/*Picard order*/, &ref); + test_picard(scn_2d, 3/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); + test_picard(scn_3d, 3/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); register_heat_paths(scn_2d, 3/*Picard order*/, stream); register_heat_paths(scn_3d, 3/*Picard order*/, stream); printf("\n"); @@ -787,8 +789,8 @@ main(int argc, char** argv) pinterf_props[BOUNDARY_pX]->Tref = 300; radenv_props->temperature = t_range[0]; radenv_props->reference = 300; - test_picard(scn_2d, 1/*Picard order*/, &ref); - test_picard(scn_3d, 1/*Picard order*/, &ref); + test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); + test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_DELTA_SPHERE, &ref); printf("\n"); /* Test picard1 with a volumic power and T4 a the reference */ @@ -801,8 +803,8 @@ main(int argc, char** argv) pinterf_props[BOUNDARY_pX]->Tref = 350; radenv_props->temperature = 280; radenv_props->reference = 280; - test_picard(scn_2d, 1/*Picard order*/, &ref); - test_picard(scn_3d, 1/*Picard order*/, &ref); + test_picard(scn_2d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); + test_picard(scn_3d, 1/*Picard order*/, SDIS_DIFFUSION_WOS, &ref); printf("\n"); /* Release memory */ diff --git a/src/test_sdis_primkey.c b/src/test_sdis_primkey.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_primkey_2d.c b/src/test_sdis_primkey_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_radiative_env.c b/src/test_sdis_radiative_env.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_scene.c b/src/test_sdis_scene.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solid_random_walk_robustness.c b/src/test_sdis_solid_random_walk_robustness.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -28,29 +28,6 @@ /******************************************************************************* * Helper functions ******************************************************************************/ -static void -compute_aabb - (const double* positions, - const size_t nvertices, - double lower[3], - double upper[3]) -{ - size_t i; - CHK(positions && nvertices && lower && upper); - - lower[0] = lower[1] = lower[2] = DBL_MAX; - upper[0] = upper[1] = upper[2] =-DBL_MAX; - - FOR_EACH(i, 0, nvertices) { - lower[0] = MMIN(lower[0], positions[i*3+0]); - lower[1] = MMIN(lower[1], positions[i*3+1]); - lower[2] = MMIN(lower[2], positions[i*3+2]); - upper[0] = MMAX(upper[0], positions[i*3+0]); - upper[1] = MMAX(upper[1], positions[i*3+1]); - upper[2] = MMAX(upper[2], positions[i*3+2]); - } -} - static double trilinear_temperature(const double pos[3]) { @@ -81,22 +58,17 @@ volumetric_temperature(const double pos[3], const double upper[3]) return temp; } -static void -print_estimation_result - (const struct sdis_estimator* estimator, const double Tref) +static const char* +algo_cstr(const enum sdis_diffusion_algorithm diff_algo) { - struct sdis_mc T = SDIS_MC_NULL; - size_t nreals; - size_t nfails; + const char* cstr = "none"; - OK(sdis_estimator_get_temperature(estimator, &T)); - OK(sdis_estimator_get_realisation_count(estimator, &nreals)); - OK(sdis_estimator_get_failure_count(estimator, &nfails)); - CHK(nfails <= Nreals * 0.0005); - printf("T = %g ~ %g +/- %g [%g, %g]; #failures = %lu / %lu\n", - Tref, T.E, T.SE, T.E - 3*T.SE, T.E + 3*T.SE, - (unsigned long)nfails, (unsigned long)Nreals); - CHK(eq_eps(T.E, Tref, T.SE*3)); + switch(diff_algo) { + case SDIS_DIFFUSION_DELTA_SPHERE: cstr = "delta sphere"; break; + case SDIS_DIFFUSION_WOS: cstr = "WoS"; break; + default: FATAL("Unreachable code.\n"); break; + } + return cstr; } /******************************************************************************* @@ -135,13 +107,55 @@ get_interface(const size_t itri, struct sdis_interface** bound, void* context) *bound = ctx->interf; } +static struct sdis_scene* +create_scene(struct sdis_device* dev, struct sdis_interface* interf) +{ + /* Star-3DUT */ + struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; + struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; + struct s3dut_mesh* msh = NULL; + + /* Stardis */ + struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; + struct sdis_scene* scn = NULL; + + /* Miscellaneous */ + struct context ctx; + + ASSERT(dev && interf); + + /* Create the solid super shape */ + f0.A = 1; f0.B = 1; f0.M = 20; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5; + f1.A = 1; f1.B = 1; f1.M = 7; f1.N0 = 1; f1.N1 = 2; f1.N2 = 5; + OK(s3dut_create_super_shape(NULL, &f0, &f1, 1, 128, 64, &msh)); + OK(s3dut_mesh_get_data(msh, &ctx.msh)); + + /*dump_mesh(stdout, ctx.msh.positions, + ctx.msh.nvertices, ctx.msh.indices, ctx.msh.nprimitives);*/ + + /* Create the scene */ + ctx.interf = interf; + scn_args.get_indices = get_indices; + scn_args.get_interface = get_interface; + scn_args.get_position = get_position; + scn_args.nprimitives = ctx.msh.nprimitives; + scn_args.nvertices = ctx.msh.nvertices; + scn_args.context = &ctx; + OK(sdis_scene_create(dev, &scn_args, &scn)); + + OK(s3dut_mesh_ref_put(msh)); + + return scn; +} + /******************************************************************************* * Interface ******************************************************************************/ enum profile { PROFILE_UNKNOWN, PROFILE_TRILINEAR, - PROFILE_VOLUMETRIC_POWER + PROFILE_VOLUMETRIC_POWER, + PROFILE_COUNT__ }; struct interf { enum profile profile; @@ -180,6 +194,31 @@ interface_get_convection_coef return ((const struct interf*)sdis_data_cget(data))->h;; } +static struct sdis_interface* +create_interface + (struct sdis_device* dev, + struct sdis_medium* front, + struct sdis_medium* back) +{ + struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; + struct sdis_interface* interf = NULL; + struct sdis_data* data = NULL; + struct interf* interf_param = NULL; + ASSERT(dev && front && back); + + OK(sdis_data_create + (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); + interf_param = sdis_data_get(data); + interf_param->h = 1; + + interf_shader.convection_coef = interface_get_convection_coef; + interf_shader.front.temperature = interface_get_temperature; + OK(sdis_interface_create(dev, front, back, &interf_shader, data, &interf)); + OK(sdis_data_ref_put(data)); + + return interf; +} + /******************************************************************************* * Fluid ******************************************************************************/ @@ -191,6 +230,19 @@ fluid_get_temperature (void)data; return Tfluid; } +static struct sdis_medium* +create_fluid(struct sdis_device* dev) +{ + struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; + struct sdis_medium* fluid = NULL; + ASSERT(dev); + + /* Create the fluid medium */ + fluid_shader.temperature = fluid_get_temperature; + OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); + + return fluid; +} /******************************************************************************* * Solid API @@ -252,51 +304,15 @@ solid_get_volumetric_power return ((const struct solid*)sdis_data_cget(data))->power; } -/******************************************************************************* - * Main function - ******************************************************************************/ -int -main(int argc, char** argv) +static struct sdis_medium* +create_solid(struct sdis_device* dev) { - struct s3dut_super_formula f0 = S3DUT_SUPER_FORMULA_NULL; - struct s3dut_super_formula f1 = S3DUT_SUPER_FORMULA_NULL; - struct s3dut_mesh* msh = NULL; - struct sdis_device* dev = NULL; - struct sdis_data* data = NULL; - struct sdis_estimator* estimator = NULL; - struct sdis_medium* fluid = NULL; - struct sdis_medium* solid = NULL; - struct sdis_interface* interf = NULL; - struct sdis_scene* scn = NULL; - struct sdis_scene_create_args scn_args = SDIS_SCENE_CREATE_ARGS_DEFAULT; - struct sdis_fluid_shader fluid_shader = DUMMY_FLUID_SHADER; struct sdis_solid_shader solid_shader = SDIS_SOLID_SHADER_NULL; - struct sdis_interface_shader interf_shader = SDIS_INTERFACE_SHADER_NULL; - struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; - struct sdis_solve_medium_args solve_mdm_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; - struct interf* interf_param = NULL; - struct solid* solid_param = NULL; - struct context ctx; - double lower[3]; - double upper[3]; - double spread; - (void)argc, (void)argv; - - OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); - - /* Create the fluid medium */ - fluid_shader.temperature = fluid_get_temperature; - OK(sdis_fluid_create(dev, &fluid_shader, NULL, &fluid)); - - /* Setup the solid shader */ - solid_shader.calorific_capacity = solid_get_calorific_capacity; - solid_shader.thermal_conductivity = solid_get_thermal_conductivity; - solid_shader.volumic_mass = solid_get_volumic_mass; - solid_shader.delta = solid_get_delta; - solid_shader.temperature = solid_get_temperature; - solid_shader.volumic_power = solid_get_volumetric_power; + struct sdis_medium* solid = NULL; + struct sdis_data* data = NULL; + struct solid* solid_param; + ASSERT(dev); - /* Create the solid medium */ OK(sdis_data_create (dev, sizeof(struct solid), ALIGNOF(struct solid), NULL, &data)); solid_param = sdis_data_get(data); @@ -306,47 +322,53 @@ main(int argc, char** argv) solid_param->rho = 1.0; solid_param->temperature = SDIS_TEMPERATURE_NONE; solid_param->power = SDIS_VOLUMIC_POWER_NONE; - OK(sdis_solid_create(dev, &solid_shader, data, &solid)); - OK(sdis_data_ref_put(data)); - /* Create the interface */ - OK(sdis_data_create - (dev, sizeof(struct interf), ALIGNOF(struct interf), NULL, &data)); - interf_param = sdis_data_get(data); - interf_param->h = 1; - interf_shader.convection_coef = interface_get_convection_coef; - interf_shader.front.temperature = interface_get_temperature; - OK(sdis_interface_create(dev, solid, fluid, &interf_shader, data, &interf)); + solid_shader.calorific_capacity = solid_get_calorific_capacity; + solid_shader.thermal_conductivity = solid_get_thermal_conductivity; + solid_shader.volumic_mass = solid_get_volumic_mass; + solid_shader.delta = solid_get_delta; + solid_shader.temperature = solid_get_temperature; + solid_shader.volumic_power = solid_get_volumetric_power; + OK(sdis_solid_create(dev, &solid_shader, data, &solid)); OK(sdis_data_ref_put(data)); - /* Create the solid super shape */ - f0.A = 1; f0.B = 1; f0.M = 20; f0.N0 = 1; f0.N1 = 1; f0.N2 = 5; - f1.A = 1; f1.B = 1; f1.M = 7; f1.N0 = 1; f1.N1 = 2; f1.N2 = 5; - OK(s3dut_create_super_shape(NULL, &f0, &f1, 1, 128, 64, &msh)); - OK(s3dut_mesh_get_data(msh, &ctx.msh)); + return solid; +} - compute_aabb(ctx.msh.positions, ctx.msh.nvertices, lower, upper); +/******************************************************************************* + * Solve functions + ******************************************************************************/ +static void +check_estimation + (const struct sdis_estimator* estimator, const double Tref) +{ + struct sdis_mc T = SDIS_MC_NULL; + size_t nreals; + size_t nfails; - /* Create the scene */ - ctx.interf = interf; - scn_args.get_indices = get_indices; - scn_args.get_interface = get_interface; - scn_args.get_position = get_position; - scn_args.nprimitives = ctx.msh.nprimitives; - scn_args.nvertices = ctx.msh.nvertices; - scn_args.context = &ctx; - OK(sdis_scene_create(dev, &scn_args, &scn)); - /*dump_mesh(stdout, ctx.msh.positions, - ctx.msh.nvertices, ctx.msh.indices, ctx.msh.nprimitives);*/ + OK(sdis_estimator_get_temperature(estimator, &T)); + OK(sdis_estimator_get_realisation_count(estimator, &nreals)); + OK(sdis_estimator_get_failure_count(estimator, &nfails)); + CHK(nfails <= Nreals * 0.0005); + printf("T = %g ~ %g +/- %g [%g, %g]; #failures = %lu / %lu\n", + Tref, T.E, T.SE, T.E - 3*T.SE, T.E + 3*T.SE, + (unsigned long)nfails, (unsigned long)Nreals); + CHK(eq_eps(T.E, Tref, T.SE*3)); +} - /* Compute the delta of the solid random walk */ - OK(sdis_scene_get_medium_spread(scn, solid, &spread)); - solid_param->delta = 0.4 / spread; /* (4V/S) / 10 */ +static void +check_probe + (struct sdis_scene* scn, + const enum profile profile, + const enum sdis_diffusion_algorithm diff_algo) +{ + struct sdis_solve_probe_args solve_args = SDIS_SOLVE_PROBE_ARGS_DEFAULT; + struct sdis_estimator* estimator = NULL; + double lower[3], upper[3]; + double Tref; + ASSERT(scn); - interf_param->upper[0] = upper[0]; - interf_param->upper[1] = upper[1]; - interf_param->upper[2] = upper[2]; - interf_param->h = 1; + printf("algo: %s\n", algo_cstr(diff_algo)); solve_args.nrealisations = Nreals; solve_args.position[0] = 0; @@ -354,38 +376,116 @@ main(int argc, char** argv) solve_args.position[2] = 0; solve_args.time_range[0] = INF; solve_args.time_range[1] = INF; + solve_args.diff_algo = diff_algo; solve_args.register_paths = SDIS_HEAT_PATH_FAILURE; - - /* Launch probe estimation with trilinear profile set at interfaces */ - interf_param->profile = PROFILE_TRILINEAR; OK(sdis_solve_probe(scn, &solve_args, &estimator)); - print_estimation_result - (estimator, trilinear_temperature(solve_args.position)); - OK(sdis_estimator_ref_put(estimator)); - /* Launch probe estimation with volumetric power profile set at interfaces */ - interf_param->profile = PROFILE_VOLUMETRIC_POWER; - solid_param->power = Pw; - OK(sdis_solve_probe(scn, &solve_args, &estimator)); - print_estimation_result - (estimator, volumetric_temperature(solve_args.position, upper)); - solid_param->power = SDIS_VOLUMIC_POWER_NONE; + switch(profile) { + case PROFILE_TRILINEAR: + Tref = trilinear_temperature(solve_args.position); + break; + case PROFILE_VOLUMETRIC_POWER: + OK(sdis_scene_get_aabb(scn, lower, upper)); + Tref = volumetric_temperature(solve_args.position, upper); + break; + default: FATAL("Unreachable code.\n"); break; + } + + check_estimation(estimator, Tref); OK(sdis_estimator_ref_put(estimator)); +} + +static void +check_medium + (struct sdis_scene* scn, + struct sdis_medium* medium, + const enum sdis_diffusion_algorithm diff_algo) +{ + struct sdis_solve_medium_args solve_mdm_args = SDIS_SOLVE_MEDIUM_ARGS_DEFAULT; + struct sdis_estimator* estimator = NULL; + ASSERT(scn); + + printf("algo: %s\n", algo_cstr(diff_algo)); - /* Launch medium integration */ - interf_param->profile = PROFILE_UNKNOWN; solve_mdm_args.nrealisations = Nreals; - solve_mdm_args.medium = solid; + solve_mdm_args.medium = medium; solve_mdm_args.time_range[0] = INF; solve_mdm_args.time_range[1] = INF; + solve_mdm_args.diff_algo = diff_algo; OK(sdis_solve_medium(scn, &solve_mdm_args, &estimator)); - print_estimation_result(estimator, Tfluid); + check_estimation(estimator, Tfluid); /*dump_heat_paths(stdout, estimator);*/ OK(sdis_estimator_ref_put(estimator)); +} + +/******************************************************************************* + * Main function + ******************************************************************************/ +int +main(int argc, char** argv) +{ + struct sdis_device* dev = NULL; + struct sdis_medium* fluid = NULL; + struct sdis_medium* solid = NULL; + struct sdis_interface* interf = NULL; + struct sdis_scene* scn = NULL; + struct interf* interf_param = NULL; + struct solid* solid_param = NULL; + double lower[3]; + double upper[3]; + double spread; + (void)argc, (void)argv; + + OK(sdis_device_create(&SDIS_DEVICE_CREATE_ARGS_DEFAULT, &dev)); + + fluid = create_fluid(dev); + solid = create_solid(dev); + interf = create_interface(dev, solid, fluid); + scn = create_scene(dev, interf); + + solid_param = sdis_data_get(sdis_medium_get_data(solid)); + interf_param = sdis_data_get(sdis_interface_get_data(interf)); + + OK(sdis_scene_get_medium_spread(scn, solid, &spread)); + OK(sdis_scene_get_aabb(scn, lower, upper)); + + /* Compute the delta of the solid random walk */ + solid_param->delta = 0.4 / spread; + + /* Setup the upper boundary required to solve the trilinear profile */ + interf_param->upper[0] = upper[0]; + interf_param->upper[1] = upper[1]; + interf_param->upper[2] = upper[2]; + + /* Launch probe estimation with trilinear profile set at interfaces */ + solid_param->delta = 0.4 / spread; + interf_param->profile = PROFILE_TRILINEAR; + check_probe(scn, PROFILE_TRILINEAR, SDIS_DIFFUSION_DELTA_SPHERE); + solid_param->delta = 1e-5 / spread; /* Make life difficult for WoS */ + check_probe(scn, PROFILE_TRILINEAR, SDIS_DIFFUSION_WOS); + + /* Launch probe estimation with volumetric power profile set at interfaces */ + solid_param->delta = 0.4 / spread; + solid_param->power = Pw; + interf_param->profile = PROFILE_VOLUMETRIC_POWER; + check_probe(scn, PROFILE_VOLUMETRIC_POWER, SDIS_DIFFUSION_DELTA_SPHERE); + solid_param->delta = 1e-5 / spread; /* Make life difficult for WoS */ + check_probe(scn, PROFILE_VOLUMETRIC_POWER, SDIS_DIFFUSION_WOS); + solid_param->power = SDIS_VOLUMIC_POWER_NONE; + + /* Launch medium integration. Do not use an analytic profile as a boundary + * condition but a Robin boundary condition */ + interf_param->profile = PROFILE_UNKNOWN; + solid_param->delta = 0.4 / spread; + check_medium(scn, solid, SDIS_DIFFUSION_DELTA_SPHERE); + + /* Contrary to previous WoS tests, don't reduce the delta parameter to avoid + * prohibitive increases in calculation time: too small a delta would trap the + * path in the solid */ + check_medium(scn, solid, SDIS_DIFFUSION_WOS); /* Release */ - OK(s3dut_mesh_ref_put(msh)); OK(sdis_device_ref_put(dev)); OK(sdis_medium_ref_put(fluid)); OK(sdis_medium_ref_put(solid)); @@ -395,4 +495,3 @@ main(int argc, char** argv) CHK(mem_allocated_size() == 0); return 0; } - diff --git a/src/test_sdis_solve_boundary.c b/src/test_sdis_solve_boundary.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_boundary_flux.c b/src/test_sdis_solve_boundary_flux.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -123,7 +123,7 @@ solid_get_delta { (void) data; CHK(vtx != NULL); - return 1.0 / 20.0; + return 1.0 / 40.0; } static double diff --git a/src/test_sdis_solve_camera.c b/src/test_sdis_solve_camera.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_medium.c b/src/test_sdis_solve_medium.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_medium_2d.c b/src/test_sdis_solve_medium_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_probe.c b/src/test_sdis_solve_probe.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_probe2.c b/src/test_sdis_solve_probe2.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -296,7 +296,7 @@ main(int argc, char** argv) /* Check the RNG state */ OK(ssp_rng_create(NULL, SSP_RNG_THREEFRY, &rng)); - OK(ssp_rng_discard(rng, 3141592653589)); /* Move the RNG state */ + OK(ssp_rng_discard(rng, 314159265358979)); /* Move the RNG state */ solve_args.rng_state = rng; solve_args.rng_type = SSP_RNG_TYPE_NULL; OK(sdis_solve_probe(scn, &solve_args, &estimator2)); diff --git a/src/test_sdis_solve_probe2_2d.c b/src/test_sdis_solve_probe2_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_probe3.c b/src/test_sdis_solve_probe3.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_probe3_2d.c b/src/test_sdis_solve_probe3_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_probe_2d.c b/src/test_sdis_solve_probe_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_probe_boundary_list.c b/src/test_sdis_solve_probe_boundary_list.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_solve_probe_list.c b/src/test_sdis_solve_probe_list.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_source.c b/src/test_sdis_source.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_unsteady.c b/src/test_sdis_unsteady.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_unsteady_1d.c b/src/test_sdis_unsteady_1d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -73,7 +73,7 @@ static const size_t nreferences = sizeof(references)/sizeof(*references); SOLID_PROP(calorific_capacity, 2000.0) /* [J/K/kg] */ SOLID_PROP(thermal_conductivity, 0.5) /* [W/m/K] */ SOLID_PROP(volumic_mass, 2500.0) /* [kg/m^3] */ -SOLID_PROP(delta, 1.0/60.0) +SOLID_PROP(delta, 1.0/80.0) #undef SOLID_PROP static double /* [K] */ diff --git a/src/test_sdis_unsteady_analytic_profile.c b/src/test_sdis_unsteady_analytic_profile.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_unsteady_analytic_profile_2d.c b/src/test_sdis_unsteady_analytic_profile_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_unsteady_atm.c b/src/test_sdis_unsteady_atm.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_utils.c b/src/test_sdis_utils.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_utils.h b/src/test_sdis_utils.h @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_volumic_power.c b/src/test_sdis_volumic_power.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 @@ -50,7 +50,7 @@ #define T0 320 #define LAMBDA 0.1 #define P0 10 -#define DELTA 1.0/55.0 +#define DELTA 1.0/60.0 /******************************************************************************* * Helper functions diff --git a/src/test_sdis_volumic_power2.c b/src/test_sdis_volumic_power2.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_volumic_power2_2d.c b/src/test_sdis_volumic_power2_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_volumic_power3_2d.c b/src/test_sdis_volumic_power3_2d.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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 diff --git a/src/test_sdis_volumic_power4.c b/src/test_sdis_volumic_power4.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2016-2024 |Méso|Star> (contact@meso-star.com) +/* Copyright (C) 2016-2025 |Méso|Star> (contact@meso-star.com) * * 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