commit ae1c884e299b19c401830dc8c5c7512e5c08faab
parent 935b2e1aa58040448edd97b0bf0f2123e989561f
Author: Vincent Forest <vincent.forest@meso-star.com>
Date: Thu, 14 Mar 2024 16:39:23 +0100
WoS algorithm now supports unsteady systems
We rely on the Star-WoS-Function library to sample the first-passage
time of a sphere or circle, depending on whether the configuration is 2D
or 3D, respectively.
Diffstat:
5 files changed, 132 insertions(+), 3 deletions(-)
diff --git a/Makefile b/Makefile
@@ -90,6 +90,8 @@ libsdis.o: $(OBJ)
echo "senc3d $(SENC3D_VERSION) not found" >&2; exit 1; fi
@if ! $(PKG_CONFIG) --atleast-version $(SSP_VERSION) star-sp; then \
echo "star-sp $(SSP_VERSION) not found" >&2; exit 1; fi
+ @if ! $(PKG_CONFIG) --atleast-version $(SWF_VERSION) swf; then \
+ echo "swf $(SWF_VERSION) not found" >&2; exit 1; fi
@echo "config done" > $@
.SUFFIXES: .c .d .o
diff --git a/config.mk b/config.mk
@@ -63,6 +63,10 @@ SSP_VERSION = 0.14
SSP_CFLAGS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags star-sp)
SSP_LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs star-sp)
+SWF_VERSION = 0.0
+SWF_CFLAGS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags swf)
+SWF_LIBS = $$($(PKG_CONFIG) $(PCFLAGS) --libs swf)
+
# For tests only
S3DUT_VERSION = 0.4
S3DUT_CFLAGS = $$($(PKG_CONFIG) $(PCFLAGS) --cflags s3dut)
@@ -75,6 +79,7 @@ DPDC_CFLAGS =\
$(SENC2D_CFLAGS)\
$(SENC3D_CFLAGS)\
$(SSP_CFLAGS)\
+ $(SWF_CFLAGS)\
$($(DISTRIB_PARALLELISM)_CFLAGS)\
-fopenmp
DPDC_LIBS =\
@@ -84,6 +89,7 @@ DPDC_LIBS =\
$(SENC2D_LIBS)\
$(SENC3D_LIBS)\
$(SSP_LIBS)\
+ $(SWF_LIBS)\
$($(DISTRIB_PARALLELISM)_LIBS)\
-lm\
-fopenmp
diff --git a/src/sdis_device.c b/src/sdis_device.c
@@ -25,6 +25,7 @@
#include <star/s2d.h>
#include <star/s3d.h>
#include <star/ssp.h>
+#include <star/swf.h>
#include <omp.h>
@@ -251,6 +252,37 @@ error:
}
static INLINE res_T
+setup_starwf(struct sdis_device* dev)
+{
+ struct swf_H_tabulate_args H2d_args = SWF_H2D_TABULATE_ARGS_DEFAULT;
+ struct swf_H_tabulate_args H3d_args = SWF_H3D_TABULATE_ARGS_DEFAULT;
+ res_T res = RES_OK;
+ ASSERT(dev);
+
+ H2d_args.allocator = dev->allocator;
+ H3d_args.allocator = dev->allocator;
+
+ res = swf_H2d_tabulate(&H2d_args, &dev->H_2d);
+ if(res != RES_OK) {
+ log_err(dev, "Unable to tabulate H2d function -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+ res = swf_H3d_tabulate(&H3d_args, &dev->H_3d);
+ if(res != RES_OK) {
+ log_err(dev, "Unable to tabulate H3d function -- %s.\n",
+ res_to_cstr(res));
+ goto error;
+ }
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
+static INLINE res_T
setup_mpi(struct sdis_device* dev, const struct sdis_device_create_args* args)
{
ASSERT(dev && args);
@@ -279,6 +311,8 @@ device_release(ref_T* ref)
dev = CONTAINER_OF(ref, struct sdis_device, ref);
if(dev->s2d_dev) S2D(device_ref_put(dev->s2d_dev));
if(dev->s3d_dev) S3D(device_ref_put(dev->s3d_dev));
+ if(dev->H_2d) SWF(tabulation_ref_put(dev->H_2d));
+ if(dev->H_3d) SWF(tabulation_ref_put(dev->H_3d));
if(dev->logger == &dev->logger__) logger_release(&dev->logger__);
ASSERT(flist_name_is_empty(&dev->interfaces_names));
ASSERT(flist_name_is_empty(&dev->media_names));
@@ -342,6 +376,8 @@ sdis_device_create
if(res != RES_OK) goto error;
res = setup_star3d(dev);
if(res != RES_OK) goto error;
+ res = setup_starwf(dev);
+ if(res != RES_OK) goto error;
res = setup_mpi(dev, args);
if(res != RES_OK) goto error;
diff --git a/src/sdis_device_c.h b/src/sdis_device_c.h
@@ -36,6 +36,7 @@
struct mutex;
struct ssp_rng;
struct ssp_rng_proxy;
+struct swf_tabulation;
struct name { FITEM; void* mem; };
#define FITEM_TYPE name
@@ -64,6 +65,9 @@ struct sdis_device {
struct s2d_device* s2d_dev;
struct s3d_device* s3d_dev;
+ struct swf_tabulation* H_2d;
+ struct swf_tabulation* H_3d;
+
ref_T ref;
};
diff --git a/src/sdis_heat_path_conductive_wos_Xd.h b/src/sdis_heat_path_conductive_wos_Xd.h
@@ -13,10 +13,13 @@
* 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_device_c.h"
#include "sdis_heat_path_conductive_c.h"
#include "sdis_medium_c.h"
#include "sdis_scene_c.h"
+#include <star/swf.h>
+
#include "sdis_Xd_begin.h"
/*******************************************************************************
@@ -49,6 +52,69 @@ error:
goto exit;
}
+static res_T
+XD(time_travel)
+ (struct sdis_scene* scn,
+ struct XD(rwalk)* rwalk,
+ struct ssp_rng* rng,
+ const double alpha, /* Diffusivity, i.e. lambda/(rho*cp) */
+ const double t0, /* Initial time [s] */
+ const double distance, /* Displacement. Must be multipled by fp_to_meter */
+ struct XD(temperature)* T)
+{
+ double dst = 0; /* Distance [m] */
+ double tau = 0; /* Time [s] */
+ double x = 0;
+ double r = 0;
+ double temperature = 0; /* [k] */
+ res_T res = RES_OK;
+ ASSERT(scn && rwalk && rng && alpha && T);
+
+ /* No displacement => no time travel */
+ if(distance == 0) goto exit;
+
+ /* Sample x = tau*alpha/distance^2 */
+ r = ssp_rng_canonical(rng);
+ x = swf_tabulation_inverse(XD(scn->dev->H), SWF_QUADRATIC, r);
+
+ /* Retrieve the time to travel */
+ dst = distance * scn->fp_to_meter;
+ tau = x / alpha * dst * dst;
+
+ /* Increment the elapsed time */
+ rwalk->elapsed_time += MMIN(tau, rwalk->vtx.time - t0);
+
+ if(IS_INF(rwalk->vtx.time)) goto exit; /* Steady computation */
+
+ /* Let's take a trip back in time */
+ rwalk->vtx.time -= MMAX(t0, rwalk->vtx.time - tau);
+
+ /* Thepath does not reach the initial condition */
+ if(rwalk->vtx.time > t0) goto exit;
+
+ /* Fethc the initial temperature */
+ temperature = medium_get_temperature(rwalk->mdm, &rwalk->vtx);
+ if(temperature < 0) {
+ log_err(scn->dev,
+ "%s:%s: the path reaches the initial condition but the "
+ "%s temperature remains unknown -- position=%g, %g, %g\n",
+ __FILE__, FUNC_NAME,
+ medium_type_to_string(sdis_medium_get_type(rwalk->mdm)),
+ SPLIT3(rwalk->vtx.P));
+ res = RES_BAD_ARG;
+ goto error;
+ }
+
+ /* Update the temperature */
+ T->value += temperature;
+ T->done = 1;
+
+exit:
+ return res;
+error:
+ goto exit;
+}
+
#if DIM == 2
static INLINE enum sdis_side
compute_hit_side_2d
@@ -276,7 +342,8 @@ static res_T
XD(sample_next_position)
(struct sdis_scene* scn,
struct XD(rwalk)* rwalk,
- struct ssp_rng* rng)
+ struct ssp_rng* rng,
+ double* distance) /* Displacement distance */
{
/* Intersection */
struct sXd(hit) hit = SXD_HIT_NULL;
@@ -289,7 +356,7 @@ XD(sample_next_position)
/* Miscellaneous */
res_T res = RES_OK;
- ASSERT(rwalk && rng);
+ ASSERT(rwalk && rng && distance);
/* Find the closest distance from the current position to the geometry */
wos_radius = (float)INF;
@@ -374,6 +441,7 @@ XD(sample_next_position)
}
exit:
+ *distance = hit.distance;
return res;
error:
goto exit;
@@ -393,6 +461,7 @@ XD(conductive_path_wos)
/* Properties */
struct solid_props props_ref = SOLID_PROPS_NULL;
struct solid_props props = SOLID_PROPS_NULL;
+ double alpha = 0; /* diffusivity, i.e. lambda/(rho*cp) */
/* Miscellaneous */
size_t ndiffusion_steps = 0; /* For debug */
@@ -418,13 +487,25 @@ XD(conductive_path_wos)
* with respect to the properties of the reinjected position. */
solid_get_properties(rwalk->mdm, &rwalk->vtx, &props_ref);
+ /* The algorithm assumes that lambda, rho and cp are constants. The
+ * diffusivity of the material (alpha) can therefore be calculated once */
+ alpha = props_ref.lambda / (props_ref.rho * props_ref.cp);
+
/* Sample a diffusive path */
for(;;) {
+ double distance = 0;
/* Find the next position of the conductive path */
- res = XD(sample_next_position)(scn, rwalk, rng);
+ res = XD(sample_next_position)(scn, rwalk, rng, &distance);
+ if(res != RES_OK) goto error;
+
+ /* Going back in time */
+ res = XD(time_travel)(scn, rwalk, rng, alpha, props_ref.t0, distance, T);
if(res != RES_OK) goto error;
+ /* The path reaches the initial condition */
+ if(T->done) break;
+
/* The path reaches a boundary */
if(!SXD_HIT_NONE(&rwalk->hit)) break;
