Skip to content

Commit

Permalink
Including the active tension as an additional state variable in the T…
Browse files Browse the repository at this point in the history 
…oRORd_Land model
  • Loading branch information
@bergolho
bergolho committed Sep 9, 2024
1 parent 70b0e5e commit d09b20a
Show file tree
Hide file tree
Showing 8 changed files with 57 additions and 44 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
58 changes: 26 additions & 32 deletions src/main_clip_mesh.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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) {
Expand All @@ -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);
}

}
Expand All @@ -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);
Expand Down
7 changes: 6 additions & 1 deletion src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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"
}
Expand Down Expand Up @@ -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"
}
1 change: 1 addition & 0 deletions src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.common.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -828,3 +828,4 @@ rDY_[45] = dCa_TRPN;
rDY_[46] = dTmBlocked;
rDY_[47] = dZETAS;
rDY_[48] = dZETAW;
rDY_[49] = Ta;
19 changes: 15 additions & 4 deletions src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
}
}
}
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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_);
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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_);
Expand Down Expand Up @@ -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"
Expand Down
12 changes: 6 additions & 6 deletions src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,15 +10,15 @@
#include "../../extra_data_library/helper_functions.h"
#include <stdlib.h>

#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]

Expand Down Expand Up @@ -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]

Expand Down
1 change: 1 addition & 0 deletions src/models_library/ToROrd/ToRORd_Land_mixed_endo_mid_epi_RL.common.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -884,3 +884,4 @@ rDY_[45] = dCa_TRPN;
rDY_[46] = dTmBlocked;
rDY_[47] = dZETAS;
rDY_[48] = dZETAW;
rDY_[49] = Ta;
2 changes: 1 addition & 1 deletion src/monodomain/monodomain_solver.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down

0 comments on commit d09b20a

Please sign in to comment.