rust snippets/crates

2 dimensional array

https://stackoverflow.com/questions/13212212/creating-two-dimensional-arrays-in-rust

let mut state = [[0u8; 4]; 6];
state[0][1] = 42;

unknown toml

https://stackoverflow.com/questions/72690767/how-to-parse-toml-in-rust-with-unknown-structure

float comparison

https://stackoverflow.com/questions/41447678/comparison-of-two-floats-in-rust-to-arbitrary-level-of-precision

criterion: benchmark/profiling lines of code

criterion https://docs.rs/criterion/latest/criterion/

parsing cli arguments

https://rust-cli.github.io/book/tutorial/cli-args.html

creating dir

use std::fs;
fs::create_dir("MyDir").expect("problem creating the directory");

csv

https://docs.rs/csv/latest/csv/

c++/rust interop

https://cxx.rs/

bindgen c / rust interop

https://rust-lang.github.io/rust-bindgen/

code dependencies

flamegraph-rs

sudo apt install linux-tools-common linux-tools-generic linux-tools-`uname -r`

rgsl

for rgsl: (debian)

sudo apt-get install libgsl0-dev

orga code building

cheatsheet

ndarray rust

readme quickstart
ndarray_for_numpy_users

todos

improve performance

rayon integration in ndarray

Zip also implements NdarrayIntoParallelIterator, and there is an extension trait so that it can use a method .par_apply directly.

https://docs.rs/ndarray/latest/ndarray/parallel/index.html

benchmarking

example in stack https://stackoverflow.com/questions/13322479/benchmarking-programs-in-rust

links given benchmark test

using ArrayView for things needing only views

using Zip! instead of indexing

https://docs.rs/ndarray/latest/ndarray/struct.Zip.html

WasserMarder advice

“Getters and setters are typically not idomatic rust but are expected to be optimized out by the compiler. Performance sensitive code should try to avoid indexing when possible but instead use iterators. The Zip/azip! functionality from ndarray should be what you are looking for.”

verify with samy

eta0

shear visc = eta0 * m * rho

eta0 is the dynamic viscosity? (changed the definition in fcts)

zeta0

eta0 * mass * rho

cin_visc = eta0?

understand why there is a jump in some simulations

vdw file with jump?

use std::io::Write; // to use “write_all” method use std::fs; // to read/write a file contents use rgsl::logarithm::log; use rgsl::exponential::exp; use crate::maths::fcts::*; use crate::maths::mystructs::*;

fn create_scalar_grid(ncol_size: i32, nrow_size: i32) -> ScalarField2D { return ScalarField2D { s: vec![vec![0.; ncol_size as usize]; nrow_size as usize], }; }

fn create_vector_grid(ncol_size: i32, nrow_size: i32) -> VectorField2D { let result = VectorField2D { x: vec![vec![0.; ncol_size as usize]; nrow_size as usize], y: vec![vec![0.; ncol_size as usize]; nrow_size as usize] }; return result; }

fn create_tensor_grid(ncol_size: i32, nrow_size: i32) -> TensorField2D { let result = TensorField2D { xx: vec![vec![0.; ncol_size as usize]; nrow_size as usize], xy: vec![vec![0.; ncol_size as usize]; nrow_size as usize], yx: vec![vec![0.; ncol_size as usize]; nrow_size as usize], yy: vec![vec![0.; ncol_size as usize]; nrow_size as usize] }; return result; }

pub fn do_sim() { let do_vdw_sim = true; // let do_vdw_sim = false; // env::set_var(“RUST_BACKTRACE”, “1”); let output_dir = “./testoutput”; let output_dir = “/home/mehdi/workdir/dossiers/ilm/these/code_simulations/rust_implementation/pfbox_git/src/testoutput”; let path = format!(“{}/log.txt”, output_dir); // Open/Create a file in write-only mode, returns `io::Result<File>` let file = match fs::File::create(&path) { Err(why) => panic!(“couldn’t create logfile {}: {}”, path, why), Ok(file) => file, }; // Open an already created file let mut file = fs::OpenOptions::new() .write(true) .append(true) // This is needed to append to file .open(&path) .unwrap();

let str_to_append = format!(“new sim\n”); // appending the string to file.write_all(&str_to_append.as_bytes());

// System initialisation let mut step = 0; let max_time_step = 1_000; // let step_count_before_save = max_time_step/20; let step_count_before_save = max_time_step/10;

let print_frequency = 20.; let mut print_percertage_threshold = 100./print_frequency;

let dt = 1e-2; let mut time = 0.;

let ncol_size = 100; let nrow_size = 2;

let rho_liq0 = 0.8; let ln_rho_liq0 = log(rho_liq0); let rho_vap0 = 0.02; let ln_rho_vap0 = log(rho_vap0); let temp0 = 0.7;

let box_info = BoxInfo{col_max: ncol_size, row_max: nrow_size};

if do_vdw_sim {

//auie physics quantities definition //////////////////////////////////////////////////////////////////////////// // Physics quantities definition //////////////////////////////////////////////////////////////////////////// // GD for grid let mut GD_rho = create_scalar_grid(ncol_size, nrow_size); let mut GD_temp = create_scalar_grid(ncol_size, nrow_size); let mut GD_pressure = create_tensor_grid(ncol_size, nrow_size); // momentum, also known as J = rho*velocity let mut GD_J = create_vector_grid(ncol_size, nrow_size); // velocity let mut GD_v = create_vector_grid(ncol_size, nrow_size);

//auie quantities used for the computations definition //////////////////////////////////////////////////////////////////////////// // Quantities used for the computations definition ////////////////////////////////////////////////////////////////////////////

let mut GD_ln_rho = create_scalar_grid(ncol_size, nrow_size); let mut GD_grad_rho = create_vector_grid(ncol_size, nrow_size); let mut GD_lap_rho = create_scalar_grid(ncol_size, nrow_size);

let mut GD_vJ = create_tensor_grid(ncol_size, nrow_size);

let mut GD_grad_v = create_tensor_grid(ncol_size, nrow_size);

let mut GD_div_v = create_scalar_grid(ncol_size, nrow_size);

let mut GD_traceless_grad_v = create_tensor_grid(ncol_size, nrow_size);

let mut GD_lap_v = create_vector_grid(ncol_size, nrow_size);

let mut GD_div_vJ = create_vector_grid(ncol_size, nrow_size); let mut GD_grad_div_v = create_vector_grid(ncol_size, nrow_size);

let mut GD_div_press = create_vector_grid(ncol_size, nrow_size);

let mut GD_ln_rho_traceless_grad_v = create_vector_grid(ncol_size, nrow_size);

let inv_cv = 1.0/(1.5*kB);

let mut GD_traceless_grad_v_dyadic_grad_v = create_scalar_grid(ncol_size, nrow_size);

let mut GD_grad_ln_rho_scalar_grad_T = create_scalar_grid(ncol_size, nrow_size);

let mut GD_grad_ln_rho = create_vector_grid(ncol_size, nrow_size);

let mut GD_v_scalar_grad_ln_rho = create_scalar_grid(ncol_size, nrow_size);

let mut GD_grad_ln_rho_traceless_grad_v = create_vector_grid(ncol_size, nrow_size);

let mut GD_grad_T = create_vector_grid(ncol_size, nrow_size);

let mut GD_lap_T = create_scalar_grid(ncol_size, nrow_size);

let mut GD_v_scalar_grad_T = create_scalar_grid(ncol_size, nrow_size);

//auie fluid initial state //////////////////////////////////////////////////////////////////////////// // Fluid initial state ////////////////////////////////////////////////////////////////////////////

for col in 0usize..ncol_size as usize { for row in 0usize..nrow_size as usize { // putting liquid in the first half if ((col as i32) < ncol_size/2){ GD_rho.set_pos(row, col, &rho_liq0); GD_ln_rho.set_pos(row, col, &ln_rho_liq0);} else {GD_rho.set_pos(row, col, &rho_vap0); GD_ln_rho.set_pos(row, col, &ln_rho_vap0);}

// setting initial temperature GD_temp.set_pos(row, col, &temp0); }}

//auie computation variables update //////////////////////////////////////////////////////////////////////////// // Computations variables update ////////////////////////////////////////////////////////////////////////////

//auie time loop for i_time_step in 0..max_time_step {

step = i_time_step; let percentage_done = 100.*(step as f64/max_time_step as f64); if (percentage_done > print_percertage_threshold) { print_percertage_threshold += 100./print_frequency; println!(“completed {percentage_done:.1}%”); }

// update of computations variables for col in 0usize..ncol_size as usize { for row in 0usize..nrow_size as usize {

let col_i32 = col as i32; let row_i32 = row as i32;

// ——————————————————- // GD_lap_rho begin update GD_lap_rho .set_pos(row, col, &laplacian(&GD_rho, row_i32, col_i32, &box_info)); // GD_lap_rho end update // ——————————————————-

// ——————————————————- // GD_lap_T begin update GD_lap_T .set_pos(row, col, &laplacian(&GD_temp, row_i32, col_i32, &box_info)); // GD_lap_T end update // ——————————————————-

// ——————————————————- // GD_grad_T begin update GD_grad_T .set_pos(row, col, &gradient(&GD_temp, row_i32, col_i32, &box_info)); // GD_grad_T end update // ——————————————————-

// ——————————————————- // GD_grad_rho begin update GD_grad_rho .set_pos(row, col, &gradient(&GD_rho, row_i32, col_i32, &box_info)); // GD_grad_rho end update // ——————————————————-

// ——————————————————- // GD_ln_rho begin update // :todo:log: let rho = GD_rho.get_pos(row, col); if (rho < 0.) { let str_to_append = format!(“step {}, col={}, row={}\n\ neg log: {}\n\ ————\n”, &step, &col, &row, &rho); // appending the string to file.write_all(&str_to_append.as_bytes()) .expect(“write failed”); println!(“error step {}:\n\ negative rho: rho = {}”, step, rho); GD_ln_rho .set_pos(row, col, &0.);} else { let ln_rho = log(rho); GD_ln_rho .set_pos(row, col, &ln_rho);} // GD_ln_rho end update // ——————————————————-

// ——————————————————- // GD_lap_v begin update GD_lap_v .set_pos(row, col, &laplacian_vector(&GD_v, row_i32, col_i32, &box_info)); // GD_lap_v end update // ——————————————————-

// ——————————————————- // GD_div_v begin update GD_div_v .set_pos(row, col, &div_vector(&GD_v, row_i32, col_i32, &box_info)); // GD_div_v end update // ——————————————————-

// ——————————————————- // GD_grad_v begin update GD_grad_v .set_pos(row, col, &gradient_vector(&GD_v, row_i32, col_i32, &box_info)); // GD_grad_v end update // ——————————————————-

// ——————————————————- // GD_grad_ln_rho begin update GD_grad_ln_rho .set_pos(row, col, &gradient(&GD_ln_rho, row_i32, col_i32, &box_info)); // GD_grad_ln_rho end update // ——————————————————-

// ——————————————————- // GD_grad_div_v begin update GD_grad_div_v .set_pos(row, col, &grad_div_vel(&GD_v, row_i32, col_i32, &box_info)); // GD_grad_div_v end update // ——————————————————-

// ——————————————————- // GD_traceless_grad_v begin update { let grad_v = GD_grad_v.get_pos(row, col); let div_v = GD_div_v.get_pos(row, col);

let traceless_grad_v = tens2D { xx: 2.*grad_v.xx - (2./(1.*dim as f64)) * div_v, xy: grad_v.xy + grad_v.yx, yx: grad_v.xy + grad_v.yx, yy: 2.*grad_v.yy - (2./(1.*dim as f64)) * div_v};

GD_traceless_grad_v.set_pos(row, col, &traceless_grad_v); } // GD_traceless_grad_v end update // ——————————————————-

// ——————————————————- // GD_vJ begin update { let v = GD_v.get_pos(row, col); let J = GD_J.get_pos(row, col); let tens_vJ = tens2D{ xx: v.x * J.x, xy: v.x * J.y, yx: v.y * J.x, yy: v.y * J.y }; GD_vJ .set_pos(row, col, &tens_vJ) } // GD_vJ end update // ——————————————————-

// ——————————————————- // GD_div_vJ begin update GD_div_vJ .set_pos(row, col, &div_tensor(&GD_vJ, row_i32, col_i32, &box_info)); // GD_div_vJ end update // ——————————————————-

// ——————————————————- // GD_v_scal_grad_T begin update GD_v_scalar_grad_T .set_pos(row, col, &scal_product(&GD_v.get_pos(row, col), &GD_grad_T.get_pos(row, col))); // GD_v_scal_grad_T end update // ——————————————————-

// ——————————————————- // GD_traceless_grad_v_dyadic_grad_v begin update GD_traceless_grad_v_dyadic_grad_v .set_pos(row, col, &dyadic_product(&GD_traceless_grad_v.get_pos(row, col), &GD_grad_v.get_pos(row, col))); // GD_traceless_grad_v_dyadic_grad_v end update // ——————————————————-

// ——————————————————- // GD_v_scalar_grad_ln_rho begin update GD_v_scalar_grad_ln_rho .set_pos(row, col, &scal_product(&GD_v.get_pos(row, col), &GD_grad_ln_rho.get_pos(row, col))); // GD_v_scalar_grad_ln_rho end update // ——————————————————-

// ——————————————————- // GD_pressure begin update GD_pressure .set_pos(row, col, &pressure(GD_rho.get_pos(row, col), &GD_grad_rho.get_pos(row, col), GD_lap_rho.get_pos(row, col), GD_temp.get_pos(row, col))); // GD_pressure end update // ——————————————————-

// ——————————————————- // GD_grad_ln_rho_traceless_grad_v begin update

GD_grad_ln_rho_traceless_grad_v .set_pos(row, col, &tens_product_vec( &GD_traceless_grad_v.get_pos(row, col), &GD_grad_ln_rho.get_pos(row, col))); // GD_grad_ln_rho_traceless_grad_v end update // ——————————————————-

// ——————————————————- // GD_grad_ln_rho_scalar_grad_T begin update GD_grad_ln_rho_scalar_grad_T .set_pos(row, col, &scal_product(&GD_grad_ln_rho.get_pos(row, col), &GD_grad_T.get_pos(row, col))); // GD_grad_ln_rho_scalar_grad_T end update // ——————————————————-

// ——————————————————- // div_press begin update GD_div_press .set_pos(row, col, &div_tensor(&GD_pressure, row_i32, col_i32, &box_info)); // div_press end update // ——————————————————-

}} // updating computations values end parenthesis

//bépo WRITING part

if (step % step_count_before_save == 0) {

let filename = format!(“{}/step”, output_dir, i_time_step); let mut file = fs::File::create(&filename) .expect(“couldn’t create log file”);

file.write_all( “# column density temperature\n”.as_bytes()) .expect(“write failed”);

let rho_profile = GD_rho.x_profile(); let temp_profile = GD_temp.x_profile();

for col_index in 0usize..ncol_size as usize { let str_to_append = format!(“{} {} {}\n”, &col_index, &rho_profile[col_index], &temp_profile[col_index]);

file.write_all(&str_to_append.as_bytes()) .expect(“write failed”); }} // let str_to_append = format!(“step {}, i={}, j={}\n\ // neg log {}\n\ // ————\n”, // &step, &i, &j, &rho); // // appending the string to // file.write_all(&str_to_append.as_bytes()) // .expect(“write failed”);

//auie main loop //////////////////////////////////////////////////////////////////////////// // Main loop ////////////////////////////////////////////////////////////////////////////

for row in 0usize..nrow_size as usize { for col in 0usize..ncol_size as usize {

let row_i32 = row as i32; let col_i32 = col as i32;

let div_vJ = GD_div_vJ.get_pos(row, col); let rho = GD_rho.get_pos(row, col); let lap_v = GD_lap_v.get_pos(row, col); let grad_div_v = GD_grad_div_v.get_pos(row, col); let grad_ln_rho_traceless_grad_v = GD_grad_ln_rho_traceless_grad_v.get_pos(row, col); let grad_ln_rho = GD_grad_ln_rho.get_pos(row, col); let div_v = GD_div_v.get_pos(row, col); let div_press = GD_div_press.get_pos(row, col); let ln_rho = GD_ln_rho.get_pos(row, col); let v_grad_ln_rho = GD_v_scalar_grad_ln_rho.get_pos(row, col); let temp = GD_temp.get_pos(row, col); let traceless_grad_v_dyadic_grad_v = GD_traceless_grad_v_dyadic_grad_v.get_pos(row, col); let grad_ln_rho_scalar_grad_T = GD_grad_ln_rho_scalar_grad_T.get_pos(row, col); let lap_T = GD_lap_T.get_pos(row, col); let v_scalar_grad_T = GD_v_scalar_grad_T.get_pos(row, col); let J = GD_J.get_pos(row, col);

//bépo MOMENTUM conservation

let mut new_J = vec2D { x: J.x + (- div_vJ.x

eta0 * rho * lap_v.x
- eta0 * (1.-2./(1.*dim as f64) + zeta0)
- rho * grad_div_v.x
eta0 * rho * grad_ln_rho_traceless_grad_v.x
- zeta0 * rho * grad_ln_rho.x * div_v
- div_press.x)
- dt,
y: J.y + (- div_vJ.y
eta0 * rho * lap_v.y
- eta0 * (1.-2./(1.*dim as f64) + zeta0)
- rho * grad_div_v.y
eta0 * rho * grad_ln_rho_traceless_grad_v.y
- zeta0 * rho * grad_ln_rho.y * div_v
- div_press.y)
- dt

};

// if you want gravity // J.y += -rho * gravity * dt;

GD_J.set_pos(row, col, &new_J);

//bépo MASS conservation

// without ln_rho : //rho[i][j] -= div_J[i][j]*dt;

let mut new_ln_rho = ln_rho - (div_v + v_grad_ln_rho) * dt; let mut new_rho = exp(new_ln_rho);

GD_ln_rho.set_pos(row, col, &new_ln_rho); GD_rho.set_pos(row, col, &new_rho);

//bépo VELOCITY from momentum GD_v.set_pos(row, col, &vec2D{x: new_J.x/new_rho, y: new_J.y/new_rho});

//bépo TEMPERATURE ENERGY conservation

// term l div_v

let mut new_T = temp + inv_cv * ( // term l div_v -kB * temp * (1. + rho * b/(1.-rho * b)) * div_v // term dissipative_stress_grad_v

eta0 * traceless_grad_v_dyadic_grad_v
zeta0 * div_v * div_v

// term laplacian T

lambda0 * (grad_ln_rho_scalar_grad_T + lap_T)

) * dt

v_scalar_grad_T * dt;

GD_temp.set_pos(row, col, &new_T);

}} // i, j loop closing parenthesis } // time step closing parenthesis } // if vdw_simu closing parenthesis } // main definition closing parenthesis

create a module for the simulation, outside of main

verify derivatives direction

et pas oublier row_dx/col_dx dans la def de BoxInfo

j’ai mis comme c’est fait dans le code de Samy, mais ça me semblait bizarre la direction dans laquelle c’était fait…
c’est dans les fcts:
partial_deriv —————————— !
grad_scalar —————————— !
gradient —————————— !
gradient_vector —————————— !
div_vector —————————— !
div_tensor —————————— !

fix the grad_scalar test unit

deal with negative log values (:todo:log:)

implementation changes

functions

for the gradient func, created partial derivative function

functions args

removed the last argument, and returns the changed argument

v_nabla_v

gets grad_v as argument, maybe not necessary!

function names

tens_vec_product

rust

tens_product_vec

c

tens_vec_product

renaming

traceless_grad_v_grad_v renamed traceless_grad_v_dyadic_grad_v
grad_ln_rho_grad_T renamed grad_ln_rho_scalar_grad_T
v_grad_T renamed v_scal_grad_T
grad_ln_rho_grad_T renamed grad_ln_rho_scalar_grad_T
v_grad_T renamed v_scal_grad_T

done

plots

adapt simplot to plot the results of my simulation on rust

State “DONE” from [2022-12-01 Thu 19:45]

physics loop construction

fix the derive problem

State “DONE” from “TODO” [2022-12-12 lun. 18:22]

before changing rows/cols

hints

with big times, we can see that the density derives…

v_x and v_y are inverted in my code and samy’s code!
seems like v_y takes a long time to tend to 0… maybe the viscosity? or something in the J conservation?
idées:
- différence d’encodage? non, les deux sont en IEEE 754 double precision
- mes x/y sont pas en accord (j’ai l’impression d’avoir bien vérifié ça tho…)

problem in pressure (coming from gradient scalar)

PROBLEM: dans pressure grad_rho n’est pas le même

pressure is updated one step just before step 0 in C also the values are good but:

grad_rho_IN_C_STEP_-1.x = grad_rho_IN_RUST_STEP_0.y

this comes from the scalar gradient computation!

problem in scalar gradient

prints

C-1: pressure debug init PPPPPPPPPPPPPPPPP P_th: -0.047407 P_iso: 2.982446 grad_rho00.x: 0.772200 grad_rho00.y: 0.000000 pressure debug end PPPPPPPPPPPPPPPPP

Cstep0: div_vJ.x: 0.000000 div_vJ.x: 0.000000 lap_v.x: 0.000000 lap_v.x: 0.000000 div_v: 0.000000 vx: 0.003279 vy: 0.000000 Jx: 0.002623 Jy: 0.000000 Jx: 0.002623 Jy: 0.000000 div_press_x: -0.262327 div_press_y: 0.000000 press_xx: 3.578739 press_xy: 0.000000 press_yx: 0.000000 press_yy: 2.982446 percentage done = 0.000000

pressure debug init PPPPPPPPPPPPPPPPP P_th: -0.047407 P_iso: 2.982446 grad_rho00.x: 0.772200 grad_rho00.y: 0.000000 pressure debug end PPPPPPPPPPPPPPPPP

Ruststep0: —————–[I] STEP 0 GD_div_vJ.x = 0.00000000 GD_div_vJ.y = 0.00000000 GD_lap_v.x = 0.00000000 GD_lap_v.y = 0.00000000 GD_div_v = 0.00000000 GD_v.x = 0.00000000 GD_v.y = 0.00000000 GD_J.x = 0.00000000 GD_J.y = 0.00000000 GD_div_press.x = 0.00000000 GD_div_press.y = 0.00000000 GD_pressure.xx = 0.00000000 GD_pressure.xy = 0.00000000 GD_pressure.yx = 0.00000000 GD_pressure.yy = 0.00000000 p_thermo = -0.04740741 p_iso = 2.98244617 grad_rho.x = 0.00000000 grad_rho.y = 0.77220000 grad_rho.x = 0.00000000 grad_rho.y = 0.77220000

verifications

fcts, DONE
ordered the update of calculation values like in samy’s code: changes nothing
looked at differences in J update: div_press is different in rust and in c
turns out this difference comes from a difference in grad_rho values, which are swapped and pressure is updated one step earlier in C (grad_rho_IN_C_STEP_-1.x = grad_rho_IN_RUST_STEP_0.y)
it’s also the case for grad_ln_rho, which means it’s a problem with the scalar gradient

SOLUTION?

problem indexing in c

when i

rho[0][99] it gives me 0.8

but for

rho[99][0] it gives me 0.002

wtf?

in c nx = 100 ny = 2

and rho[x_indice][y_indice]

so it’s the number of rows!!!!! fuccccckkkkkkkk

solution

take a convention in the matrixes, and see if i need to invert the position in samy’s code

Matrix organisation choice:

the rows are the y direction
the columns are the x direction

this way, we don’t need to transpose or change the matrixes if we want to plot them

To index our objects (ScalarField, VectorField ect…) we use ScalarField.get_pos(x, y) with x, y the indices in respectively the x and y direction, because it’s usually intuitive to have f(x,y) in physics. But internally, the matrixes are accessible through Matrix[row][column], thus: Matrix[y][x]

another problem

problem

div_press is not good

compare with div_vJ between the two ?

verifications

pressure is good (sure) it’s exactly the same in rust and c at step 0rust results step 0:

GD_pressure.xx [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025], [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025]

GD_pressure.xy [0.0, 0.0, 0.0, 0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0]

GD_pressure.yx [0.0, 0.0, 0.0, 0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0]

GD_pressure.yy = [2.9824461725925957, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 2.9824461725925957, -0.3787229884051399, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, -0.3787229884051399], [2.9824461725925957, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 2.9824461725925957, -0.3787229884051399, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, -0.3787229884051399]

c results step 0:

STEP 0 ———————————(!) pressure print tens_xx [3.578739, -0.047407, -0.047407, -0.047407, 3.578739, 0.217570, 0.002623, 0.002623, 0.002623, 0.217570]

print tens_xy 0

print tens_yx 0

print tens_yy [2.982446, -0.047407, -0.047407, -0.047407, 2.982446, -0.378723, 0.002623, 0.002623, 0.002623, -0.378723]

div_tensor is good (sure)analysis div_tensor gives the good result, the problem lies elsewhere..?proof divergence of tensor: tensor.xx [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025] [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025] tensor.yx=xy=0 tensor.yy [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025] [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025]

should be: (c code) divt.x [-0.262327, -3.589885, 0.0, 3.589885, 0.262327, -3.540354, -0.212797, 0.0, 0.212797, 3.540354], [-0.262327, -3.589885, 0.0, 3.589885, 0.262327, -3.540354, -0.212797, 0.0, 0.212797, 3.540354], divt.y: 0

test result of rust div_tensor application on tensor:

divt.x: [-0.26232748641224496, -3.589884955800003, 0.0, 3.589884955800003, 0.26232748641224496, -3.540354425187758, -0.2127969558000001, 0.0, 0.2127969558000001, 3.540354425187758], [-0.26232748641224496, -3.589884955800003, 0.0, 3.589884955800003, 0.26232748641224496, -3.540354425187758, -0.2127969558000001, 0.0, 0.2127969558000001, 3.540354425187758]], divt.y: 0

GOOD!

pressure value (again)in rust and in clistRUST RIST press.xx [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025], [3.578739012592596, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 3.578739012592596, 0.21756985159486025, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, 0.21756985159486025],

press.yx=xy=0

press.yy [2.9824461725925957, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 2.9824461725925957, -0.3787229884051399, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, -0.3787229884051399], [2.9824461725925957, -0.04740740740740734, -0.04740740740740734, -0.04740740740740734, 2.9824461725925957, -0.3787229884051399, 0.002623431594860167, 0.002623431594860167, 0.002623431594860167, -0.3787229884051399]

C

print tens_xx [3.578739, -0.047407, -0.047407, -0.047407, 3.578739, 0.217570, 0.002623, 0.002623, 0.002623, 0.217570], [3.578739, -0.047407, -0.047407, -0.047407, 3.578739, 0.217570, 0.002623, 0.002623, 0.002623, 0.217570]

print tens_yy [2.982446, -0.047407, -0.047407, -0.047407, 2.982446, -0.378723, 0.002623, 0.002623, 0.002623, -0.378723]

print i: 0, j: 0, x -0.262327 i: 1, j: 0, x -3.589885 i: 2, j: 0, x 0.000000 i: 3, j: 0, x 3.589885 i: 4, j: 0, x 0.262327 i: 5, j: 0, x -3.540354 i: 6, j: 0, x -0.212797 i: 7, j: 0, x 0.000000 i: 8, j: 0, x 0.212797 i: 9, j: 0, x 3.540354div press in rust div_press.x is not the same in the two lines! [0.0, -2.3619677483111134, 0.031288888888888845, 0.031288888888888845, 0.031288888888888845, -2.3619677483111134, -0.14359610205260778, -0.0017314648526077104, -0.0017314648526077104, 2.3602362834585056], [-0.08744249547074831, -3.5585960669111145, 0.031288888888888845, 1.2279172074888898, 0.11873138435963716, -3.542085890040366, -0.21452842065260783, -0.0017314648526077104, 0.06920085374739232, 3.540354425187758]div press in c div_press.x -0.262327, -3.589885, 0.000000, 3.589885, 0.262327, -3.540354, -0.212797, 0.000000, 0.212797, 3.540354

div_press.y 0

printing stuff near the problematic line

the real problem

we need to update entirely independant stuff, and then each the other…

in my previous code:

bad code

i,j loop / ( update t1 ( update t2 ( update t3(t1, t2) \ end i, j loop

good code

i,j loop / ( update t1 ( update t2 \ end i, j loop

i,j loop / ( update t3(t1, t2) \ end i, j loop

stuff I learned

printing in the first steps

to debug a difference between two simulation programs suposed to do exactly the same thing, print the results (obviously with the same input/initial/computation parameters) at the first steps, to see what term is not equal to the other simulation

do comment code blocks

like the one I did with :

// -------------------------------------------------------
// GD_traceless_grad_v_dyadic_grad_v begin update            
GD_traceless_grad_v_dyadic_grad_v
    .set_pos(SOME CODE);
// GD_traceless_grad_v_dyadic_grad_v end update
// -------------------------------------------------------

make functions to print/write debug info

on csv files for instance

if something really strange happen, check the code logic

for instance in this case, the problem was:

bad code:

  i,j loop
 / 
(  update t1
(  update t2
(  update t3(t1, t2) // t3 depends on t2 and t1
 \ end i, j loop
  -----------------------

good code:

  i,j loop
 / 
(  update t1
(  update t2
 \ end i, j loop
  -----------------------

  i,j loop
 / 
(  update t3(t1, t2)
 \ end i, j loop
  -----------------------

write the code logic when possible

and update it.

Do it in an org file!

for instance:

  defining_all_variables

  fluid_initialisation

  i, j loop
 / 
(  updating GD_rho(i,j)
(  ...
(  ...
 \ 
  -----------------------

dichotomy is key

to find where the problem comes from, dichotomise you research!!!

compare my simulation outputs and samy’s

State “DONE” from “TODO” [2022-12-05 lun. 11:56]

they are almost the same! but they have a little shift in the x direction

idea

idea:

print the parameters just before the c simulation, to see it’s actually all from the input file, and to check if the parameters are actually equal in my rust simu and the c simu
compare results with exactly the same parameters

c code parameters (used)

constants

used:

zeta0
eta0
lambda0
lambda
kB
NY
dy
NX
dx
w
dim
aa
b
m

unused:

Tc
rhom_c
Pc
G
DeBroglie0
inv_m
Jev
hlv
forcex
flux
j_wall_bot
j_wall_top
rho_wall
Tw
nsteps_eq_heat
rho_min

input file

used:

T0
rho_liq
rho_vap
HISTO_SAVE
FINAL_TIME
dt

unused:

STEP_EQ
HISTO_FREQ
T1

used in computations

eta0
m
zeta0
lambda
dx
dy
inv_cv
lambda0
dim
b

pressure

aa
w
kB

fluid initialization

T0
rho_liq
rho_vap

main loop

FINAL_TIME

logging

HISTO_SAVE

(cahn hilliard)

kB
aa
w
lambda0
Tc

understand why it is unstable

State “DONE” from “TODO” [2022-12-03 Sat 16:54]

why ? the time step was not tiny enough, and I needed to put the derivatives like it is done in samy’s code, even though it’s counter intuitive

the grad x/y in samy is not in the same direction as mine

State “DONE” from “TOCHECK” [2022-12-03 Sat 16:59]

in his grad.x it’s a gradient through the row axis (so the y axis!)

update verification

State “DONE” from “TODO” [2022-12-03 Sat 13:00]

things updated

lap rho ok
lap T ok
grad T ok
grad rho ok
ln rho ok
lap v ok
div v ok
grad v ok
grad ln rho ok
grad div v ok
traceless grad v ok
vJ ok
div vJ ok
v scalar grad T ok
traceless grad v dyadic grad ok
v scalar grad ln rho ok
pressure ok
grad ln rho traceless grad v ok
grad ln rho scalar grad T ok
div press ok

conservation eq verification

State “DONE” from “TODO” [2022-12-03 Sat 16:59]

momentum

State “DONE” from [2022-12-03 Sat 13:12]

mass

State “DONE” from [2022-12-03 Sat 13:16]

thermal energy

State “DONE” from [2022-12-03 Sat 16:59]

functions verification

State “DONE” from “TODO” [2022-12-03 Sat 16:59]

shear_viscosity ok
bulk_viscosity ok
dissipative_stress ok
v_nabla_v ok
scal_product ok
tens_product_vec ok
dyadic_product ok
partial_deriv ok
grad_scalar —————————— !
gradient —————————— !
gradient_vector —————————— !
div_vector —————————— !
div_tensor —————————— !
lap_scalar ok
laplacian ok
laplacian_vector ok
grad_div_vel ok
pressure ok

don’t forget to update v with updated J at the end

State “DONE” from “TODO” [2022-12-01 Thu 19:45]

check indexing

State “DONE” from “TODO” [2022-12-01 Thu 19:44]

I think I have problems with indexing:

i, j, x_size, y_size…
the position in [ ][ ] for my stuff…
fix ALL THAT

write conservation equations

State “DONE” from [2022-12-01 Thu 16:21]

compute all the terms needed for the main loop

State “DONE” from [2022-12-01 Thu 15:07]

details

traceless_grad_v

[ ] TensorField2D grad_v
[ ] ScalarField2D div_v

momentum eq

[ ] f32 dt (no computation needed)
[ ] VectorField2D lap_v
[ ] VectorField2D div_vJ = div(vJ)
[ ] VectorField2D grad_div_v = grad_div_function on velocity
[ ] VectorField2D grad_ln_rho_traceless_grad_v = tens_vec_product(traceless_grad_v, grad_ln_rho)
[ ] VectorField2D grad_ln_rho
[ ] VectorField2D div_press = div_tensor(press)
[ ] ScalarField2D div_v

already done

rho

constants

eta0
zeta0

thermal energy

const (not defined) cv = 1.5*kB
1/cv
div_v
traceless_grad_v_grad_v
grad_ln_rho_grad_T
lap_T
v_grad_T

constants

eta0
zeta0
lambda0

mass

ScalarField2D div_v
f32 v_grad_ln_rho = scalar_product(v, grad_ln_rho)
VectorField2D grad_ln_rho

functions:

exp (check gsl lib)

computed

[X] div_press
[X] grad_ln_rho_grad_T renamed grad_ln_rho_scalar_grad_T
[X] grad_ln_rho_traceless_grad_v
[X] press
[X] v_grad_ln_rho
[X] traceless_grad_v_grad_v renamed GD_traceless_grad_v_dyadic_grad_v
[X] v_grad_T renamed GD_v_scal_grad_T
[X] div_vJ
[X] vJ
[X] grad_div_v
[X] grad_ln_rho
[X] grad_v
[X] div_v
[X] lap_v
[X] ln_rho
[X] grad_rho
[X] lap_rho
[X] lap_T
[X] grad_T

dependencies list

name	dependencies
div_press	5
grad_ln_rho_grad_T	5
grad_ln_rho_traceless_grad_v	4
press	4
v_grad_ln_rho	4
traceless_grad_v_grad_v	3
v_grad_T	3
div_vJ	2
vJ	2
grad_div_v	2
grad_ln_rho	2
grad_v	1
div_v	1
lap_v	1
ln_rho	1
grad_rho	1
lap_rho	1
lap_T	1
grad_T	1

define all the terms needed for the main loop

State “DONE” from [2022-12-01 Thu 15:02]

all variables checked

State “DONE” from [2022-12-01 Thu 15:02]

[X] TensorField2D grad_v
[X] ScalarField2D div_v
[X] f32 dt (no computation needed)
[X] VectorField2D lap_v
[X] VectorField2D div_vJ = div(vJ)
[X] VectorField2D grad_div_v = grad_div_function on velocity
[X] VectorField2D grad_ln_rho_traceless_grad_v = tens_vec_product(traceless_grad_v, grad_ln rho)
[X] VectorField2D grad_ln_rho
[X] VectorField2D div_press = div_tensor(press)
[X] ScalarField2D div_v
[X] const (not defined) cv = 1.5*kB
[X] 1/cv
[X] div_v
[X] traceless_grad_v_grad_v = dyadic_product(traceless_grad_v,grad_v) traceless_grad_v_grad_v renamed traceless_grad_v_dyadic_grad_v
[X] grad_ln_rho_grad_T = scal_product(grad_ln_rho,grad_T) grad_ln_rho_grad_T renamed grad_ln_rho_scalar_grad_T
[X] lap_T
[X] v_grad_T renamed v_scal_grad_T

traceless_grad_v

State “DONE” from [2022-11-30 Wed 17:23]

[X] TensorField2D grad_v
[X] ScalarField2D div_v

momentum eq

State “DONE” from [2022-11-30 Wed 20:07]

[X] f32 dt (no computation needed)
[X] VectorField2D lap_v
[X] VectorField2D div_vJ = div(vJ)
[X] VectorField2D grad_div_v = grad_div_function on velocity
[X] VectorField2D grad_ln_rho_traceless_grad_v = tens_vec_product(traceless_grad_v, grad_ln rho)
[X] VectorField2D grad_ln_rho
[X] VectorField2D div_press = div_tensor(press)
[X] ScalarField2D div_v

already done

rho

constants

eta0
zeta0

thermal energy

State “DONE” from [2022-12-01 Thu 15:01]

[X] const (not defined) cv = 1.5*kB
[X] 1/cv
[X] div_v
[X] traceless_grad_v_grad_v = dyadic_product(traceless_grad_v,grad_v) traceless_grad_v_grad_v renamed traceless_grad_v_dyadic_grad_v
[X] grad_ln_rho_grad_T = scal_product(grad_ln_rho,grad_T) grad_ln_rho_grad_T renamed grad_ln_rho_scalar_grad_T
[X] lap_T
[X] v_grad_T renamed v_scal_grad_T

constants

eta0
zeta0
lambda0

ln rho

State “DONE” from [2022-12-01 Thu 15:01]

ScalarField2D div_v
f32 v_grad_ln_rho = scalar_product(v, grad_ln_rho)
VectorField2D grad_ln_rho

functions:

exp (check gsl lib)

performance

before v0.3.0

compare c/rust performances on jemac

State “DONE” from [2022-12-09 ven. 11:55]

org table

timesteps	rust	c
10000	1.363861436	3.829171
50000	6.623272585	19.006658
100000	13.228002208	38.029916
500000	65.768119034	190.004333
1000000	131.536727645	380.368056

python

timesteps [10000, 50000, 100000, 500000, 1000000]

rust [1.363861436, 6.623272585, 13.228002208, 65.768119034, 131.536727645]

c [3.829171, 19.006658, 38.029916, 190.004333, 380.368056]

compare c/rust performances on personal computer

State “DONE” from “TODO” [2022-12-08 Thu 23:17]

org table

nb time steps	rust (time in seconds)	c (time in seconds)
10000	0.636702614	1.605196
50000	3.105668876	8.001673
100000	6.165691205	15.999283
500000	30.774588656	80.354180
1000000	61.782680750	159.692652

python

rust_time_list = np.array([0.636702614, 3.105668876, 6.165691205, 30.774588656, 61.782680750])
c_time_list = np.array([1.605196, 8.001673, 15.999283, 80.354180, 159.692652])
nb_step_list = np.array([10000, 50000, 100000, 500000, 1000000])

change Vector2 to Array2

State “DONE” from “TODO” [2022-12-07 Wed 21:44]

to see

https://docs.rs/ndarray/latest/ndarray/index.html https://docs.rs/ndarray/latest/ndarray/doc/ndarray_for_numpy_users/index.html

reddit answer

my message Hey everyone, I have two questions

context: I’m doing a 2d fluid simulation, and I’m using 2D Vectors to store variables that have a value for each position in my 2D plane. I wrap these Vectors in a structure which has “get_position(i: usize, j: usize)” method implemented, which allows me to change the way these 2D information are stored without changing the code of the simulation accessing these values for computations

Is it way better to use arrays instead of vectors for performance? I’m using vectors because the size of the simulation should be given as an input of the program, and thus is not always known at compile time (but this can be changed is the speed gain is huge)
Is it very slow to wrap my big 2D vectors in structures with setters/getters? this makes my code more flexible but I can also change it if it speeds it very badly

Thank you very much! If you are curious, here is my github, feel free to ask more questions on my project, any help or curiosity is welcomed!

answer by WasserMarder

Vectors of vectors are typically bad for performance because they tend do give band caching performance. My advice is to use Array2 from the ndarray crate. It is also somewhat compatible with numpy if you need python interoperability.

You can implement your TensorField2D ontop of Array3.

Use ArrayView and ArrayViewMut for functions that only need views. This makes them usable independent of data layout i.e. you can create a view with shape (n, m) from any array with shape (…, n, …, m, …).

measuring c execution time

State “DONE” from [2022-12-07 Wed 23:08]

I used this: https://www.geeksforgeeks.org/how-to-measure-time-taken-by-a-program-in-c/

compare execution time/performances rust/c

State “DONE” from [2022-12-07 Wed 22:59]

n_step	c_time (s)	rust_time (s)
100000	17.6	6.7
1000000	172.172902	66.774802945S

257.84% faster on the 1 million test

api

make the program take arguments from input file

State “DONE” from “TODO” [2022-12-07 Wed 16:58]

and choose which arguments to fix and which to get in input!

function that creates a default/template SimConfig

State “DONE” from [2022-11-25 ven. 15:40]

cancelled

CANCELED other error?

Note taken on [2022-12-12 lun. 18:19]
not a problem. I didn’t understand that
rho[i][j]

was the rho table in the i index on x, and j on y
State “CANCELED” from “TODO” [2022-12-12 lun. 18:18]

in samy’s code: div_vector: derivative on x component with respect to rows, so y?

TO FIX ON MY CODE TOO

tohos

TOHO [#2] when getting the info col_dx, because col is y, it’s ambiguous

in lap_scalar and partial_deriv col_dx: dx,

TOHO add log entries, readable in csv

this is human readable, not csv readable

let str_to_append = format!("IMPOSSIBLE COMPUTATION: NEGATIVE LOG\n\
                             step {}, col={}, row={}\n\
                             negative log of rho avoided!\n\
                             rho value: {}\n\
                             ------------\n",
                            &step, &yi, &xi, &rho);
logproblem_counts += 1;
if print_logproblems {
    println!("negative log detected ! check log for more info")};
logfile.write_all(&str_to_append.as_bytes())
    .expect("write failed");
println!("error step {}:\n\
          negative rho: rho = {}", step, rho);

TOHO change gradient: so that it increases in the intuitive sense?

right now, the gradient is positive when the scalar field is bigger in the bottow part of the matrix

or maybe keep it this way to not have overhead with matrices indices

TOHO change tensor.xx .xy .yx .yy notation

if possible?

confusing with accessing x/y positions !

same with vector!

TOHO as usize conversions in fcts

in partial_deriv

let (i, j) = (i as usize, j as usize);

TOHO compilation optimisation

build config, performance book

with cpu-specific optimisations

RUSTFLAGS=”-C target-cpu=native” cargo build –release

with release

cargo build –release cargo run –release

???

-C opt-level

If you are unsure whether -C target-cpu=native is working optimally, compare the output of rustc –print cfg and rustc –print cfg -C target-cpu=native to see if the CPU features are being detected correctly in the latter case. If not, you can use -C target-feature to target specific features.

TOHO separate constants from functions

TOHO colorful outputs

https://docs.rs/colored/latest/colored/

TOHO to verify/test

div tensor

in fcts

laplacian

in fcts

seems to returns 0 a lot…

TOHO Cahn_Hilliard

not done because it seems it’s not used in the code

TOHO move fn create objects in another place than main

TOHO remove things like “as i32” or “as usize” if possible

TOHO harmoniser les fonctions qui prennent VecVec et MyStructs

exemple:

gradient/grad_scalar

mais aussi:

laplacian

TOCHECK traceless_grad_v doesn’t appear in dependencies py program?

TOHO verify if all pub is

is the good way to go in “./src/configfile/cfg_struct.rs”

TOHO remove all the “allow”

unused_variables

TOHO remove & in setter for scalar fields

discussion samy thierry

parallélisation “à la main” Thierry

découper le système en plusieurs cases, et faire les calculs

comparer performances en c et en rust

compiler le rust en –release
faire le parrall avec rayon
faire le openMP sur C (voir ce que samy a envoyé)

tension de surface + dimensions

potentiel chimique à coex et de la pression
puis calculer avec l’intégrale
voir courbe tension superficielle en fct de la température
calcul avec l’intégrale
finir de dimensionner
(densité = g/cm3, correspondant à ce que j’avais pris comme paramètres la dernière fois)

writing

Matrix convention

Matrix organisation choice:

the rows are the y direction
the columns are the x direction

this way, we don’t need to transpose or change the matrixes if we want to plot them

To index our objects (ScalarField, VectorField ect…) we use ScalarField.get_pos(x, y) with x, y the indices in respectively the x and y direction, because it’s usually intuitive to have f(x,y) in physics. But internally, the matrixes are accessible through Matrix[row][column], thus: Matrix[y][x]

stuff I learned before v0.3.0

printing in the first steps

to debug a difference between two simulation programs suposed to do exactly the same thing, print the results (obviously with the same input/initial/computation parameters) at the first steps, to see what term is not equal to the other simulation

do comment code blocks

like the one I did with :

// -------------------------------------------------------
// GD_traceless_grad_v_dyadic_grad_v begin update            
GD_traceless_grad_v_dyadic_grad_v
    .set_pos(SOME CODE);
// GD_traceless_grad_v_dyadic_grad_v end update
// -------------------------------------------------------

make functions to print/write debug info

on csv files for instance

if something really strange happen, check the code logic

for instance in this case, the problem was:

bad code:

  i,j loop
 / 
(  update t1
(  update t2
(  update t3(t1, t2) // t3 depends on t2 and t1
 \ end i, j loop
  -----------------------

good code:

  i,j loop
 / 
(  update t1
(  update t2
 \ end i, j loop
  -----------------------

  i,j loop
 / 
(  update t3(t1, t2)
 \ end i, j loop
  -----------------------

write the code logic when possible

and update it.

Do it in an org file!

for instance:

  defining_all_variables

  fluid_initialisation

  i, j loop
 / 
(  updating GD_rho(i,j)
(  ...
(  ...
 \ 
  -----------------------

dichotomy is key

to find where the problem comes from, dichotomise you research!!!

Files

project.org

Latest commit

History

project.org

File metadata and controls

rust snippets/crates

2 dimensional array

unknown toml

float comparison

criterion: benchmark/profiling lines of code

parsing cli arguments

creating dir

csv

c++/rust interop

bindgen c / rust interop

code dependencies

flamegraph-rs

rgsl

orga code building

cheatsheet

ndarray rust

todos

improve performance

rayon integration in ndarray

benchmarking

using ArrayView for things needing only views

using Zip! instead of indexing

verify with samy

eta0

zeta0

understand why there is a jump in some simulations

vdw file with jump?

create a module for the simulation, outside of main

verify derivatives direction

fix the grad_scalar test unit

deal with negative log values (:todo:log:)

next

make a function that takes in a SimConfig and runs the simulation

implementation changes

functions

functions

for the gradient func, created partial derivative function

functions args

v_nabla_v

function names

tens_vec_product

rust

c

renaming

done

done

plots

adapt simplot to plot the results of my simulation on rust

physics loop construction

fix the derive problem

before changing rows/cols

hints

problem in pressure (coming from gradient scalar)

problem in scalar gradient

prints

verifications

SOLUTION?

problem indexing in c

solution

another problem

problem

verifications

the real problem

stuff I learned

printing in the first steps

do comment code blocks

make functions to print/write debug info

if something really strange happen, check the code logic

write the code logic when possible

dichotomy is key

compare my simulation outputs and samy’s

idea

c code parameters (used)

constants