Merge pull request #76 from bergolho/master

Changes to reproduce T-wave paper and minor fixes
rsachetto · Oct 17, 2024 · 5fc8906 · 5fc8906
2 parents a7fc88a + 1c6761b
commit 5fc8906
Show file tree

Hide file tree

Showing 23 changed files with 1,647 additions and 741 deletions.
diff --git a/src/3dparty/ini_parser/ini.h b/src/3dparty/ini_parser/ini.h
@@ -100,7 +100,7 @@ int ini_parse_string(const char* string, ini_handler handler, void* user);
 
 /* Maximum line length for any line in INI file. */
 #ifndef INI_MAX_LINE
-#define INI_MAX_LINE 200
+#define INI_MAX_LINE 2048
 #endif
 
 #ifdef __cplusplus

diff --git a/src/config/ecg_config.h b/src/config/ecg_config.h
@@ -11,7 +11,7 @@
 #include "../monodomain/monodomain_solver.h"
 #include "config_common.h"
 
-#define CALC_ECG(name) void name(struct time_info *time_info, struct config *config, struct grid *the_grid)
+#define CALC_ECG(name) void name(struct time_info *time_info, struct config *config, struct grid *the_grid, char *output_dir)
 typedef CALC_ECG(calc_ecg_fn);
 
 #define INIT_CALC_ECG(name) void name(struct config *config, struct ode_solver *the_ode_solver, struct grid *the_grid, char *output_dir)

diff --git a/src/domains_library/mesh_info_data.h b/src/domains_library/mesh_info_data.h
@@ -12,17 +12,12 @@ struct default_fibrotic_mesh_info {
 
 // TODO: Move this struct and its macros to "custom_mesh_info_data.h"!
 struct dti_mesh_info {
-	enum dti_transmurality_labels {
-		DTI_FAST_ENDO,
-		DTI_ENDO,
-		DTI_MID,
-		DTI_EPI
-	} dti_transmurality_labels;
-	float transmurality;
-    float apicobasal;
-	float base_apex_heterogeneity;
-    float sf_IKs;
-    int fast_endo;
+		float dti_transmurality_labels;
+		float transmurality;
+    	float apicobasal;
+		float base_apex_heterogeneity;
+    	float sf_IKs;
+    	int fast_endo;
 };
 
 #define FIBROTIC_INFO(grid_cell) (struct default_fibrotic_mesh_info *)(grid_cell)->mesh_extra_info

diff --git a/src/ecg_library/ecg.c b/src/ecg_library/ecg.c
@@ -482,12 +482,12 @@ CALC_ECG(pseudo_bidomain) {
     GET_PARAMETER_BOOLEAN_VALUE_OR_USE_DEFAULT(gpu, config, "use_gpu");
     if(gpu) {
 #ifdef COMPILE_CUDA
-        pseudo_bidomain_gpu(time_info, config, the_grid);
+        pseudo_bidomain_gpu(time_info, config, the_grid, output_dir);
 #else
-        pseudo_bidomain_cpu(time_info, config, the_grid);
+        pseudo_bidomain_cpu(time_info, config, the_grid, output_dir);
 #endif
     } else {
-        pseudo_bidomain_cpu(time_info, config, the_grid);
+        pseudo_bidomain_cpu(time_info, config, the_grid, output_dir);
     }
 }
 
@@ -504,3 +504,316 @@ END_CALC_ECG(end_pseudo_bidomain) {
         end_pseudo_bidomain_cpu(config);
     }
 }
+
+INIT_CALC_ECG(init_pseudo_bidomain_with_diffusive_current_cpu) {
+    config->persistent_data = CALLOC_ONE_TYPE(struct pseudo_bidomain_persistent_data);
+
+    uint32_t nlen_output_dir = strlen(output_dir);
+    char *filename = (char*)malloc(sizeof(char)*nlen_output_dir+10);
+    sprintf(filename, "%s/ecg.txt", output_dir);
+
+    char *dir = get_dir_from_path(filename);
+    create_dir(dir);
+    free(dir);
+
+    PSEUDO_BIDOMAIN_DATA->output_file = fopen(filename, "w");
+
+    if(PSEUDO_BIDOMAIN_DATA->output_file == NULL) {
+        log_error_and_exit("init_pseudo_bidomain - Unable to open file %s!\n", filename);
+    }
+
+    real_cpu sigma_b = 1.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_b, config, "sigma_b");
+    uint32_t diff_curr_rate = 100;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(uint32_t, diff_curr_rate, config, "diff_curr_rate");
+    PSEUDO_BIDOMAIN_DATA->diff_curr_rate = diff_curr_rate;
+    real_cpu diff_curr_max_time = 100.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, diff_curr_max_time, config, "diff_curr_max_time");
+    PSEUDO_BIDOMAIN_DATA->diff_curr_max_time = diff_curr_max_time;
+
+    if(sigma_b == 0.0) {
+        log_error_and_exit("init_pseudo_bidomain - sigma_b can't be 0!\n");
+    }
+
+    PSEUDO_BIDOMAIN_DATA->scale_factor = 1.0 / (4.0 * M_PI * sigma_b);
+
+    get_leads(config);
+    PSEUDO_BIDOMAIN_DATA->n_leads = arrlen(PSEUDO_BIDOMAIN_DATA->leads);
+
+    uint32_t n_active = the_grid->num_active_cells;
+    struct cell_node **ac = the_grid->active_cells;
+
+    PSEUDO_BIDOMAIN_DATA->distances = MALLOC_ARRAY_OF_TYPE(real, PSEUDO_BIDOMAIN_DATA->n_leads * n_active);
+
+    PSEUDO_BIDOMAIN_DATA->beta_im = MALLOC_ARRAY_OF_TYPE(real, n_active);
+
+    // calc the distances from each volume to each electrode (r)
+    for(uint32_t i = 0; i < PSEUDO_BIDOMAIN_DATA->n_leads; i++) {
+
+        struct point_3d lead = PSEUDO_BIDOMAIN_DATA->leads[i];
+
+        OMP(parallel for)
+        for(int j = 0; j < n_active; j++) {
+            uint32_t index = i * n_active + j;
+            struct point_3d center = ac[j]->center;
+            PSEUDO_BIDOMAIN_DATA->distances[index] = EUCLIDIAN_DISTANCE(lead, center);
+        }
+    }
+
+    assembly_divergent(the_grid);
+
+    free(filename);
+}
+
+CALC_ECG(pseudo_bidomain_with_diffusive_current_cpu) {
+    // use the equation described in https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3378475/#FD7
+    uint32_t n_active = the_grid->num_active_cells;
+    struct cell_node **ac = the_grid->active_cells;
+
+    OMP(parallel for)
+    for(uint32_t i = 0; i < n_active; i++) {
+        struct element *cell_elements = ac[i]->elements;
+        size_t max_el = arrlen(cell_elements);
+
+        struct point_3d d = ac[i]->discretization;
+        real_cpu volume = d.x * d.y * d.z;
+
+        PSEUDO_BIDOMAIN_DATA->beta_im[i] = 0;
+
+        for(size_t el = 0; el < max_el; el++) {
+            PSEUDO_BIDOMAIN_DATA->beta_im[i] += cell_elements[el].value_ecg * cell_elements[el].cell->v / volume;
+        }
+    }
+
+    // DIFFUSIVE CURRENT PLOT
+    if (time_info->iteration % PSEUDO_BIDOMAIN_DATA->diff_curr_rate == 0 && time_info->current_t <= PSEUDO_BIDOMAIN_DATA->diff_curr_max_time) {
+        char *filename_diffusive_current[200];
+        sprintf(filename_diffusive_current, "%s/diffusive_current_t_%g_ms.txt", output_dir, time_info->current_t);
+        FILE *output_diffusive_current = fopen(filename_diffusive_current, "w");
+        for(uint32_t i = 0; i < n_active; i++) {
+            fprintf(output_diffusive_current, "%g\n", PSEUDO_BIDOMAIN_DATA->beta_im[i]);
+        }
+        fclose(output_diffusive_current);
+    }
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", time_info->current_t);
+
+    for(uint32_t i = 0; i < PSEUDO_BIDOMAIN_DATA->n_leads; i++) {
+        real_cpu local_sum = 0.0;
+
+        OMP(parallel for reduction(+:local_sum))
+        for(uint32_t j = 0; j < n_active; j++) {
+            struct point_3d d = ac[j]->discretization;
+            real_cpu volume = d.x * d.y * d.z;
+
+            uint32_t index = i * n_active + j;
+            local_sum += ((PSEUDO_BIDOMAIN_DATA->beta_im[j] / PSEUDO_BIDOMAIN_DATA->distances[index])) * volume;
+        }
+
+        fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", -PSEUDO_BIDOMAIN_DATA->scale_factor * local_sum);
+    }
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "\n");
+}
+
+END_CALC_ECG(end_pseudo_bidomain_with_diffusive_current_cpu) {
+    fclose(PSEUDO_BIDOMAIN_DATA->output_file);
+    free(PSEUDO_BIDOMAIN_DATA->beta_im);
+    arrfree(PSEUDO_BIDOMAIN_DATA->leads);
+}
+
+#ifdef COMPILE_CUDA
+
+INIT_CALC_ECG(init_pseudo_bidomain_with_diffusive_current_gpu) {
+
+    init_pseudo_bidomain_with_diffusive_current_cpu(config, NULL, the_grid, output_dir);
+
+    if(!the_ode_solver->gpu) {
+        log_warn("The current implementation of pseudo_bidomain_gpu only works when the odes are also being solved using the GPU! Falling back to CPU version\n");
+        shput(config->config_data, "use_gpu", "false");
+        return;
+    }
+
+    //This is allocated when using the CPU code, but we do not need it in the gpu version
+    free(PSEUDO_BIDOMAIN_DATA->beta_im);
+
+    check_cublas_error(cusparseCreate(&(PSEUDO_BIDOMAIN_DATA->cusparseHandle)));
+    check_cublas_error(cublasCreate(&(PSEUDO_BIDOMAIN_DATA->cublasHandle)));
+
+    int_array I = NULL, J = NULL;
+    f32_array val = NULL;
+
+    grid_to_csr_for_ecg(the_grid, &val, &I, &J, false, true);
+
+    int nz = arrlen(val);
+    uint32_t N = the_grid->num_active_cells;
+
+    // DIFFUSIVE CURRENT
+    // Allocate memory for the CPU beta_im
+    PSEUDO_BIDOMAIN_DATA->beta_im_cpu = MALLOC_ARRAY_OF_TYPE(real, N);
+
+    // We need to do that because val is originally a float value. Maybe we should make it real.
+    real *new_val = NULL;
+    arrsetlen(new_val, nz);
+    for(int i = 0; i < nz; i++) {
+        new_val[i] = (real)val[i];
+    }
+
+    PSEUDO_BIDOMAIN_DATA->nz = nz;
+
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_col), nz * sizeof(int)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_row), (N + 1) * sizeof(int)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_val), nz * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->beta_im), N * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_distances), PSEUDO_BIDOMAIN_DATA->n_leads * N * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->d_volumes), N * sizeof(real)));
+    check_cuda_error(cudaMalloc((void **)&(PSEUDO_BIDOMAIN_DATA->tmp_data), N * sizeof(real)));
+
+#if CUBLAS_VER_MAJOR >= 11
+    check_cuda_error(cusparseCreateCsr(&(PSEUDO_BIDOMAIN_DATA->matA), N, N, nz, PSEUDO_BIDOMAIN_DATA->d_row, PSEUDO_BIDOMAIN_DATA->d_col,
+                                       PSEUDO_BIDOMAIN_DATA->d_val, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUBLAS_SIZE));
+    check_cuda_error(cusparseCreateDnVec(&(PSEUDO_BIDOMAIN_DATA->vec_beta_im), N, PSEUDO_BIDOMAIN_DATA->beta_im, CUBLAS_SIZE));
+    check_cuda_error(cusparseCreateDnVec(&(PSEUDO_BIDOMAIN_DATA->vec_vm), N, the_ode_solver->sv, CUBLAS_SIZE));
+#else
+    check_cuda_error((cudaError_t)cusparseCreateMatDescr(&(PSEUDO_BIDOMAIN_DATA->descr)));
+    cusparseSetMatType(PSEUDO_BIDOMAIN_DATA->descr, CUSPARSE_MATRIX_TYPE_GENERAL);
+    cusparseSetMatIndexBase(PSEUDO_BIDOMAIN_DATA->descr, CUSPARSE_INDEX_BASE_ZERO);
+    PSEUDO_BIDOMAIN_DATA->local_sv = the_ode_solver->sv;
+#endif
+
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_col, J, nz * sizeof(int), cudaMemcpyHostToDevice));        // JA
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_row, I, (N + 1) * sizeof(int), cudaMemcpyHostToDevice));   // IA
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_val, new_val, nz * sizeof(real), cudaMemcpyHostToDevice)); // A
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_distances, PSEUDO_BIDOMAIN_DATA->distances, PSEUDO_BIDOMAIN_DATA->n_leads * N * sizeof(real),
+                                cudaMemcpyHostToDevice));
+
+    PSEUDO_BIDOMAIN_DATA->volumes = MALLOC_ARRAY_OF_TYPE(real, N);
+
+    struct cell_node **ac = the_grid->active_cells;
+    OMP(parallel for)
+    for(int i = 0; i < N; i += 1) {
+        struct point_3d d = ac[i]->discretization;
+        PSEUDO_BIDOMAIN_DATA->volumes[i] = d.x * d.y * d.z;
+    }
+
+    check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->d_volumes, PSEUDO_BIDOMAIN_DATA->volumes, N * sizeof(real), cudaMemcpyHostToDevice));
+
+#if CUSPARSE_VER_MAJOR >= 11
+    real alpha = 1.0;
+    real beta = 0.0;
+
+    check_cuda_error(cusparseSpMV_bufferSize(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, PSEUDO_BIDOMAIN_DATA->matA,
+                                             PSEUDO_BIDOMAIN_DATA->vec_vm, &beta, PSEUDO_BIDOMAIN_DATA->vec_beta_im, CUBLAS_SIZE, CUSPARSE_ALG,
+                                             &(PSEUDO_BIDOMAIN_DATA->bufferSize)));
+
+    check_cuda_error(cudaMalloc(&(PSEUDO_BIDOMAIN_DATA->buffer), PSEUDO_BIDOMAIN_DATA->bufferSize));
+#endif
+
+    arrfree(I);
+    arrfree(J);
+    arrfree(val);
+    arrfree(new_val);
+
+    free(PSEUDO_BIDOMAIN_DATA->volumes);
+    free(PSEUDO_BIDOMAIN_DATA->distances);
+}
+
+CALC_ECG(pseudo_bidomain_with_diffusive_current_gpu) {
+
+    uint32_t n_active = the_grid->num_active_cells;
+
+    // VM is correct
+    real alpha = 1.0;
+    real beta = 0.0;
+#if CUSPARSE_VER_MAJOR >= 11
+    check_cublas_error(cusparseSpMV(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, PSEUDO_BIDOMAIN_DATA->matA,
+                                    PSEUDO_BIDOMAIN_DATA->vec_vm, &beta, PSEUDO_BIDOMAIN_DATA->vec_beta_im, CUBLAS_SIZE, CUSPARSE_ALG,
+                                    PSEUDO_BIDOMAIN_DATA->buffer));
+#else
+
+#ifdef CELL_MODEL_REAL_DOUBLE
+    cusparseDcsrmv(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, n_active, n_active, PSEUDO_BIDOMAIN_DATA->nz, &alpha,
+                   PSEUDO_BIDOMAIN_DATA->descr, PSEUDO_BIDOMAIN_DATA->d_val, PSEUDO_BIDOMAIN_DATA->d_row, PSEUDO_BIDOMAIN_DATA->d_col,
+                   PSEUDO_BIDOMAIN_DATA->local_sv, &beta, PSEUDO_BIDOMAIN_DATA->beta_im);
+#else
+    cusparseScsrmv(PSEUDO_BIDOMAIN_DATA->cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, n_active, n_active, PSEUDO_BIDOMAIN_DATA->nz, &alpha,
+                   PSEUDO_BIDOMAIN_DATA->descr, PSEUDO_BIDOMAIN_DATA->d_val, PSEUDO_BIDOMAIN_DATA->d_row, PSEUDO_BIDOMAIN_DATA->d_col,
+                   PSEUDO_BIDOMAIN_DATA->local_sv, &beta, PSEUDO_BIDOMAIN_DATA->beta_im);
+#endif
+
+#endif
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", time_info->current_t);
+
+    gpu_vec_div_vec(PSEUDO_BIDOMAIN_DATA->beta_im, PSEUDO_BIDOMAIN_DATA->d_volumes, PSEUDO_BIDOMAIN_DATA->beta_im, n_active);
+
+    // DIFFUSIVE CURRENT
+    if (time_info->iteration % PSEUDO_BIDOMAIN_DATA->diff_curr_rate == 0 && time_info->current_t <= PSEUDO_BIDOMAIN_DATA->diff_curr_max_time) {
+        char *filename_diffusive_current[200];
+        sprintf(filename_diffusive_current, "%s/diffusive_current_t_%g_ms.txt", output_dir, time_info->current_t);
+        FILE *output_diffusive_current = fopen(filename_diffusive_current, "w");
+        check_cuda_error(cudaMemcpy(PSEUDO_BIDOMAIN_DATA->beta_im_cpu, PSEUDO_BIDOMAIN_DATA->beta_im, n_active * sizeof(real), cudaMemcpyDeviceToHost));
+        for(uint32_t i = 0; i < n_active; i++) {
+            fprintf(output_diffusive_current, "%g\n", PSEUDO_BIDOMAIN_DATA->beta_im_cpu[i]);
+        }
+        fclose(output_diffusive_current);
+    }
+
+    for(int i = 0; i < PSEUDO_BIDOMAIN_DATA->n_leads; i++) {
+
+        // beta_im / distance
+        gpu_vec_div_vec(PSEUDO_BIDOMAIN_DATA->beta_im, PSEUDO_BIDOMAIN_DATA->d_distances + n_active * i, PSEUDO_BIDOMAIN_DATA->tmp_data, n_active);
+        real local_sum;
+
+#ifdef CELL_MODEL_REAL_DOUBLE
+        check_cublas_error(
+            cublasDdot(PSEUDO_BIDOMAIN_DATA->cublasHandle, n_active, PSEUDO_BIDOMAIN_DATA->tmp_data, 1, PSEUDO_BIDOMAIN_DATA->d_volumes, 1, &local_sum));
+#else
+        check_cublas_error(
+            cublasSdot(PSEUDO_BIDOMAIN_DATA->cublasHandle, n_active, PSEUDO_BIDOMAIN_DATA->tmp_data, 1, PSEUDO_BIDOMAIN_DATA->d_volumes, 1, &local_sum));
+#endif
+
+        fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "%lf ", -PSEUDO_BIDOMAIN_DATA->scale_factor * local_sum);
+    }
+
+    fprintf(PSEUDO_BIDOMAIN_DATA->output_file, "\n");
+
+}
+
+END_CALC_ECG(end_pseudo_bidomain_with_diffusive_current_gpu) {
+
+    struct pseudo_bidomain_persistent_data *persistent_data = (struct pseudo_bidomain_persistent_data *)config->persistent_data;
+
+    if(!persistent_data)
+        return;
+
+    check_cuda_error((cudaError_t)cusparseDestroy(persistent_data->cusparseHandle));
+    check_cuda_error((cudaError_t)cublasDestroy(persistent_data->cublasHandle));
+
+#if CUBLAS_VER_MAJOR >= 11
+    if(persistent_data->matA) {
+        check_cuda_error(cusparseDestroySpMat(persistent_data->matA));
+    }
+    if(persistent_data->vec_beta_im) {
+        check_cuda_error(cusparseDestroyDnVec(persistent_data->vec_beta_im));
+    }
+    if(persistent_data->vec_vm) {
+        check_cuda_error(cusparseDestroyDnVec(persistent_data->vec_vm));
+    }
+#else
+    check_cuda_error((cudaError_t)cusparseDestroyMatDescr(persistent_data->descr));
+#endif
+    check_cuda_error(cudaFree(persistent_data->d_col));
+    check_cuda_error(cudaFree(persistent_data->d_row));
+    check_cuda_error(cudaFree(persistent_data->d_val));
+    check_cuda_error(cudaFree(persistent_data->beta_im));
+    check_cuda_error(cudaFree(persistent_data->tmp_data));
+    check_cuda_error(cudaFree(persistent_data->d_distances));
+    check_cuda_error(cudaFree(persistent_data->d_volumes));
+
+    free(persistent_data->beta_im_cpu);
+    fclose(persistent_data->output_file);
+
+    free(persistent_data);
+}
+#endif
diff --git a/src/ecg_library/ecg.h b/src/ecg_library/ecg.h
@@ -13,6 +13,8 @@ struct pseudo_bidomain_persistent_data {
     FILE *output_file;
     real_cpu scale_factor;
     uint32_t n_leads;
+    uint32_t diff_curr_rate;
+    real_cpu diff_curr_max_time;
 
 #ifdef COMPILE_CUDA
     cusparseHandle_t cusparseHandle;
@@ -23,6 +25,7 @@ struct pseudo_bidomain_persistent_data {
     real *volumes;
     real *tmp_data;
     real *d_val;
+    real *beta_im_cpu;
 
     size_t bufferSize;
     void *buffer;