dem_granular 0.1.3

Granular physics for MDDEM: Hertz normal contact, Mindlin tangential friction, rotational dynamics
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

use std::{
    fs::{self, OpenOptions},
    io::Write,
};

use mddem_app::prelude::*;
use mddem_scheduler::prelude::*;

use mddem_core::{Atom, CommResource, Input, RunConfig, RunState};

/// Writes granular temperature (average KE per particle) to a data file each thermo interval.
pub struct GranularTempPlugin;

impl Plugin for GranularTempPlugin {
    fn build(&self, app: &mut App) {
        app.add_update_system(print_granular_temperature, ScheduleSet::PreExchange);
    }
}

pub fn print_granular_temperature(
    atoms: Res<Atom>,
    run_state: Res<RunState>,
    comm: Res<CommResource>,
    run_config: Res<RunConfig>,
    scheduler_manager: Res<SchedulerManager>,
    input: Res<Input>,
) {
    let index = scheduler_manager.index;
    let thermo_interval = run_config.current_stage(index).thermo;
    if !run_state.total_cycle.is_multiple_of(thermo_interval) {
        return;
    }
    let nlocal = atoms.nlocal as usize;
    let mut local_mv_x = 0.0;
    let mut local_mv_y = 0.0;
    let mut local_mv_z = 0.0;
    let mut local_mass = 0.0;
    for i in 0..nlocal {
        local_mv_x += atoms.mass[i] * atoms.vel[i][0];
        local_mv_y += atoms.mass[i] * atoms.vel[i][1];
        local_mv_z += atoms.mass[i] * atoms.vel[i][2];
        local_mass += atoms.mass[i];
    }
    let global_mv_x = comm.all_reduce_sum_f64(local_mv_x);
    let global_mv_y = comm.all_reduce_sum_f64(local_mv_y);
    let global_mv_z = comm.all_reduce_sum_f64(local_mv_z);
    let global_mass = comm.all_reduce_sum_f64(local_mass);
    let avg_vx = global_mv_x / global_mass;
    let avg_vy = global_mv_y / global_mass;
    let avg_vz = global_mv_z / global_mass;
    let mut vel_diff = 0.0;
    for i in 0..nlocal {
        vel_diff += atoms.mass[i] * (atoms.vel[i][0] - avg_vx).powi(2)
            + atoms.mass[i] * (atoms.vel[i][1] - avg_vy).powi(2)
            + atoms.mass[i] * (atoms.vel[i][2] - avg_vz).powi(2);
    }
    let vel_diff_sum = comm.all_reduce_sum_f64(vel_diff);
    let granular_temperature = vel_diff_sum / (3.0 * global_mass);

    // KE computation — all ranks must participate in allreduce
    let mut local_ke = 0.0;
    for i in 0..nlocal {
        local_ke += 0.5
            * atoms.mass[i]
            * (atoms.vel[i][0].powi(2) + atoms.vel[i][1].powi(2) + atoms.vel[i][2].powi(2));
    }
    let global_ke = comm.all_reduce_sum_f64(local_ke);

    // Total momentum magnitude (should be ~conserved)
    let mom_mag = (global_mv_x * global_mv_x + global_mv_y * global_mv_y + global_mv_z * global_mv_z).sqrt();

    if comm.rank() != 0 {
        return;
    }
    let physical_time = run_state.total_cycle as f64 * atoms.dt;
    let base_dir = match input.output_dir.as_deref() {
        Some(dir) => format!("{}/data", dir),
        None => "data".to_string(),
    };
    let result = fs::create_dir_all(&base_dir);
    if let Err(_error) = result {
        println!("Could not create file directory {}", base_dir)
    }
    let data_path = format!("{}/GranularTemp.txt", base_dir);
    let mut file = if run_state.total_cycle == 0 {
        OpenOptions::new()
            .create(true)
            .write(true)
            .truncate(true)
            .open(&data_path)
            .unwrap()
    } else {
        OpenOptions::new()
            .create(true)
            .append(true)
            .open(&data_path)
            .unwrap()
    };

    writeln!(
        &mut file,
        "{} {:.6e} {:.10e} {:.10e} {:.10e}",
        run_state.total_cycle, physical_time, granular_temperature, global_ke, mom_mag
    )
    .unwrap();
}