tset default density and add temporal functionality - granular - granular dynamics simulation
 (HTM) git clone git://src.adamsgaard.dk/granular
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
 (DIR) commit b56ddb27d47ed549245aedb47bbed51949c8c36e
 (DIR) parent b5116979a34a77b651a34b8e51665749e20dd2f0
 (HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
       Date:   Thu, 22 Apr 2021 13:48:09 +0200
       
       set default density and add temporal functionality
       
       Diffstat:
         M grain.c                             |      27 +++++++++++++++++++++++++++
         M grain.h                             |       4 ++++
         M granular.c                          |       1 +
         M simulation.c                        |      50 +++++++++++++++++++++++++++++--
         M simulation.h                        |       4 +++-
       
       5 files changed, 83 insertions(+), 3 deletions(-)
       ---
 (DIR) diff --git a/grain.c b/grain.c
       t@@ -40,6 +40,7 @@ grain_defaults(struct grain *g)
                        g->gridpos[i]
                        = g->acc_lock[i] = 0;
                }
       +        g->density = 2.5e3;
                g->fixed = 0;
                g->rotating = 1;
                g->enabled = 1;
       t@@ -375,3 +376,29 @@ grain_register_contact(struct grain *g, size_t i, size_t j,
                                }
                        }
        }
       +
       +double
       +grain_stiffness_normal(const struct grain *g)
       +{
       +        return 2.0 / 3.0 * g->youngs_modulus / (1.0 - pow(g->poissons_ratio, 2.0));
       +}
       +
       +double
       +grain_stiffness_tangential(const struct grain *g)
       +{
       +        return 2.0 * g->youngs_modulus
       +               / ((1.0 + g->poissons_ratio) * (2.0 - g->poissons_ratio));
       +}
       +
       +double
       +grain_collision_time(const struct grain *g)
       +{
       +        double m;
       +
       +        m = grain_mass(g);
       +
       +        return fmin(PI * pow(2.0 * grain_stiffness_normal(g) / m
       +                                         - pow(g->damping_n, 2.0) / 4.0, -0.5),
       +                                PI * pow(2.0 * grain_stiffness_tangential(g) / m
       +                                         - pow(g->damping_t, 2.0) / 4.0, -0.5));
       +}
 (DIR) diff --git a/grain.h b/grain.h
       t@@ -60,4 +60,8 @@ void grain_temporal_integration(struct grain *g, double dt);
        void grain_register_contact(struct grain *g, size_t i, size_t j,
                                    double centerdist[3], double overlap);
        
       +double grain_stiffness_normal(const struct grain *g);
       +double grain_stiffness_tangential(const struct grain *g);
       +double grain_collision_time(const struct grain *g);
       +
        #endif
 (DIR) diff --git a/granular.c b/granular.c
       t@@ -63,6 +63,7 @@ main(int argc, char *argv[])
                        usage();
        
                sim_read_grains(&sim, stdin);
       +        sim_set_timestep(&sim);
                if (sim.t < sim.t_end)
                        sim_run_time_loop(&sim);
                sim_print_grains(stdout, &sim);
 (DIR) diff --git a/simulation.c b/simulation.c
       t@@ -6,6 +6,7 @@
        #include <stdlib.h>
        #include <err.h>
        #include <string.h>
       +#include <math.h>
        #include "grain.h"
        #include "grains.h"
        #include "simulation.h"
       t@@ -23,9 +24,11 @@ sim_new(void)
        void
        sim_defaults(struct simulation *sim)
        {
       -        int i;
       +        int i, ret;
        
       -        snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME);
       +        ret = snprintf(sim->name, sizeof(sim->name), DEFAULT_SIMULATION_NAME);
       +        if (ret < 0 || (size_t)ret >= sizeof(sim->name))
       +                errx(1, "%s: could not write simulation name", __func__);
                sim->ng = 0;
                for (i = 0; i < 3; i++) {
                        sim->constacc[i] = 0.0;
       t@@ -89,6 +92,24 @@ sim_print_grains_vtk(FILE *stream, const struct simulation *sim)
                grains_print_vtk(stream, sim->grains, sim->ng);
        }
        
       +void
       +sim_write_output_files(struct simulation *sim)
       +{
       +        int ret;
       +        char outfile[256];
       +        FILE *fp;
       +
       +        ret = snprintf(outfile, sizeof(outfile), "%s.grains.%05zu.tsv",
       +                       sim->name, sim->n_file++);
       +        if (ret < 0 || (size_t)ret >= sizeof(outfile))
       +                errx(1, "%s: could not write grains filename", __func__);
       +        if ((fp = fopen(outfile, "w")) != NULL) {
       +                sim_print_grains(fp, sim);
       +                fclose(fp);
       +        } else
       +                err(1, "%s: fopen: %s", __func__, outfile);
       +}
       +
        double
        sim_grain_pair_overlap(struct simulation *sim, size_t i, size_t j)
        {
       t@@ -123,6 +144,24 @@ sim_check_add_contact(struct simulation *sim, size_t i, size_t j)
                        grain_register_contact(&sim->grains[i], j, i, centerdist, overlap);
        }
        
       +/* Silbert et al. 2001 */
       +void
       +sim_set_timestep(struct simulation *sim)
       +{
       +        int ret = 0;
       +        size_t i;
       +        double t_col, t_col_min = INFINITY;
       +
       +        for (i = 0; i < sim->ng; i++) {
       +                t_col = grain_collision_time(&sim->grains[i]);
       +                if (t_col_min > t_col)
       +                        t_col_min = t_col;
       +        }
       +
       +        check_float("t_col_min", t_col_min, &ret);
       +        sim->dt = t_col_min / 50.0;
       +}
       +
        /* TODO: add grid sorting */
        void
        sim_detect_contacts(struct simulation *sim)
       t@@ -162,6 +201,13 @@ sim_step_time(struct simulation *sim)
                for (i = 0; i < sim->ng; i++)
                        grain_temporal_integration(&sim->grains[i], sim->dt);
        
       +        if ((sim->t_file >= sim->dt_file || sim->iter == 0) &&
       +                strncmp(sim->name, DEFAULT_SIMULATION_NAME,
       +                                sizeof(sim->name)) != 0) {
       +                sim_write_output_files(sim);
       +                sim->t_file = 0.0;
       +        }
       +
                sim->t += sim->dt;
                sim->t_file += sim->dt;
                sim->iter++;
 (DIR) diff --git a/simulation.h b/simulation.h
       t@@ -1,7 +1,7 @@
        #ifndef GRANULAR_SIMULATION_
        #define GRANULAR_SIMULATION_
        
       -#define DEFAULT_SIMULATION_NAME "unnamed_simulation"
       +#define DEFAULT_SIMULATION_NAME "unnamed"
        
        #include <stdio.h>
        #include "grain.h"
       t@@ -43,7 +43,9 @@ void sim_read_grains(struct simulation *sim, FILE *stream);
        
        void sim_print_grains(FILE *stream, const struct simulation *sim);
        void sim_print_grains_vtk(FILE *stream, const struct simulation *sim);
       +void sim_write_output_files(struct simulation *sim);
        
       +void sim_set_timestep(struct simulation *sim);
        void sim_detect_contacts(struct simulation *sim);
        void sim_step_time(struct simulation *sim);
        void sim_run_time_loop(struct simulation *sim);