From d09b20a4477bc599d4c94717983c403258d005de Mon Sep 17 00:00:00 2001 From: bergolho Date: Mon, 9 Sep 2024 12:43:10 +0100 Subject: [PATCH] Including the active tension as an additional state variable in the ToRORd_Land model --- .gitignore | 1 + src/main_clip_mesh.c | 58 +++++++++---------- .../ToROrd/ToRORd_Land_mixed_endo_mid_epi.c | 7 ++- .../ToRORd_Land_mixed_endo_mid_epi.common.c | 1 + .../ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu | 19 ++++-- .../ToROrd/ToRORd_Land_mixed_endo_mid_epi.h | 12 ++-- ...ToRORd_Land_mixed_endo_mid_epi_RL.common.c | 1 + src/monodomain/monodomain_solver.c | 2 +- 8 files changed, 57 insertions(+), 44 deletions(-) diff --git a/.gitignore b/.gitignore index a7d422ec..3548927b 100644 --- a/.gitignore +++ b/.gitignore @@ -169,6 +169,7 @@ scripts/tuneCV/outputs/* private_models/ src/*/custom_functions.c src/*/custom_functions_BACKUP.c +src/*/custom_functions_luana.c src/*/custom_functions.h bsbash/find_functions_with_mpi.sh diff --git a/src/main_clip_mesh.c b/src/main_clip_mesh.c index 77ce2143..41d49096 100644 --- a/src/main_clip_mesh.c +++ b/src/main_clip_mesh.c @@ -87,22 +87,17 @@ static void expand_scar(const char *input, const char *output, int n_rows, char* struct config *domain_config; domain_config = alloc_and_init_config_data(); - char *discretization = "300"; - //char *discretization = "500"; + char *discretization = "500"; - shput_dup_value(domain_config->config_data, "start_discretization", strdup(discretization)); - shput_dup_value(domain_config->config_data, "maximum_discretization", strdup(discretization)); + shput_dup_value(domain_config->config_data, "original_discretization", strdup(discretization)); + shput_dup_value(domain_config->config_data, "desired_discretization", strdup(discretization)); shput_dup_value(domain_config->config_data, "mesh_file", strdup(input)); - domain_config->main_function_name = strdup("initialize_grid_with_dti_mesh"); - shput_dup_value(domain_config->config_data, "name", "Oxford DTI004 with Transmurality and Fiber orientation"); - - //shput(domain_config->config_data, "side_length_x", strdup(discretization)); - //shput(domain_config->config_data, "side_length_y", strdup(discretization)); - //shput(domain_config->config_data, "side_length_z", strdup(discretization)); + domain_config->main_function_name = strdup("initialize_grid_with_dti_mesh_with_dense_and_sparse_regions_twave"); + shput_dup_value(domain_config->config_data, "name", "Oxford DTI004 with Transmurality, Fiber orientation and Dense/Sparse regions"); + shput(domain_config->config_data, "num_volumes", strdup(n_volumes)); - //shput(domain_config->config_data, "num_extra_fields", "5"); - + init_config_functions(domain_config, "./shared_libs/libdefault_domains.so", "domain"); int success = ((set_spatial_domain_fn*)domain_config->main_function)(domain_config, grid); @@ -117,15 +112,16 @@ static void expand_scar(const char *input, const char *output, int n_rows, char* real_cpu dx, dy, dz; real_cpu min_x, max_x, min_y, max_y, min_z, max_z; real *f, *s, *n; - real_cpu transmurality, base_apex_heterogeneity, apicobasal; - uint32_t transmurality_labels; + real_cpu transmurality, transmurality_labels, base_apex; + real_cpu dense_sparse_tag, healthy, sf_IKs; - min_x=76803.6; - min_y=51176.3; - min_z=8740.9; - max_x=96803.6; - max_y=71176.3; - max_z=28740.9; + // Box bounds + min_x=67108; + min_y=52945; + min_z=14308; + max_x=77108; + max_y=62945; + max_z=24308; FILE *out = fopen(output, "w"); FOR_EACH_CELL(grid) { @@ -141,22 +137,20 @@ static void expand_scar(const char *input, const char *output, int n_rows, char* extra_data = (struct dti_mesh_info *)cell->mesh_extra_info; transmurality = extra_data->transmurality; - base_apex_heterogeneity = extra_data->base_apex_heterogeneity; - apicobasal = extra_data->apicobasal; transmurality_labels = extra_data->dti_transmurality_labels; - - // Change the FAST_ENDO tag to ENDO - if (transmurality_labels == 3) { - transmurality_labels = 0; - } + base_apex = extra_data->base_apex_heterogeneity; + dense_sparse_tag = extra_data->fast_endo; + healthy = 1.0; + sf_IKs = extra_data->sf_IKs; ignore_cell = cell->center.x < min_x || cell->center.x > max_x || cell->center.y < min_y || cell->center.y > max_y || cell->center.z < min_z || cell->center.z > max_z; if(!ignore_cell) { - fprintf(out, "%g,%g,%g,%g,%g,%g,%g,%g,%g,%u,%g,%g,%g,%g,%g,%g,%g,%g,%g\n", cell->center.x, cell->center.y, cell->center.z, dx, dy, dz, \ - transmurality, base_apex_heterogeneity, apicobasal, \ - transmurality_labels, \ - f[0], f[1], f[2], s[0], s[1], s[2], n[0], n[1], n[2]); + fprintf(out, "%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g\n", cell->center.x, cell->center.y, cell->center.z, dx, dy, dz, \ + transmurality, transmurality_labels, base_apex, \ + dense_sparse_tag, healthy, \ + f[0], f[1], f[2], s[0], s[1], s[2], n[0], n[1], n[2], \ + sf_IKs); } } @@ -183,7 +177,7 @@ int main(int argc, char **argv) { } if(!output) { - output = "expanded_mesh.alg"; + output = "clipped_mesh.alg"; } expand_scar(input, output, options->n_rows, options->n_volumes); diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c index bf32adbd..44dbe570 100644 --- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c @@ -117,6 +117,7 @@ SET_ODE_INITIAL_CONDITIONS_CPU(set_model_initial_conditions_cpu) { sv[46] = 9.993734e-01; sv[47] = 0.000000e+00; sv[48] = 0.000000e+00; + sv[49] = 0.000000e+00; // Default initial conditions for MID cell (from original Matlab script) //sv[0] = -8.953800e+01; @@ -265,7 +266,7 @@ SOLVE_MODEL_ODES(solve_model_odes_cpu) { else { // Default: initialize all current modifiers for (uint32_t i = 0; i < num_extra_parameters; i++) { - if (i == 9) + if (i == 10) extra_par[i] = 0.0; else extra_par[i] = 1.0; @@ -368,6 +369,7 @@ void solve_model_ode_cpu(real dt, real *sv, real stim_current, real transmuralit SOLVE_EQUATION_EULER_CPU(46); // TmBlocked SOLVE_EQUATION_EULER_CPU(47); // ZETAS SOLVE_EQUATION_EULER_CPU(48); // ZETAW + SOLVE_EQUATION_CONSTANT_CPU(49); // Ta } void solve_forward_euler_cpu_adpt(real *sv, real stim_curr, real transmurality, real final_time, int sv_id, struct ode_solver *solver, real const *extra_params) { @@ -583,6 +585,7 @@ void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, real transm real TmBlocked = sv[46]; real ZETAS = sv[47]; real ZETAW = sv[48]; + real TA = sv[49]; #include "ToRORd_Land_mixed_endo_mid_epi.common.c" } @@ -678,5 +681,7 @@ void RHS_RL_cpu(real *a_, real *b_, const real *sv, real *rDY_, real stim_curren real ZETAS = sv[47]; real ZETAW = sv[48]; + real TA = sv[49]; + #include "ToRORd_Land_mixed_endo_mid_epi_RL.common.c" } diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c index 065a8c08..e371e395 100644 --- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c @@ -828,3 +828,4 @@ rDY_[45] = dCa_TRPN; rDY_[46] = dTmBlocked; rDY_[47] = dZETAS; rDY_[48] = dZETAW; +rDY_[49] = Ta; diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu index 5b0f64c7..d048e054 100644 --- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu @@ -57,12 +57,13 @@ __global__ void kernel_set_model_initial_conditions(real *sv, int num_volumes, s *((real * )((char *) sv + pitch * 46) + threadID) = 9.993734e-01; *((real * )((char *) sv + pitch * 47) + threadID) = 0.000000e+00; *((real * )((char *) sv + pitch * 48) + threadID) = 0.000000e+00; + *((real * )((char *) sv + pitch * 49) + threadID) = 0.000000e+00; } if(use_adpt_dt) { - *((real *)((char *)sv + pitch * 49) + threadID) = min_dt; // dt - *((real *)((char *)sv + pitch * 50) + threadID) = 0.0; // time_new - *((real *)((char *)sv + pitch * 51) + threadID) = 0.0; // previous dt + *((real *)((char *)sv + pitch * 50) + threadID) = min_dt; // dt + *((real *)((char *)sv + pitch * 51) + threadID) = 0.0; // time_new + *((real *)((char *)sv + pitch * 52) + threadID) = 0.0; // previous dt } } } @@ -230,7 +231,7 @@ extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) { else { // Default: initialize all current modifiers for (uint32_t i = 0; i < num_extra_parameters; i++) { - if (i == 9) + if (i == 10) extra_par[i] = 0.0; else extra_par[i] = 1.0; @@ -334,6 +335,7 @@ __global__ void solve_gpu(real cur_time, real dt, real *sv, real *stim_currents, SOLVE_EQUATION_EULER_GPU(46); // TmBlocked SOLVE_EQUATION_EULER_GPU(47); // ZETAS SOLVE_EQUATION_EULER_GPU(48); // ZETAW + SOLVE_EQUATION_CONSTANT_GPU(49); // Ta } } else { solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], 0.0, extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); @@ -416,6 +418,7 @@ __global__ void solve_endo_mid_epi_gpu(real cur_time, real dt, real *sv, real *s SOLVE_EQUATION_EULER_GPU(46); // TmBlocked SOLVE_EQUATION_EULER_GPU(47); // ZETAS SOLVE_EQUATION_EULER_GPU(48); // ZETAW + SOLVE_EQUATION_CONSTANT_GPU(49); // Ta } } else { solve_forward_euler_gpu_adpt(sv, stim_currents[threadID], transmurality[threadID], extra_params, cur_time + max_dt, sv_id, pitch, abstol, reltol, dt, max_dt); @@ -645,6 +648,7 @@ inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, \ real ZETAS; real ZETAW; + real TA; if (use_adpt_dt) { v = sv[0]; nai = sv[1]; @@ -701,6 +705,7 @@ inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, \ TmBlocked = sv[46]; ZETAS = sv[47]; ZETAW = sv[48]; + TA = sv[49]; } else { v = *((real *)((char *)sv + pitch * 0) + threadID_); nai = *((real *)((char *)sv + pitch * 1) + threadID_); @@ -758,6 +763,8 @@ inline __device__ void RHS_gpu(real *sv, real *rDY_, real stim_current, \ TmBlocked = (*((real *)((char *)sv + pitch * 46) + threadID_)); ZETAS = (*((real *)((char *)sv + pitch * 47) + threadID_)); ZETAW = (*((real *)((char *)sv + pitch * 48) + threadID_)); + + TA = (*((real *)((char *)sv + pitch * 49) + threadID_)); } #include "ToRORd_Land_mixed_endo_mid_epi.common.c" @@ -855,6 +862,8 @@ inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real real ZETAS; real ZETAW; + real TA; + if (use_adpt_dt) { v = sv[0]; nai = sv[1]; @@ -911,6 +920,7 @@ inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real TmBlocked = sv[46]; ZETAS = sv[47]; ZETAW = sv[48]; + TA = sv[49]; } else { v = *((real *)((char *)sv + pitch * 0) + threadID_); nai = *((real *)((char *)sv + pitch * 1) + threadID_); @@ -968,6 +978,7 @@ inline __device__ void RHS_RL_gpu(real *a_, real *b_, real *sv, real *rDY_, real TmBlocked = (*((real *)((char *)sv + pitch * 46) + threadID_)); ZETAS = (*((real *)((char *)sv + pitch * 47) + threadID_)); ZETAW = (*((real *)((char *)sv + pitch * 48) + threadID_)); + TA = (*((real *)((char *)sv + pitch * 49) + threadID_)); } #include "ToRORd_Land_mixed_endo_mid_epi_RL.common.c" diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h index a2a874e3..4205d8f6 100644 --- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h @@ -10,15 +10,15 @@ #include "../../extra_data_library/helper_functions.h" #include -#define NEQ 49 +#define NEQ 50 #define INITIAL_V (-8.863699e+01) -#define ENDO 0.0 -#define MID 1.0 -#define EPI 2.0 +#define ENDO 1.0 +#define MID 2.0 +#define EPI 3.0 // CPU macros -#define SOLVE_EQUATION_CONSTANT_CPU(id) sv[id] = rY[id] +#define SOLVE_EQUATION_CONSTANT_CPU(id) sv[id] = rDY[id] #define SOLVE_EQUATION_EULER_CPU(id) sv[id] = dt * rDY[id] + rY[id] @@ -49,7 +49,7 @@ greatestError = (auxError > greatestError) ? auxError : greatestError // GPU macros -#define SOLVE_EQUATION_CONSTANT_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = *((real *)((char *)sv + pitch * id) + sv_id) +#define SOLVE_EQUATION_CONSTANT_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = rDY[id] #define SOLVE_EQUATION_EULER_GPU(id) *((real *)((char *)sv + pitch * id) + sv_id) = *((real *)((char *)sv + pitch * id) + sv_id) + dt * rDY[id] diff --git a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c index 66cb582a..331c446f 100644 --- a/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c +++ b/src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c @@ -884,3 +884,4 @@ rDY_[45] = dCa_TRPN; rDY_[46] = dTmBlocked; rDY_[47] = dZETAS; rDY_[48] = dZETAW; +rDY_[49] = Ta; diff --git a/src/monodomain/monodomain_solver.c b/src/monodomain/monodomain_solver.c index 9b0778ad..f1d885fe 100644 --- a/src/monodomain/monodomain_solver.c +++ b/src/monodomain/monodomain_solver.c @@ -1599,7 +1599,7 @@ void write_pmj_delay(struct grid *the_grid, struct config *config, struct termin // fprintf(output_file,"ERROR! Probably there was a block on the anterograde direction!\n"); has_block = true; cur_pulse = 0; - // return; + return; } mean_tissue_lat += activation_times_array_tissue[cur_pulse]; if(activation_times_array_tissue[cur_pulse] < min_tissue_lat) {