diff --git a/ToDo b/ToDo index 7b7cfaa0..e736f162 100644 --- a/ToDo +++ b/ToDo @@ -19,3 +19,12 @@ The logger symbols are only exported to an executable if an static library linked to a shared library uses then. For now this is ok. But I think it will be a pain in future releases. The GPU linear system solver is not working for purkinje-only simulations When the minimum number of PMJs is not reached the solver will be in an infinite loop + +Lucas ToDo list: + +0) Revise all the examples inside the example_configs/ folder. + +1) Purkinje coupling: + - Find a better way to write the PMJ delay into a file for the propagation block cases. + + \ No newline at end of file diff --git a/example_configs/cable_mesh_with_ToRORd_Land_endo_mid_epi.ini b/example_configs/cable_mesh_with_ToRORd_Land_endo.ini similarity index 100% rename from example_configs/cable_mesh_with_ToRORd_Land_endo_mid_epi.ini rename to example_configs/cable_mesh_with_ToRORd_Land_endo.ini diff --git a/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini b/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini index 11100129..b659960b 100644 --- a/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini +++ b/example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini @@ -1,6 +1,10 @@ # ==================================================================== # Author: Lucas Berg # Description: Simple simulation to test the ToRORd_fkatp model in a cable. +# With the [extra_data] section we consider: +# 45% of the cable ENDO +# 25% of the cable MID +# 30% of the cable EPI # When no [extra_data] information is provided all cells are # considered to be control ENDO. # ==================================================================== diff --git a/networks/simple_his_bundle_pmj_location.vtk b/networks/simple_his_bundle_pmj_location.vtk new file mode 100644 index 00000000..47327b75 --- /dev/null +++ b/networks/simple_his_bundle_pmj_location.vtk @@ -0,0 +1,10 @@ +# vtk DataFile Version 4.2 +vtk output +ASCII +DATASET POLYDATA +POINTS 2 float +10000 15000 0 +10000 5000 0 +VERTICES 2 4 +1 0 +1 1 diff --git a/scripts/error_calculator/main.cpp b/scripts/error_calculator/main.cpp index cded0d87..afa50a65 100644 --- a/scripts/error_calculator/main.cpp +++ b/scripts/error_calculator/main.cpp @@ -71,6 +71,10 @@ int main (int argc, char *argv[]) cerr << "[-] ERROR! The number of cells are different" << endl; exit(EXIT_FAILURE); } + else + { + cout << "[+] Cell number are equal " << total_num_cells_1 << " " << total_num_cells_2 << endl; + } vtkSmartPointer array_1 = vtkFloatArray::SafeDownCast(unstructured_grid_1->GetCellData()->GetArray(array_name.c_str())); diff --git a/src/domains_library/domain.c b/src/domains_library/domain.c index 3f6b1187..b296d4f2 100644 --- a/src/domains_library/domain.c +++ b/src/domains_library/domain.c @@ -439,8 +439,11 @@ SET_SPATIAL_DOMAIN(initialize_grid_with_cuboid_and_sphere_fibrotic_mesh_with_con real_cpu phi = 0.0; GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, phi, config, "phi"); - real_cpu plain_center = 0.0; - GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, plain_center, config, "plain_center"); + real_cpu plain_center_x = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, plain_center_x, config, "plain_center_x"); + + real_cpu plain_center_y = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, plain_center_y, config, "plain_center_y"); real_cpu sphere_radius = 0.0; GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sphere_radius, config, "sphere_radius"); @@ -459,7 +462,7 @@ SET_SPATIAL_DOMAIN(initialize_grid_with_cuboid_and_sphere_fibrotic_mesh_with_con //set_square_mesh(config, the_grid); set_cuboid_domain_mesh(the_grid, start_dx, start_dy, start_dz, side_length_x, side_length_y, side_length_z); - set_cuboid_sphere_fibrosis_with_conic_path(the_grid, phi, plain_center, sphere_radius, border_zone_size, border_zone_radius, seed, conic_slope); + set_cuboid_sphere_fibrosis_with_conic_path(the_grid, phi, plain_center_x, plain_center_y, sphere_radius, border_zone_size, border_zone_radius, seed, conic_slope); return 1; -} \ No newline at end of file +} diff --git a/src/domains_library/domain_helpers.c b/src/domains_library/domain_helpers.c index 0ee43c8a..0a79d335 100644 --- a/src/domains_library/domain_helpers.c +++ b/src/domains_library/domain_helpers.c @@ -1010,8 +1010,8 @@ int calc_num_refs(real_cpu start_h, real_cpu desired_h) { return num_refs; } -void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu phi, real_cpu plain_center, real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius, - unsigned fib_seed, real_cpu cone_slope) { +void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu phi, real_cpu plain_center_x, real_cpu plain_center_y, \ + real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius, unsigned fib_seed, real_cpu cone_slope) { log_info("Making %.2lf %% of cells inactive\n", phi * 100.0f); @@ -1028,23 +1028,22 @@ void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu grid_cell = the_grid->first_cell; while(grid_cell != 0) { - //Calcula distância da célula para o centro da malha - real_cpu distance = pow(grid_cell->center.x - plain_center, 2.0) + pow(grid_cell->center.y - plain_center, 2.0); - real_cpu h_distance = abs(grid_cell->center.x - plain_center); + // Calculate distance to the center of the mesh + real_cpu distance = pow(grid_cell->center.x - plain_center_x, 2.0) + pow(grid_cell->center.y - plain_center_y, 2.0); + real_cpu h_distance = abs(grid_cell->center.y - plain_center_y); if(grid_cell->active) { INITIALIZE_FIBROTIC_INFO(grid_cell); if(distance <= bz_radius_2) { - //Dentro da border zone + // Inside border zone if(distance <= sphere_radius_2) { - //dentro da "esfera" - if(h_distance < cone_slope*grid_cell->center.y*plain_center) //abs(cone_slope*grid_cell->center.y) + // Inside the sphere + if(h_distance < cone_slope*abs(the_grid->mesh_side_length.x - grid_cell->center.x)*plain_center_x) //abs(cone_slope*grid_cell->center.y) { + // Inside the cone FIBROTIC(grid_cell) = true; - - //Dentro do cone } else{ grid_cell->active = false; @@ -1069,8 +1068,8 @@ void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu grid_cell->active = false; grid_cell->can_change = false; } else if(BORDER_ZONE(grid_cell)) { - real_cpu distance_from_center = sqrt((grid_cell->center.x - plain_center) * (grid_cell->center.x - plain_center) + - (grid_cell->center.y - plain_center) * (grid_cell->center.y - plain_center)); + real_cpu distance_from_center = sqrt((grid_cell->center.x - plain_center_x) * (grid_cell->center.x - plain_center_x) + + (grid_cell->center.y - plain_center_y) * (grid_cell->center.y - plain_center_y)); distance_from_center = (distance_from_center - sphere_radius) / bz_size; real_cpu phi_local = phi - phi * distance_from_center; real_cpu p = (real_cpu)(rand()) / (RAND_MAX); diff --git a/src/domains_library/domain_helpers.h b/src/domains_library/domain_helpers.h index 14feb812..ebefea48 100644 --- a/src/domains_library/domain_helpers.h +++ b/src/domains_library/domain_helpers.h @@ -56,9 +56,9 @@ uint32_t set_custom_mesh_from_file(struct grid *the_grid, const char *mesh_file, void set_cube_sphere_fibrosis(struct grid *the_grid, real_cpu phi, real_cpu sphere_center[3], real_cpu sphere_radius, unsigned fib_seed); -void set_cuboid_sphere_fibrosis_with_conic_path (struct grid *the_grid, real_cpu phi, real_cpu plain_center, real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius, - unsigned fib_seed, real_cpu cone_slope); - +void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu phi, real_cpu plain_center_x, real_cpu plain_center_y, \ + real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius, unsigned fib_seed, real_cpu cone_slope); + int calc_num_refs(real_cpu start_h, real_cpu desired_h); #endif // MONOALG3D_DOMAIN_HELPERS_H diff --git a/src/domains_library/mesh_info_data.h b/src/domains_library/mesh_info_data.h index 53b9f196..06c74704 100644 --- a/src/domains_library/mesh_info_data.h +++ b/src/domains_library/mesh_info_data.h @@ -10,6 +10,7 @@ struct default_fibrotic_mesh_info { int tissue_type; }; +// TODO: Move this struct and its macros to "custom_mesh_info_data.h"! struct dti_mesh_info { enum dti_transmurality_labels { DTI_FAST_ENDO, @@ -20,6 +21,7 @@ struct dti_mesh_info { float transmurality; float apicobasal; float base_apex_heterogeneity; + float sf_IKs; int fast_endo; }; @@ -45,5 +47,6 @@ struct dti_mesh_info { #define DTI_MESH_BASE_APEX_HETEROGENEITY(grid_cell) (DTI_MESH_INFO(grid_cell))->base_apex_heterogeneity #define DTI_MESH_APICOBASAL(grid_cell) (DTI_MESH_INFO(grid_cell))->apicobasal #define DTI_MESH_FAST_ENDO(grid_cell) (DTI_MESH_INFO(grid_cell))->fast_endo +#define DTI_MESH_SCALE_FACTOR_IKS(grid_cell) (DTI_MESH_INFO(grid_cell))->sf_IKs #endif /* __MESH_INFO_DATA_H */ diff --git a/src/extra_data_library/helper_functions.h b/src/extra_data_library/helper_functions.h index ab812114..eca7df6c 100644 --- a/src/extra_data_library/helper_functions.h +++ b/src/extra_data_library/helper_functions.h @@ -31,11 +31,16 @@ struct extra_data_for_tt3 { real GCaL_multiplicator; real INaCa_multiplicator; real Ikatp_multiplicator; + real GKr_multiplicator; + real GKs_multiplicator; + real alpha_multiplicator; real Vm_modifier; real *fibrosis; real *transmurality; }; + + struct extra_data_for_torord { real INa_Multiplier; real ICaL_Multiplier; @@ -54,6 +59,8 @@ struct extra_data_for_torord { real IClb_Multiplier; real Jrel_Multiplier; real Jup_Multiplier; + real aCaMK_Multiplier; + real taurelp_Multiplier; real *initial_ss_endo; real *initial_ss_epi; real *initial_ss_mid; @@ -120,11 +127,78 @@ struct extra_data_for_trovato { //real *purkinje_tags; }; +struct extra_data_for_torord_land_twave { + real INa_Multiplier; + real INaL_Multiplier; + real INaCa_Multiplier; + real INaK_Multiplier; + real INab_Multiplier; + real Ito_Multiplier; + real IKr_Multiplier; + real IKs_Multiplier; + real IK1_Multiplier; + real IKb_Multiplier; + real IKCa_Multiplier; + real ICaL_Multiplier; + real ICab_Multiplier; + real IpCa_Multiplier; + real ICaCl_Multiplier; + real IClb_Multiplier; + real Jrel_Multiplier; + real Jup_Multiplier; + real aCaMK_Multiplier; + real taurelp_Multiplier; + real *initial_ss_endo; + real *initial_ss_epi; + real *initial_ss_mid; + real *transmurality; + real *sf_IKs; + real *basetoapex; +}; + +struct extra_data_for_torord_general { + real INa_Multiplier; + real INaL_Multiplier; + real INaCa_Multiplier; + real INaK_Multiplier; + real INab_Multiplier; + real Ito_Multiplier; + real IKr_Multiplier; + real IKs_Multiplier; + real IK1_Multiplier; + real IKb_Multiplier; + real IKCa_Multiplier; + real ICaL_Multiplier; + real ICab_Multiplier; + real IpCa_Multiplier; + real ICaCl_Multiplier; + real IClb_Multiplier; + real Jrel_Multiplier; + real Jup_Multiplier; + real aCaMK_Multiplier; + real taurelp_Multiplier; + real *initial_ss_endo; + real *initial_ss_epi; + real *initial_ss_mid; + real *transmurality; + real *sf_IKs; + real *basetoapex; +}; + struct extra_data_for_fibrosis * set_common_schemia_data(struct config *config, uint32_t num_cells); struct extra_data_for_tt3 * set_common_tt3_data (struct config *config, uint32_t num_cells); struct extra_data_for_torord * set_common_torord_data (struct config *config, uint32_t num_cells); struct extra_data_for_torord * set_common_torord_dyncl_data (struct config *config, uint32_t num_cells); struct extra_data_for_torord_land * set_common_torord_Land_data (struct config *config, uint32_t num_cells); +struct extra_data_for_torord_land_twave * set_common_torord_Land_twave_data (struct config *config, uint32_t num_cells); struct extra_data_for_trovato * set_common_trovato_data (struct config *config, uint32_t num_cells); +struct extra_data_for_torord_general * set_common_torord_general (struct config *config, uint32_t num_cells, uint32_t model_id); + +void configure_torord_extra_data_with_user_input (struct extra_data_for_torord_general *extra_data, struct config *config, uint32_t model_id); +void configure_torord_extra_data_initial_conditions (struct extra_data_for_torord_general *extra_data, uint32_t model_id); +void configure_torord_fkatp_initial_conditions (struct extra_data_for_torord_general *extra_data); +void configure_torord_dynCl_initial_conditions (struct extra_data_for_torord_general *extra_data); +void configure_torord_Land_initial_conditions (struct extra_data_for_torord_general *extra_data); + #endif // MONOALG3D_C_EXTRA_DATA_HELPER_FUNCTIONS_H diff --git a/src/monodomain/monodomain_solver.c b/src/monodomain/monodomain_solver.c index dcdeeefb..01a92307 100644 --- a/src/monodomain/monodomain_solver.c +++ b/src/monodomain/monodomain_solver.c @@ -668,45 +668,73 @@ int solve_monodomain(struct monodomain_solver *the_monodomain_solver, struct ode } } - if(purkinje_config) { + // REACTION TERM + // REACTION: Purkinje + if (purkinje_config) { start_stop_watch(&stop_watch); - // REACTION: Purkinje solve_all_volumes_odes(the_purkinje_ode_solver, cur_time, purkinje_stimuli_configs, configs->purkinje_ode_extra_config); - purkinje_ode_total_time += stop_stop_watch(&stop_watch); - - start_stop_watch(&stop_watch); - - // UPDATE: Purkinje ((update_monodomain_fn *)update_monodomain_config->main_function)(&time_info, update_monodomain_config, the_grid, the_monodomain_solver, the_grid->purkinje->num_active_purkinje_cells, the_grid->purkinje->purkinje_cells, the_purkinje_ode_solver, original_num_purkinje_cells); purkinje_ode_total_time += stop_stop_watch(&stop_watch); - + } + // REACTION: Myocardium + if (domain_config) { start_stop_watch(&stop_watch); - // TODO: show the purkinje fibers in the visualization tool - // #ifdef COMPILE_GUI - // if (show_gui) { - // omp_set_lock(&gui_config->draw_lock); - // } - // #endif + solve_all_volumes_odes(the_ode_solver, cur_time, stimuli_configs, configs->ode_extra_config); + + ((update_monodomain_fn *)update_monodomain_config->main_function)(&time_info, update_monodomain_config, the_grid, the_monodomain_solver, + the_grid->num_active_cells, the_grid->active_cells, the_ode_solver, + original_num_cells); + + ode_total_time += stop_stop_watch(&stop_watch); + + #ifdef COMPILE_GUI + if(show_gui) { + omp_set_lock(&gui_config->draw_lock); + } + #endif + } + + // PURKINJE-MUSCLE-JUNCTION COUPLING: + // TODO: Compute the PMJ coupling time ... + if(purkinje_config && domain_config) { + // Calculate the PMJ current from the Myocardium to the Purkinje (retrograde direction) + if (calc_retropropagation) { + compute_pmj_current_myocardium_to_purkinje(the_purkinje_ode_solver, the_ode_solver, the_grid, the_terminals); + } + // Calculate the PMJ current from the Purkinje to the Myocardium (anterograde direction) + compute_pmj_current_purkinje_to_myocardium(the_purkinje_ode_solver, the_ode_solver, the_grid, the_terminals); + } - // COUPLING: Calculate the PMJ current from the Tissue to the Purkinje - if(domain_config && calc_retropropagation) - compute_pmj_current_tissue_to_purkinje(the_purkinje_ode_solver, the_grid, the_terminals); + // DIFFUSION TERM: + // DIFFUSION: Purkinje + if (purkinje_config) { - // DIFUSION: Purkinje - if(purkinje_linear_system_solver_config) // Purkinje-coupled + start_stop_watch(&stop_watch); + + // Purkinje-coupled + if(purkinje_linear_system_solver_config) { ((linear_system_solver_fn *)purkinje_linear_system_solver_config->main_function)( &time_info, purkinje_linear_system_solver_config, the_grid, the_grid->purkinje->num_active_purkinje_cells, the_grid->purkinje->purkinje_cells, &purkinje_solver_iterations, &purkinje_solver_error); - else // Only-Purkinje + if(isnan(solver_error)) { + log_error("\nSimulation stopped due to NaN on time %lf. This is probably a problem with the Purkinje cellular model solver.\n.", cur_time); + + return SIMULATION_FINISHED; + } + } + // Only-Purkinje + else { ((linear_system_solver_fn *)linear_system_solver_config->main_function)( &time_info, linear_system_solver_config, the_grid, the_grid->purkinje->num_active_purkinje_cells, the_grid->purkinje->purkinje_cells, &purkinje_solver_iterations, &purkinje_solver_error); + } + purkinje_cg_partial = stop_stop_watch(&stop_watch); @@ -714,42 +742,22 @@ int solve_monodomain(struct monodomain_solver *the_monodomain_solver, struct ode purkinje_total_cg_it += purkinje_solver_iterations; } - - if(domain_config) { - - start_stop_watch(&stop_watch); - - // REACTION - solve_all_volumes_odes(the_ode_solver, cur_time, stimuli_configs, configs->ode_extra_config); - ((update_monodomain_fn *)update_monodomain_config->main_function)(&time_info, update_monodomain_config, the_grid, the_monodomain_solver, - the_grid->num_active_cells, the_grid->active_cells, the_ode_solver, - original_num_cells); - - ode_total_time += stop_stop_watch(&stop_watch); + // DIFFUSION: Myocardium + if (domain_config) { start_stop_watch(&stop_watch); - -#ifdef COMPILE_GUI - if(show_gui) { - omp_set_lock(&gui_config->draw_lock); - } -#endif - - // COUPLING: Calculate the PMJ current from the Purkinje to the Tissue - if(purkinje_config) { - compute_pmj_current_purkinje_to_tissue(the_ode_solver, the_grid, the_terminals); - } - - // DIFUSION: Tissue + ((linear_system_solver_fn *)linear_system_solver_config->main_function)( &time_info, linear_system_solver_config, the_grid, the_grid->num_active_cells, the_grid->active_cells, &solver_iterations, &solver_error); if(isnan(solver_error)) { - log_error("\nSimulation stopped due to NaN on time %lf. This is probably a problem with the cellular model solver.\n.", cur_time); -#ifdef COMPILE_GUI + log_error("\nSimulation stopped due to NaN on time %lf. This is probably a problem with the myocardium cellular model solver.\n.", cur_time); + + #ifdef COMPILE_GUI if(show_gui) { omp_unset_lock(&gui_config->draw_lock); } -#endif + #endif + return SIMULATION_FINISHED; } @@ -1493,6 +1501,230 @@ void compute_pmj_current_tissue_to_purkinje(struct ode_solver *the_purkinje_ode_ } } +void compute_pmj_current_purkinje_to_myocardium(struct ode_solver *the_purkinje_ode_solver, struct ode_solver *the_myocardium_ode_solver,\ + struct grid *the_grid, struct terminal *the_terminals) { + assert(the_myocardium_ode_solver); + assert(the_purkinje_ode_solver); + assert(the_grid); + assert(the_grid->purkinje); + assert(the_terminals); + + // Myocardium solution + struct cell_node **ac_myo = the_grid->active_cells; + uint32_t num_active_myo = the_grid->num_active_cells; + real *sv_myo = the_myocardium_ode_solver->sv; + uint32_t num_odes_myo = the_myocardium_ode_solver->model_data.number_of_ode_equations; + + // Purkinje solution + struct cell_node **ac_purk = the_grid->purkinje->purkinje_cells; + uint32_t num_active_purk = the_grid->purkinje->num_active_purkinje_cells; + real *sv_purk = the_purkinje_ode_solver->sv; + uint32_t num_odes_purk = the_purkinje_ode_solver->model_data.number_of_ode_equations; + + // Purkinje coupling parameters + real rpmj = the_grid->purkinje->network->rpmj; + real pmj_scale = the_grid->purkinje->network->pmj_scale; + + real Gpmj = 1.0 / rpmj; + + // TODO: Consider other combinations for the device that is solving the ODE + // Example: Purkinje=CPU, Myocardium=GPU || Purkinje=GPU, Myocardium=CPU + if(the_myocardium_ode_solver->gpu && the_purkinje_ode_solver->gpu) { + #ifdef COMPILE_CUDA + + real *vms_myo; + real *vms_purk; + uint32_t max_number_of_cells_myo = the_myocardium_ode_solver->original_num_cells; + uint32_t max_number_of_cells_purk = the_purkinje_ode_solver->original_num_cells; + size_t mem_size_myo = max_number_of_cells_myo * sizeof(real); + size_t mem_size_purk = max_number_of_cells_purk * sizeof(real); + + vms_myo = (real*)malloc(mem_size_myo); + vms_purk = (real*)malloc(mem_size_purk); + + check_cuda_error(cudaMemcpy(vms_myo, sv_myo, mem_size_myo, cudaMemcpyDeviceToHost)); + check_cuda_error(cudaMemcpy(vms_purk, sv_purk, mem_size_purk, cudaMemcpyDeviceToHost)); + + if(the_grid->adaptive) { + printf("[-] ERROR! The PMJ coupling is not yet implemented when using space adaptivity.\n"); + exit(EXIT_FAILURE); + } + + uint32_t num_of_purkinje_terminals = the_grid->purkinje->network->number_of_terminals; + for(uint32_t i = 0; i < num_of_purkinje_terminals; i++) { + + // Compute the PMJ current + real Ipmj = 0.0; + uint32_t num_tissue_cells = arrlen(the_terminals[i].tissue_cells); + uint32_t purkinje_index = the_terminals[i].purkinje_cell->id; + rpmj = the_terminals[i].purkinje_cell->rpmj; + Gpmj = 1.0 / rpmj; + for(uint32_t j = 0; j < num_tissue_cells; j++) { + uint32_t tissue_index = the_terminals[i].tissue_cells[j]->sv_position; + Ipmj += (vms_myo[tissue_index] - vms_purk[purkinje_index]); + } + Ipmj *= (Gpmj / pmj_scale); + + // Add this current to the RHS from each tissue cell + for(uint32_t j = 0; j < num_tissue_cells; j++) { + uint32_t tissue_index = the_terminals[i].tissue_cells[j]->sv_position; + + ac_myo[tissue_index]->b -= Ipmj; + } + } + + free(vms_myo); + free(vms_purk); + + #endif + } else if (!the_myocardium_ode_solver->gpu && !the_purkinje_ode_solver->gpu) { + uint32_t num_of_purkinje_terminals = the_grid->purkinje->network->number_of_terminals; + for(uint32_t i = 0; i < num_of_purkinje_terminals; i++) { + + // Compute the PMJ current + real Ipmj = 0.0; + uint32_t num_tissue_cells = arrlen(the_terminals[i].tissue_cells); + uint32_t purkinje_index = the_terminals[i].purkinje_cell->id; + rpmj = the_terminals[i].purkinje_cell->rpmj; + Gpmj = 1.0 / rpmj; + for(uint32_t j = 0; j < num_tissue_cells; j++) { + uint32_t tissue_index = the_terminals[i].tissue_cells[j]->sv_position; + Ipmj += (sv_myo[tissue_index * num_odes_myo] - sv_purk[purkinje_index * num_odes_purk]); + } + Ipmj *= (Gpmj / pmj_scale); + + // Add this current to the RHS from each tissue cell + for(uint32_t j = 0; j < num_tissue_cells; j++) { + uint32_t tissue_index = the_terminals[i].tissue_cells[j]->sv_position; + + ac_myo[tissue_index]->b -= Ipmj; + } + } + } + else { + printf("[purkinje coupling] ERROR! You are using the cellular models in different devices!\n"); + if (!the_myocardium_ode_solver->gpu) + printf("[purkinje coupling] Myocardium cellular model = GPU\n"); + else + printf("[purkinje coupling] Myocardium cellular model = CPU\n"); + if (!the_purkinje_ode_solver->gpu) + printf("[purkinje coupling] Purkinje cellular model = GPU\n"); + else + printf("[purkinje coupling] Purkinje cellular model = CPU\n"); + printf("[purkinje coupling] Please change your configuration file, so that all cellular models are on the same device!\n"); + exit(EXIT_FAILURE); + } +} + +void compute_pmj_current_myocardium_to_purkinje(struct ode_solver *the_purkinje_ode_solver, struct ode_solver *the_myocardium_ode_solver,\ + struct grid *the_grid, struct terminal *the_terminals) { + assert(the_purkinje_ode_solver); + assert(the_myocardium_ode_solver); + assert(the_grid); + assert(the_grid->purkinje); + assert(the_terminals); + + // Myocardium solution + struct cell_node **ac = the_grid->active_cells; + uint32_t num_active_myo = the_grid->num_active_cells; + real *sv_myo = the_myocardium_ode_solver->sv; + uint32_t num_odes_myo = the_myocardium_ode_solver->model_data.number_of_ode_equations; + + // Purkinje solution + struct cell_node **ac_purk = the_grid->purkinje->purkinje_cells; + uint32_t num_active_purk = the_grid->purkinje->num_active_purkinje_cells; + real *sv_purk = the_purkinje_ode_solver->sv; + uint32_t num_odes_purk = the_purkinje_ode_solver->model_data.number_of_ode_equations; + + // Purkinje coupling parameters + real rpmj = the_grid->purkinje->network->rpmj; + real pmj_scale = the_grid->purkinje->network->pmj_scale; + real asymm_ratio = the_grid->purkinje->network->asymm_ratio; + + real Gpmj = 1.0 / rpmj; + + // TODO: Consider other combinations for the device that is solving the ODE + // Example: Purkinje=CPU, Myocardium=GPU || Purkinje=GPU, Myocardium=CPU + if(the_myocardium_ode_solver->gpu && the_purkinje_ode_solver->gpu) { + #ifdef COMPILE_CUDA + + real *vms_myo; + real *vms_purk; + uint32_t max_number_of_cells_myo = the_myocardium_ode_solver->original_num_cells; + uint32_t max_number_of_cells_purk = the_purkinje_ode_solver->original_num_cells; + size_t mem_size_myo = max_number_of_cells_myo * sizeof(real); + size_t mem_size_purk = max_number_of_cells_purk * sizeof(real); + + vms_myo = (real*)malloc(mem_size_myo); + vms_purk = (real*)malloc(mem_size_purk); + + check_cuda_error(cudaMemcpy(vms_myo, sv_myo, mem_size_myo, cudaMemcpyDeviceToHost)); + check_cuda_error(cudaMemcpy(vms_purk, sv_purk, mem_size_purk, cudaMemcpyDeviceToHost)); + + if(the_grid->adaptive) { + printf("[-] ERROR! The PMJ coupling is not yet implemented when using space adaptivity.\n"); + exit(EXIT_FAILURE); + } + + uint32_t num_of_purkinje_terminals = the_grid->purkinje->network->number_of_terminals; + for(uint32_t i = 0; i < num_of_purkinje_terminals; i++) { + + // Compute the PMJ current + real Ipmj = 0.0; + uint32_t num_tissue_cells = arrlen(the_terminals[i].tissue_cells); + uint32_t purkinje_index = the_terminals[i].purkinje_cell->id; + rpmj = the_terminals[i].purkinje_cell->rpmj; + Gpmj = 1.0 / rpmj; + for(uint32_t j = 0; j < num_tissue_cells; j++) { + uint32_t tissue_index = the_terminals[i].tissue_cells[j]->sv_position; + Ipmj += (vms_purk[purkinje_index] - vms_myo[tissue_index]); + } + // Asymmetry of conduction across the PMJ + Ipmj *= (Gpmj / (pmj_scale * asymm_ratio)); + + // Add this current to the RHS of the Purkinje cell + ac_purk[purkinje_index]->b -= Ipmj; + } + + free(vms_myo); + free(vms_purk); + + #endif + } else if (!the_myocardium_ode_solver->gpu && !the_purkinje_ode_solver->gpu) { + uint32_t num_of_purkinje_terminals = the_grid->purkinje->network->number_of_terminals; + for(uint32_t i = 0; i < num_of_purkinje_terminals; i++) { + + // Compute the PMJ current + real Ipmj = 0.0; + uint32_t num_tissue_cells = arrlen(the_terminals[i].tissue_cells); + uint32_t purkinje_index = the_terminals[i].purkinje_cell->id; + rpmj = the_terminals[i].purkinje_cell->rpmj; + Gpmj = 1.0 / rpmj; + for(uint32_t j = 0; j < num_tissue_cells; j++) { + uint32_t tissue_index = the_terminals[i].tissue_cells[j]->sv_position; + Ipmj += (sv_purk[purkinje_index * num_odes_purk] - sv_myo[tissue_index * num_odes_myo]); + } + Ipmj *= (Gpmj / (pmj_scale * asymm_ratio)); + + // Add this current to the RHS of the Purkinje cell + ac_purk[purkinje_index]->b -= Ipmj; + } + } + else { + printf("[purkinje coupling] ERROR! You are using the cellular models in different devices!\n"); + if (!the_myocardium_ode_solver->gpu) + printf("[purkinje coupling] Myocardium cellular model = GPU\n"); + else + printf("[purkinje coupling] Myocardium cellular model = CPU\n"); + if (!the_purkinje_ode_solver->gpu) + printf("[purkinje coupling] Purkinje cellular model = GPU\n"); + else + printf("[purkinje coupling] Purkinje cellular model = CPU\n"); + printf("[purkinje coupling] Please change your configuration file, so that all cellular models are on the same device!\n"); + exit(EXIT_FAILURE); + } +} + // TODO: Maybe write a library to the PMJ coupling ... void write_pmj_delay(struct grid *the_grid, struct config *config, struct terminal *the_terminals) { assert(the_grid); @@ -1593,13 +1825,14 @@ void write_pmj_delay(struct grid *the_grid, struct config *config, struct termin // Check if the number of activations from the tissue and Purkinje cell are equal if(n_activations_purkinje > n_activations_tissue) { - // log_error("[purkinje_coupling] ERROR! The number of activations of the tissue and Purkinje cells are different!\n"); - // log_error("[purkinje_coupling] Probably there was a block on the anterograde direction!\n"); - // log_error("[purkinje_coupling] Consider only the result from the second pulse! (retrograde direction)!\n"); - // fprintf(output_file,"ERROR! Probably there was a block on the anterograde direction!\n"); has_block = true; cur_pulse = 0; - // return; + log_error("[purkinje_coupling] ERROR! The number of activations of the tissue and Purkinje cells are different!\n"); + log_error("[purkinje_coupling] Probably there was a block on the anterograde direction!\n"); + log_error("[purkinje_coupling] Consider only the result from the second pulse! (retrograde direction)!\n"); + fprintf(output_file,"%d,%u,%g,%g,%g,%d,%d\n", k, term_id, purkinje_lat, mean_tissue_lat, 0.0, (int)is_terminal_active, (int)has_block); + fclose(output_file); + return; } mean_tissue_lat += activation_times_array_tissue[cur_pulse]; if(activation_times_array_tissue[cur_pulse] < min_tissue_lat) { diff --git a/src/monodomain/monodomain_solver.h b/src/monodomain/monodomain_solver.h index 1331fe3a..cef5a33a 100644 --- a/src/monodomain/monodomain_solver.h +++ b/src/monodomain/monodomain_solver.h @@ -56,6 +56,11 @@ bool update_ode_state_vector_and_check_for_activity(real_cpu vm_threshold, struc void compute_pmj_current_purkinje_to_tissue (struct ode_solver *the_ode_solver, struct grid *the_grid, struct terminal *the_terminals); void compute_pmj_current_tissue_to_purkinje (struct ode_solver *the_purkinje_ode_solver, struct grid *the_grid, struct terminal *the_terminals); +void compute_pmj_current_purkinje_to_myocardium(struct ode_solver *the_purkinje_ode_solver, struct ode_solver *the_myocardium_ode_solver,\ + struct grid *the_grid, struct terminal *the_terminals); +void compute_pmj_current_myocardium_to_purkinje(struct ode_solver *the_purkinje_ode_solver, struct ode_solver *the_myocardium_ode_solver,\ + struct grid *the_grid, struct terminal *the_terminals); + void write_pmj_delay (struct grid *the_grid, struct config *config, struct terminal *the_terminals); void write_terminals_info (struct grid *the_grid, struct config *config, struct terminal *the_terminals);