diff --git a/.gitignore b/.gitignore index 15d042c5..02ba9797 100644 --- a/.gitignore +++ b/.gitignore @@ -170,6 +170,8 @@ src/*/custom_functions.c src/*/custom_functions_BACKUP.c src/*/custom_functions.h +bsbash/find_functions_with_mpi.sh + .cache/ .ccls-cache/ diff --git a/example_configs/introduction_to_monoalg3d/EX01_plain_wave.ini b/example_configs/introduction_to_monoalg3d/EX01_plain_wave.ini new file mode 100644 index 00000000..94a5fafb --- /dev/null +++ b/example_configs/introduction_to_monoalg3d/EX01_plain_wave.ini @@ -0,0 +1,90 @@ +# Version: 12/01/2024 +# =============================================================================================== +# Author: Lucas Berg (@bergolho) +# Last update: 12/01/2024 +# Description: Plain wave simulation using a slab (5,5cm x 5,5cm) using a +# space discretization of 200um. +# Stimulus: +# - Two pulses with a Basic Cycle Length (BCL) equal to 1000ms +# Cellular model: +# - Ten & Tusscher 3 +# ECG: +# - Two electrodes positioned on each side of the slab. +# +# ______ x = electrodes +# | | +# x | | x +# |______| +# ----------------------------------------------------------------------------------------------- +# Execute:> ./bin/MonoAlg3D -c example_configs/intro_to_monoalg3d/EX01_plain_wave.ini +# Visualize:> ./bin/MonoAlg3D_visualizer ./outputs/EX01_IntroMonoAlg_plain_mesh_healthy_200um +# - The simulation can be open on Paraview as well! +# =============================================================================================== + +[main] +num_threads=6 +dt_pde=0.02 ; miliseconds +simulation_time=2000.0 ; miliseconds +abort_on_no_activity=false +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=200 +output_dir=./outputs/EX01_IntroMonoAlg_plain_mesh_healthy_200um +add_timestamp=false +binary=true +main_function=save_as_ensight +remove_older_simulation=true + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.00005336 ; mS/um { ~44 cm/s } +sigma_y=0.00005336 ; mS/um { ~44 cm/s } +sigma_z=0.00005336 ; mS/um { ~44 cm/s } +main_function=homogeneous_sigma_assembly_matrix + +[linear_system_solver] +tolerance=1e-16 +use_preconditioner=no +max_iterations=200 +use_gpu=true +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient +main_function=conjugate_gradient + +[domain] +name=Plain Mesh +num_layers=1 +start_dx=200.0 ; micrometers +start_dy=200.0 ; micrometers +start_dz=200.0 ; micrometers +side_length=55000.0 ; micrometers +main_function=initialize_grid_with_square_mesh + +[ode_solver] +dt=0.02 ; miliseconds +use_gpu=yes +gpu_id=0 +library_file=shared_libs/libten_tusscher_3_endo.so + +[stim_plain] +start = 0.0 +duration = 2.0 ; miliseconds +period = 1000.0 ; miliseconds +current = -38.0 +x_limit = 500.0 ; micrometers +main_function=stim_if_x_less_than + +[calc_ecg] +main_function=pseudo_bidomain +init_function=init_pseudo_bidomain +end_function=end_pseudo_bidomain +calc_rate=10 +lead1=-5000,27500,50 ; micrometers +lead2=60000,27500,50 ; micrometers +sigma_b=20 +use_gpu=true +filename=./outputs/EX01_IntroMonoAlg_plain_mesh_healthy_200um/ecg.txt diff --git a/example_configs/introduction_to_monoalg3d/EX02_S1S2_protocol.ini b/example_configs/introduction_to_monoalg3d/EX02_S1S2_protocol.ini new file mode 100644 index 00000000..e67c2dca --- /dev/null +++ b/example_configs/introduction_to_monoalg3d/EX02_S1S2_protocol.ini @@ -0,0 +1,103 @@ +# Version: 12/01/2024 +# =============================================================================================== +# Author: Lucas Berg (@bergolho) +# Last update: 12/01/2024 +# Description: Simulation using a slab (5,5cm x 5,5cm) considering the S1S2 protocol with a +# space discretization of 200um. +# Stimulus: +# - First stimulus (S1=0ms) in the left border of the slab. +# - Second stimulus (S2=360ms) in the square [0,27500]x[0,27500] of the slab. +# Cellular model: +# - Ten & Tusscher 3 +# ECG: +# - Two electrodes positioned on each side of the slab. +# +# ______ x = electrodes +# | | +# x | | x +# |______| +# ----------------------------------------------------------------------------------------------- +# Execute:> ./bin/MonoAlg3D -c example_configs/intro_to_monoalg3d/EX02_S1S2_protocol.ini +# Visualize:> ./bin/MonoAlg3D_visualizer ./outputs/EX02_IntroMonoAlg_plain_mesh_S1S2_protocol_200um +# - The simulation can be open on Paraview as well! +# =============================================================================================== + +[main] +num_threads=6 +dt_pde=0.02 +simulation_time=2000.0 +abort_on_no_activity=true +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=200 +output_dir=./outputs/EX02_IntroMonoAlg_plain_mesh_S1S2_protocol_200um +add_timestamp=false +binary=true +main_function=save_as_ensight +remove_older_simulation=true + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.00005336 +sigma_y=0.00005336 +sigma_z=0.00005336 +main_function=homogeneous_sigma_assembly_matrix + +[linear_system_solver] +tolerance=1e-16 +use_preconditioner=no +max_iterations=200 +use_gpu=true +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient +main_function=conjugate_gradient + +[domain] +name=Plain Mesh S1S2 Protocol +num_layers=1 +start_dx=200.0 +start_dy=200.0 +start_dz=200.0 +side_length=55000 +main_function=initialize_grid_with_square_mesh + +[ode_solver] +dt=0.02 +use_gpu=yes +gpu_id=0 +library_file=shared_libs/libten_tusscher_3_endo.so + +; First stimulus S1 +[stim_plain_s1] +start = 0.0 +duration = 2.0 +current = -38.0 +x_limit = 500.0 +main_function=stim_if_x_less_than + +; Second stimulus S2 +[stim_plain_s2] +start = 360.0 +duration = 2.0 +current = -38.0 +min_x = 0.0 +min_y = 0.0 +max_x = 27550.0 +max_y = 27550.0 +main_function=stim_x_y_limits + +[calc_ecg] +main_function=pseudo_bidomain +init_function=init_pseudo_bidomain +end_function=end_pseudo_bidomain +calc_rate=10 +lead1=-5000,27500,50 +lead2=60000,27500,50 +sigma_b=20 +use_gpu=true +filename=./outputs/EX02_IntroMonoAlg_plain_mesh_S1S2_protocol_200um/ecg.txt + diff --git a/example_configs/introduction_to_monoalg3d/EX03_ischemia_modeling.ini b/example_configs/introduction_to_monoalg3d/EX03_ischemia_modeling.ini new file mode 100644 index 00000000..564d8645 --- /dev/null +++ b/example_configs/introduction_to_monoalg3d/EX03_ischemia_modeling.ini @@ -0,0 +1,123 @@ +# Version: 12/01/2024 +# =============================================================================================== +# Author: Lucas Berg (@bergolho) +# Last update: 12/01/2024 +# Description: Simulation using a slab (5,5cm x 5,5cm) considering ischemia remodeling around a +# circle centered at (2.75, 2.75) with a radius equal to 1.5cm and a border zone +# size of 0.75cm. Space discretization equal to 200um. +# Stimulus: +# - First stimulus (S1=0ms) on the left border of the slab. +# - Second stimulus (S2=375ms) on the left border of the slab. +# Cellular model: +# - Ten & Tusscher 3 +# Ischemia: +# - Change the action potential behaviour of the TT3 model inside the ischemic zone using the +# parameters: +# - Hyperkelemia = [Ko] +# - Hypoxia = [Ikatp_mod] +# - Acidosis = [INa_mod, ICaL_mod] +# ECG: +# - Two electrodes positioned on each side of the slab. +# +# ______ x = electrodes +# | | +# x | | x +# |______| +# ----------------------------------------------------------------------------------------------- +# Execute:> ./bin/MonoAlg3D -c example_configs/intro_to_monoalg3d/EX03_IntroMonoAlg_plain_mesh_ischemia_modeling_200um +# Visualize:> ./bin/MonoAlg3D_visualizer ./outputs/EX03_IntroMonoAlg_plain_mesh_ischemia_modeling_200um +# - The simulation can be open on Paraview as well! +# =============================================================================================== + +[main] +num_threads=6 +dt_pde=0.02 +simulation_time=2000.0 +abort_on_no_activity=true +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=100 +output_dir=./outputs/EX03_IntroMonoAlg_plain_mesh_ischemia_modeling_200um +add_timestamp=false +binary=true +main_function=save_as_ensight +remove_older_simulation=true + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.00005336 +sigma_y=0.00005336 +sigma_z=0.00005336 +main_function=homogeneous_sigma_assembly_matrix + +[linear_system_solver] +tolerance=1e-16 +use_preconditioner=no +max_iterations=200 +use_gpu=true +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient +main_function=conjugate_gradient + +[domain] +name=Plain Mesh with fibrosis and circle +start_dx=200.0 +start_dy=200.0 +start_dz=200.0 +num_layers=1 +side_length=55000 +seed=1508201274 +phi=0.0 +plain_center=27550.0 +sphere_radius=15000.0 +border_zone_radius=22500.0 +border_zone_size=7500.0 +main_function=initialize_grid_with_plain_and_sphere_fibrotic_mesh + +[ode_solver] +dt=0.02 +use_gpu=yes +gpu_id=0 +library_file=shared_libs/libten_tusscher_3_endo.so + +[stim_plain_s1] +start = 0.0 +duration = 2.0 +current = -38.0 +x_limit = 2500.0 +main_function=stim_if_x_less_than + +[stim_plain_s2] +start = 375.0 +duration = 2.0 +current = -38.0 +x_limit = 2500.0 +main_function=stim_if_x_less_than + +; Ischemia parameters +; Only the cells inside the ischemic region will change! +[extra_data] +Ko=9.0 ; Hyperkelemia +GNa_multiplicator=0.75 ; Ischemia +GCaL_multiplicator=0.75 ; Ischemia +Ikatp_multiplicator=0.01 ; Hypoxia +plain_center=27550.0 +sphere_radius=15000.0 +border_zone_radius=22500.0 +border_zone_size=7500.0 +main_function=set_extra_data_for_fibrosis_sphere + +[calc_ecg] +main_function=pseudo_bidomain +init_function=init_pseudo_bidomain +end_function=end_pseudo_bidomain +calc_rate=10 +lead1=-5000,27500,50 +lead2=60000,27500,50 +sigma_b=20 +use_gpu=true +filename=./outputs/EX03_IntroMonoAlg_plain_mesh_ischemia_modeling_200um/ecg.txt diff --git a/example_configs/introduction_to_monoalg3d/EX04_fibrosis_modeling.ini b/example_configs/introduction_to_monoalg3d/EX04_fibrosis_modeling.ini new file mode 100644 index 00000000..244edcb4 --- /dev/null +++ b/example_configs/introduction_to_monoalg3d/EX04_fibrosis_modeling.ini @@ -0,0 +1,107 @@ +# Version: 12/01/2024 +# =============================================================================================== +# Author: Lucas Berg (@bergolho) +# Last update: 12/01/2024 +# Description: Simulation using a slab (4cm x 4cm) ischemia remodeling and fibrosis around a +# circle centered at (2cm, 2cm) with a radius equal to 1.4cm and a border zone +# size of 0.2cm. Space discretization equal to 100um. +# Stimulus: +# - First stimulus (S1=0ms) on the left border of the slab. +# Cellular model: +# - Ten & Tusscher 3 +# ECG: +# - Two electrodes positioned on each side of the slab. +# +# ______ x = electrodes +# | | +# x | | x +# |______| +# Fibrosis parameters taken from this paper: +# - Sachetto, Rafael, Sergio Alonso, and Rodrigo Weber Dos Santos. +# "Killing many birds with two stones: hypoxia and fibrosis can generate ectopic beats in a human ventricular model." +# Frontiers in Physiology 9 (2018): 764. +# ----------------------------------------------------------------------------------------------- +# Execute:> ./bin/MonoAlg3D -c example_configs/intro_to_monoalg3d/EX04_fibrosis_modeling.ini +# Visualize:> ./bin/MonoAlg3D_visualizer ./outputs/EX04_IntroMonoAlg_plain_mesh_fibrosis_modeling_100um +# - The simulation can be open on Paraview as well! +# =============================================================================================== +[main] +num_threads=6 +dt_pde=0.02 +simulation_time=2000.0 +abort_on_no_activity=true +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=200 +output_dir=./outputs/EX04_IntroMonoAlg_plain_mesh_fibrosis_modeling_100um +add_timestamp=false +binary=true +main_function=save_as_ensight +remove_older_simulation=true + +[assembly_matrix] +sigma_x=0.00005336 +sigma_y=0.00005336 +sigma_z=0.00005336 +main_function=homogeneous_sigma_assembly_matrix +init_function=set_initial_conditions_fvm + +[linear_system_solver] +tolerance=1e-16 +use_preconditioner=no +max_iterations=200 +use_gpu=true +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient +main_function=conjugate_gradient + +[domain] +name=Plain Mesh with fibrosis and circle +start_dx=100.0 +start_dy=100.0 +start_dz=100.0 +main_function=initialize_grid_with_plain_and_sphere_fibrotic_mesh +num_layers=1 +side_length=40000 +phi=0.39 +plain_center=20050.0 +sphere_radius=14000.0 +border_zone_radius=16000.0 +border_zone_size=2000.0 +seed=1562001492 + +[ode_solver] +dt=0.02 +use_gpu=yes +gpu_id=0 +library_file=shared_libs/libten_tusscher_3_endo.so + +[stim_plain] +start = 0.0 +duration = 2.0 +current = -38.0 +x_limit = 500.0 +main_function=stim_if_x_less_than + +[extra_data] +atpi=2.0 +plain_center=20050.0 +sphere_radius=14000.0 +border_zone_radius=16000.0 +border_zone_size=2000.0 +main_function=set_extra_data_for_fibrosis_sphere + +[calc_ecg] +main_function=pseudo_bidomain +init_function=init_pseudo_bidomain +end_function=end_pseudo_bidomain +calc_rate=10 +lead1=-5000,20000,50 +lead2=45000,20000,50 +sigma_b=20 +use_gpu=true +filename=./outputs/EX04_IntroMonoAlg_plain_mesh_fibrosis_modeling_100um/ecg.txt diff --git a/example_configs/introduction_to_monoalg3d/EX05_transmurality_modeling_torord_land.ini b/example_configs/introduction_to_monoalg3d/EX05_transmurality_modeling_torord_land.ini new file mode 100644 index 00000000..260ed183 --- /dev/null +++ b/example_configs/introduction_to_monoalg3d/EX05_transmurality_modeling_torord_land.ini @@ -0,0 +1,96 @@ +# Version: 12/01/2024 +# =============================================================================================== +# Author: Lucas Berg (@bergolho) +# Last update: 12/01/2024 +# Description: Plain wave simulation using a slab (5,5cm x 5,5cm) using a +# space discretization of 200um. Here, we model transmurality across +# the domain using three types of cells (ENDO/MCELL/EPI). +# Stimulus: +# - Two pulses with a Basic Cycle Length (BCL) equal to 1000ms +# Cellular model: +# - ToRORd_Land +# ECG: +# - Two electrodes positioned on each side of the slab. +# +# ______ x = electrodes +# | | +# x | | x +# |______| +# ----------------------------------------------------------------------------------------------- +# Execute:> ./bin/MonoAlg3D -c example_configs/intro_to_monoalg3d/EX05_transmurality_modeling_torord_land.ini +# Visualize:> ./bin/MonoAlg3D_visualizer ./outputs/EX05_IntroMonoAlg_plain_mesh_healthy_transmurality_200um +# - The simulation can be open on Paraview as well! +# =============================================================================================== + +[main] +num_threads=6 +dt_pde=0.01 +simulation_time=2000.0 +abort_on_no_activity=false +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=200 +output_dir=./outputs/EX05_IntroMonoAlg_plain_mesh_healthy_transmurality_200um +add_timestamp=false +binary=true +main_function=save_as_ensight +remove_older_simulation=true + +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.00005336 +sigma_y=0.00005336 +sigma_z=0.00005336 +main_function=homogeneous_sigma_assembly_matrix + +[linear_system_solver] +tolerance=1e-16 +use_preconditioner=no +max_iterations=200 +use_gpu=true +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient +main_function=conjugate_gradient + +[domain] +name=Plain Mesh +num_layers=1 +start_dx=200.0 +start_dy=200.0 +start_dz=200.0 +side_length=55000.0 +main_function=initialize_grid_with_square_mesh + +[ode_solver] +adaptive=false +dt=0.01 +use_gpu=yes +gpu_id=0 +library_file= shared_libs/libToRORd_Land_mixed_endo_mid_epi.so + +[stim_plain] +start = 0.0 +duration = 1.0 +period = 1000.0 +current = -53.0 +x_limit = 500.0 +main_function=stim_if_x_less_than + +[extra_data] +main_function=set_extra_data_mixed_torord_Land_epi_mid_endo + +[calc_ecg] +main_function=pseudo_bidomain +init_function=init_pseudo_bidomain +end_function=end_pseudo_bidomain +calc_rate=10 +lead1=-5000,27500,50 +lead2=60000,27500,50 +sigma_b=20 +use_gpu=true +filename=./outputs/EX05_IntroMonoAlg_plain_mesh_healthy_transmurality_200um/ecg.txt + diff --git a/example_configs/introduction_to_monoalg3d/EX06_transmurality_modeling_tt3.ini b/example_configs/introduction_to_monoalg3d/EX06_transmurality_modeling_tt3.ini new file mode 100644 index 00000000..86f49f8c --- /dev/null +++ b/example_configs/introduction_to_monoalg3d/EX06_transmurality_modeling_tt3.ini @@ -0,0 +1,94 @@ +# Version: 12/01/2024 +# =============================================================================================== +# Author: Lucas Berg (@bergolho) +# Last update: 12/01/2024 +# Description: Plain wave simulation using a slab (5,5cm x 5,5cm) using a +# space discretization of 200um. Here, we model transmurality across +# the domain using three types of cells (ENDO/MCELL/EPI). +# Stimulus: +# - Two pulses with a Basic Cycle Length (BCL) equal to 1000ms +# Cellular model: +# - Ten & Tusscher 3 +# ECG: +# - Two electrodes positioned on each side of the slab. +# +# ______ x = electrodes +# | | +# x | | x +# |______| +# ----------------------------------------------------------------------------------------------- +# Execute:> ./bin/MonoAlg3D -c example_configs/intro_to_monoalg3d/EX06_transmurality_modeling_tt3.ini +# Visualize:> ./bin/MonoAlg3D_visualizer ./outputs/EX06_transmurality_modeling_tt3 +# - The simulation can be open on Paraview as well! +# =============================================================================================== +[main] +num_threads=6 +dt_pde=0.02 +simulation_time=500.0 +abort_on_no_activity=false +use_adaptivity=false + +[update_monodomain] +main_function=update_monodomain_default + +[save_result] +print_rate=200 +output_dir=./outputs/EX06_transmurality_modeling_tt3 +add_timestamp=false +binary=true +save_ode_state_variables=false +main_function=save_as_ensight + +; Anisotropic +[assembly_matrix] +init_function=set_initial_conditions_fvm +sigma_x=0.0001334 +sigma_y=0.0000176 +sigma_z=0.0000176 +main_function=homogeneous_sigma_assembly_matrix + +[linear_system_solver] +tolerance=1e-16 +use_preconditioner=no +max_iterations=500 +library_file=shared_libs/libdefault_linear_system_solver.so +use_gpu=yes +main_function=conjugate_gradient +init_function=init_conjugate_gradient +end_function=end_conjugate_gradient + +[domain] +name=Plain Mesh Transmurality +num_layers=1 +start_dx=200.0 +start_dy=200.0 +start_dz=200.0 +side_length=55000 +main_function=initialize_grid_with_square_mesh + +[ode_solver] +dt=0.02 +use_gpu=yes +gpu_id=0 +library_file=shared_libs/libten_tusscher_tt3_mixed_endo_mid_epi.so + +[extra_data] +main_function=set_extra_data_mixed_tt3 + +[stim_plain_s1] +start = 0.0 +duration = 2.0 +current = -38.0 +x_limit = 500.0 +main_function=stim_if_x_less_than + +[calc_ecg] +main_function=pseudo_bidomain +init_function=init_pseudo_bidomain +end_function=end_pseudo_bidomain +calc_rate=10 +lead1=-5000,27500,50 +lead2=60000,27500,50 +sigma_b=20 +use_gpu=true +filename=./outputs/EX06_transmurality_modeling_tt3/ecg.txt diff --git a/example_configs/introduction_to_monoalg3d/plot_ecg.py b/example_configs/introduction_to_monoalg3d/plot_ecg.py new file mode 100644 index 00000000..972db683 --- /dev/null +++ b/example_configs/introduction_to_monoalg3d/plot_ecg.py @@ -0,0 +1,61 @@ +''' +==================================================================================================================================== +Version: 12/01/2024 +Author: Lucas Berg (@bergolho) +==================================================================================================================================== +Script to read and plot the ECG data coming from a MonoAlg3D simulation using Matplotlib. +------------------------------------------------------------------------------------------------------------------------------------ +How to use:> python plot_ecg.py /ecg.txt +==================================================================================================================================== +''' + +import sys +import numpy as np +import matplotlib.pyplot as plt + +def read_ecg_readings (input_file): + data = np.genfromtxt(input_file) + nlin, ncol = np.shape(data) + timesteps = data[:,0] + currents = data[:,1:] + num_leads = ncol-1 + + return timesteps, currents, num_leads + +def plot_ecg_readings (t, currents, nleads): + + # Plot 3 signals => 1 = Lead A, 2 = Lead B, 3 = (Lead B + Lead A) + fig, axs = plt.subplots(nleads+1, 1, sharex=True, sharey=False) + fig.text(0.5, 0.04, 'Time (ms)', ha='center', fontsize=10) + fig.text(0.02, 0.5, 'Current (mA)', va='center', rotation='vertical', fontsize=10) + # Signals 1 and 2 + for i in range(nleads): + axs[i].plot(t, currents[:,i], label="lead_%d" % (i), c="black", linewidth=2.0) + axs[i].set_title("ECG reading - Lead %d" % (i),fontsize=10) + # Signal 3 + axs[nleads].plot(t, currents[:,1]+currents[:,0], label="lead_%d" % (nleads), c="black", linewidth=2.0) + axs[nleads].set_title("ECG reading - Lead %d" % (nleads),fontsize=10) + + #plt.show() + #plt.savefig("ecg.pdf") + plt.savefig("ecg.png", dpi=300) + print("[+] ECG readings saved in the current folder with the name 'ecg.png'!") + +def main(): + + if len(sys.argv) != 2: + print("-------------------------------------------------------------------------") + print("Usage:> python %s " % sys.argv[0]) + print("-------------------------------------------------------------------------") + print(" = Input file with the ECG reading from each timestep") + print("-------------------------------------------------------------------------") + return 1 + + input_file = sys.argv[1] + + t, currents, num_leads = read_ecg_readings(input_file) + + plot_ecg_readings(t,currents,num_leads) + +if __name__ == "__main__": + main() diff --git a/src/common_types/common_types.h b/src/common_types/common_types.h index caa9c61e..09f6049f 100644 --- a/src/common_types/common_types.h +++ b/src/common_types/common_types.h @@ -67,6 +67,13 @@ struct fiber_coords { real_cpu n[3]; }; +struct fiber_coords_scale { + real_cpu f[3]; + real_cpu s[3]; + real_cpu n[3]; + real_cpu x[3]; +}; + struct condutivity { real_cpu x, y, z, xy, xz, yz; struct fiber_coords fibers; diff --git a/src/domains_library/domain.c b/src/domains_library/domain.c index e9fc51a6..cd92da15 100644 --- a/src/domains_library/domain.c +++ b/src/domains_library/domain.c @@ -361,3 +361,45 @@ SET_SPATIAL_DOMAIN(initialize_grid_with_custom_mesh) { return ret; } + +SET_SPATIAL_DOMAIN(initialize_grid_with_square_mesh_and_source_sink_fibrotic_region) { + + real_cpu phi = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, phi, config, "phi"); + + unsigned seed = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(unsigned, seed, config, "seed"); + + real_cpu min_x = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, min_x, config, "region_min_x"); + + real_cpu max_x = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, max_x, config, "region_max_x"); + + real_cpu min_y = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, min_y, config, "region_min_y"); + + real_cpu max_y = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, max_y, config, "region_max_y"); + + real_cpu min_z = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, min_z, config, "region_min_z"); + + real_cpu max_z = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, max_z, config, "region_max_z"); + + real_cpu source_sink_min_x = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, source_sink_min_x, config, "source_sink_min_x"); + + real_cpu source_sink_max_x = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, source_sink_max_x, config, "source_sink_max_x"); + + real_cpu side_length = 0.0; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, side_length, config, "side_length"); + + set_square_mesh(config, the_grid); + + set_plain_fibrosis_source_sink_region(the_grid, phi, seed, min_x, max_x, min_y, max_y, min_z, max_z, source_sink_min_x, source_sink_max_x, side_length); + + 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 bfbb8f67..6e825f68 100644 --- a/src/domains_library/domain_helpers.c +++ b/src/domains_library/domain_helpers.c @@ -1009,3 +1009,45 @@ int calc_num_refs(real_cpu start_h, real_cpu desired_h) { return num_refs; } + +void set_plain_fibrosis_source_sink_region (struct grid *the_grid, real_cpu phi, unsigned fib_seed, const double min_x, const double max_x, const double min_y, + const double max_y, const double min_z, const double max_z, + real_cpu source_sink_min_x, real_cpu source_sink_max_x, real_cpu side_length) { + log_info("Making %.2lf %% of cells inside the region inactive\n", phi * 100.0); + + struct cell_node *grid_cell; + + if(fib_seed == 0) + fib_seed = (unsigned)time(NULL) + getpid(); + + srand(fib_seed); + + log_info("Using %u as seed\n", fib_seed); + + real_cpu a1 = (2.0*side_length) / (side_length - 2*source_sink_min_x); + real_cpu b1 = -source_sink_min_x*a1; + real_cpu a2 = (2.0*side_length) / (side_length - 2*source_sink_max_x); + real_cpu b2 = -source_sink_max_x*a2; + + grid_cell = the_grid->first_cell; + while(grid_cell != 0) { + real center_x = grid_cell->center.x; + real center_y = grid_cell->center.y; + real center_z = grid_cell->center.z; + + if(center_x >= min_x && center_x <= max_x && center_y >= min_y && center_y <= max_y && center_z >= min_z && center_z <= max_z + && (center_y > a1*center_x + b1 || center_y > a2*center_x + b2)) { + if(grid_cell->active) { + real_cpu p = (real_cpu)(rand()) / (RAND_MAX); + if(p < phi) { + grid_cell->active = false; + } + + INITIALIZE_FIBROTIC_INFO(grid_cell); + FIBROTIC(grid_cell) = true; + } + } + + grid_cell = grid_cell->next; + } +} \ No newline at end of file diff --git a/src/domains_library/domain_helpers.h b/src/domains_library/domain_helpers.h index 98a7d281..9eb3d203 100644 --- a/src/domains_library/domain_helpers.h +++ b/src/domains_library/domain_helpers.h @@ -46,6 +46,10 @@ void set_plain_fibrosis_using_file(struct grid *the_grid, const char filename[]) void set_plain_fibrosis_inside_region(struct grid *the_grid, real_cpu phi, unsigned fib_seed, double min_x, double max_x, double min_y, double max_y, double min_z, double max_z); +void set_plain_fibrosis_source_sink_region (struct grid *the_grid, real_cpu phi, unsigned fib_seed, const double min_x, const double max_x, const double min_y, + const double max_y, const double min_z, const double max_z, + real_cpu source_sink_min_x, real_cpu source_sink_max_x, real_cpu side_length); + uint32_t set_custom_mesh_from_file(struct grid *the_grid, const char *mesh_file, uint32_t num_volumes, double start_h, uint8_t num_extra_fields, set_custom_data_for_mesh_fn set_custom_data_for_mesh); diff --git a/src/extra_data_library/extra_data.c b/src/extra_data_library/extra_data.c index 1dc64b8e..9d8a7e7d 100644 --- a/src/extra_data_library/extra_data.c +++ b/src/extra_data_library/extra_data.c @@ -429,6 +429,48 @@ SET_EXTRA_DATA (set_mixed_model_purkinje_and_tissue) return (void*)mapping; } +// 'libten_tusscher_tt3_mixed_endo_mid_epi.so' with transmurality and fibrosis (all cells healthy) +SET_EXTRA_DATA (set_extra_data_mixed_tt3) { + uint32_t num_active_cells = the_grid->num_active_cells; + real side_length = the_grid->mesh_side_length.x; + struct cell_node ** ac = the_grid->active_cells; + + // + struct extra_data_for_tt3 *extra_data = NULL; + extra_data = set_common_tt3_data(config, num_active_cells); + + // Divide the domain in three sections (ENDO/MID/EPI) + // The percentages were taken from the ToRORd paper (Transmural experiment) + real side_length_endo = side_length*0.45; + real side_length_mid = side_length_endo + side_length*0.25; + real side_length_epi = side_length_mid + side_length*0.3; + + int i; + + // Transmurality and fibrosis tags + OMP(parallel for) + for (int i = 0; i < num_active_cells; i++) { + + real center_x = ac[i]->center.x; + + // Tag the model transmurality + // ENDO=0, MID=1, EPI=2 + if (center_x < side_length_endo) + extra_data->transmurality[i] = 0.0; + else if (center_x >= side_length_endo && center_x < side_length_mid) + extra_data->transmurality[i] = 1.0; + else + extra_data->transmurality[i] = 2.0; + + // Tag the fibrosis region + extra_data->fibrosis[i] = 1.0; + } + + SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_tt3)); + + return (void*)extra_data; +} + // Initial condition - 'libToRORd_fkatp_mixed_endo_mid_epi.so' + transmurality + current modifiers (plain and cuboid) SET_EXTRA_DATA(set_extra_data_mixed_torord_fkatp_epi_mid_endo) { diff --git a/src/extra_data_library/helper_functions.c b/src/extra_data_library/helper_functions.c index 8d709992..729490fb 100644 --- a/src/extra_data_library/helper_functions.c +++ b/src/extra_data_library/helper_functions.c @@ -32,14 +32,60 @@ struct extra_data_for_fibrosis* set_common_schemia_data(struct config *config, u real Vm_modifier = 0.0f; GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Vm_modifier, config, "Vm_modifier"); + real Ikatp_multiplicator = 1.0f; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Ikatp_multiplicator, config, "Ikatp_multiplicator"); + + extra_data->atpi = atpi; + extra_data->Ko = Ko; + extra_data->Ki = Ki; + extra_data->GNa_multiplicator = GNa_multiplicator; + extra_data->GCaL_multiplicator = GCaL_multiplicator; + extra_data->INaCa_multiplicator = INaCa_multiplicator; + extra_data->Ikatp_multiplicator = Ikatp_multiplicator; + extra_data->Vm_modifier = Vm_modifier; + extra_data->fibrosis = MALLOC_ARRAY_OF_TYPE(real, num_cells); + + return extra_data; +} + +struct extra_data_for_tt3* set_common_tt3_data (struct config *config, uint32_t num_cells) { + + struct extra_data_for_tt3 *extra_data = MALLOC_ONE_TYPE(struct extra_data_for_tt3); + + real atpi = 6.8; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, atpi, config, "atpi"); + + real Ko = 5.4; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Ko, config, "Ko"); + + real Ki = 138.3; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Ki, config, "Ki"); + + real GNa_multiplicator = 1.0f; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, GNa_multiplicator, config, "GNa_multiplicator"); + + real GCaL_multiplicator = 1.0f; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, GCaL_multiplicator, config, "GCaL_multiplicator"); + + real INaCa_multiplicator = 1.0f; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, INaCa_multiplicator, config, "INaCa_multiplicator"); + + real Vm_modifier = 0.0f; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Vm_modifier, config, "Vm_modifier"); + + real Ikatp_multiplicator = 1.0f; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, Ikatp_multiplicator, config, "Ikatp_multiplicator"); + extra_data->atpi = atpi; extra_data->Ko = Ko; extra_data->Ki = Ki; extra_data->GNa_multiplicator = GNa_multiplicator; extra_data->GCaL_multiplicator = GCaL_multiplicator; extra_data->INaCa_multiplicator = INaCa_multiplicator; + extra_data->Ikatp_multiplicator = Ikatp_multiplicator; extra_data->Vm_modifier = Vm_modifier; extra_data->fibrosis = MALLOC_ARRAY_OF_TYPE(real, num_cells); + extra_data->transmurality = MALLOC_ARRAY_OF_TYPE(real, num_cells); return extra_data; } diff --git a/src/extra_data_library/helper_functions.h b/src/extra_data_library/helper_functions.h index b15fdadb..ab812114 100644 --- a/src/extra_data_library/helper_functions.h +++ b/src/extra_data_library/helper_functions.h @@ -18,10 +18,24 @@ struct extra_data_for_fibrosis { real GNa_multiplicator; real GCaL_multiplicator; real INaCa_multiplicator; + real Ikatp_multiplicator; real Vm_modifier; real *fibrosis; }; +struct extra_data_for_tt3 { + real atpi; + real Ko; + real Ki; + real GNa_multiplicator; + real GCaL_multiplicator; + real INaCa_multiplicator; + real Ikatp_multiplicator; + real Vm_modifier; + real *fibrosis; + real *transmurality; +}; + struct extra_data_for_torord { real INa_Multiplier; real ICaL_Multiplier; @@ -107,6 +121,7 @@ struct extra_data_for_trovato { }; 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); diff --git a/src/matrix_assembly_library/assembly_common.c b/src/matrix_assembly_library/assembly_common.c index 0178ebed..c844e05e 100644 --- a/src/matrix_assembly_library/assembly_common.c +++ b/src/matrix_assembly_library/assembly_common.c @@ -1185,3 +1185,43 @@ static struct fiber_coords *read_fibers(char *fiber_file_path, bool normalize_ve return fibers; } + +// Albert`s code +static struct fiber_coords_scale *read_fibers_scale(char *fiber_file_path) { + + FILE *fibers_file = open_file_or_exit(fiber_file_path, "r"); + + struct fiber_coords_scale *fibers = NULL; + char *line = NULL; + size_t len; + + while((getline(&line, &len, fibers_file)) != -1) { + + int split_count; + sds *points = sdssplit(line, " ", &split_count); + struct fiber_coords_scale f_coords; + + f_coords.f[0] = strtod(points[0], NULL); + f_coords.f[1] = strtod(points[1], NULL); + f_coords.f[2] = strtod(points[2], NULL); + + f_coords.s[0] = strtod(points[3], NULL); + f_coords.s[1] = strtod(points[4], NULL); + f_coords.s[2] = strtod(points[5], NULL); + + f_coords.n[0] = strtod(points[6], NULL); + f_coords.n[1] = strtod(points[7], NULL); + f_coords.n[2] = strtod(points[8], NULL); + + f_coords.x[0] = strtod(points[9], NULL); + f_coords.x[1] = strtod(points[10], NULL); + f_coords.x[2] = strtod(points[11], NULL); + + arrput(fibers, f_coords); + sdsfreesplitres(points, split_count); + } + + free(line); + + return fibers; +} diff --git a/src/models_library/ten_tusscher/build.sh b/src/models_library/ten_tusscher/build.sh index ad4888bd..867b96ab 100644 --- a/src/models_library/ten_tusscher/build.sh +++ b/src/models_library/ten_tusscher/build.sh @@ -53,3 +53,11 @@ COMMON_HEADERS="ten_tusscher_2004_mixed_endo_mid_epi.h" COMPILE_MODEL_LIB "ten_tusscher_2004_mixed_endo_mid_epi" "$MODEL_FILE_CPU" "$MODEL_FILE_GPU" "$COMMON_HEADERS" ############################################################################################ +############## TENTUSSCHER 3 MIXED ENDO MID EPI ##################################################### +MODEL_FILE_CPU="ten_tusscher_tt3_mixed_endo_mid_epi.c" +MODEL_FILE_GPU="ten_tusscher_tt3_mixed_endo_mid_epi.cu" +COMMON_HEADERS="ten_tusscher_tt3_mixed_endo_mid_epi.h" + +COMPILE_MODEL_LIB "ten_tusscher_tt3_mixed_endo_mid_epi" "$MODEL_FILE_CPU" "$MODEL_FILE_GPU" "$COMMON_HEADERS" +############################################################################################ + diff --git a/src/models_library/ten_tusscher/ten_tusscher_3_RS_CPU.c b/src/models_library/ten_tusscher/ten_tusscher_3_RS_CPU.c index 20b11708..61c1adad 100644 --- a/src/models_library/ten_tusscher/ten_tusscher_3_RS_CPU.c +++ b/src/models_library/ten_tusscher/ten_tusscher_3_RS_CPU.c @@ -66,7 +66,7 @@ SOLVE_MODEL_ODES(solve_model_odes_cpu) { real dt = ode_solver->min_dt; uint32_t num_steps = ode_solver->num_steps; - int num_extra_parameters = 7; + int num_extra_parameters = 8; real extra_par[num_extra_parameters]; real fibs_size = num_cells_to_solve*sizeof(real); @@ -82,6 +82,7 @@ SOLVE_MODEL_ODES(solve_model_odes_cpu) { extra_par[4] = extra_data_from_solver->GNa_multiplicator; extra_par[5] = extra_data_from_solver->GCaL_multiplicator; extra_par[6] = extra_data_from_solver->INaCa_multiplicator; + extra_par[7] = extra_data_from_solver->Ikatp_multiplicator; } else { // Default values for a healthy cell /////////// @@ -92,6 +93,7 @@ SOLVE_MODEL_ODES(solve_model_odes_cpu) { extra_par[4] = 1.0f; extra_par[5] = 1.0f; extra_par[6] = 1.0f; + extra_par[7] = 1.0f; fibrosis = (real*) malloc(fibs_size); diff --git a/src/models_library/ten_tusscher/ten_tusscher_3_RS_GPU.cu b/src/models_library/ten_tusscher/ten_tusscher_3_RS_GPU.cu index cdc108e0..3203964b 100644 --- a/src/models_library/ten_tusscher/ten_tusscher_3_RS_GPU.cu +++ b/src/models_library/ten_tusscher/ten_tusscher_3_RS_GPU.cu @@ -66,7 +66,7 @@ extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) { real *fibrosis_device; real *fibs = NULL; - int num_extra_parameters = 7; + int num_extra_parameters = 8; real extra_par[num_extra_parameters]; size_t extra_parameters_size = num_extra_parameters * sizeof(real); @@ -87,6 +87,7 @@ extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) { extra_par[4] = extra_data_from_cpu->GNa_multiplicator; extra_par[5] = extra_data_from_cpu->GCaL_multiplicator; extra_par[6] = extra_data_from_cpu->INaCa_multiplicator; + extra_par[7] = extra_data_from_cpu->Ikatp_multiplicator; } else { extra_par[0] = 6.8f; @@ -96,6 +97,7 @@ extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) { extra_par[4] = 1.0f; extra_par[5] = 1.0f; extra_par[6] = 1.0f; + extra_par[7] = 1.0f; fibs = (real*) malloc(fibs_size); diff --git a/src/models_library/ten_tusscher/ten_tusscher_3_RS_common.inc b/src/models_library/ten_tusscher/ten_tusscher_3_RS_common.inc index c8abfe46..53c5550c 100644 --- a/src/models_library/ten_tusscher/ten_tusscher_3_RS_common.inc +++ b/src/models_library/ten_tusscher/ten_tusscher_3_RS_common.inc @@ -5,7 +5,7 @@ //Linear changing of atpi depending on the fibrosis and distance from the center of the scar (only for border zone cells) real atpi = extra_parameters[0]; real atpi_change = 6.8f - atpi; - atpi = atpi +atpi_change*fibrosis; + atpi = atpi + atpi_change*fibrosis; //Extracellular potassium concentration was elevated //from its default value of 5.4 mM to values between 6.0 and 8.0 mM @@ -32,6 +32,10 @@ real INaCa_multplicator = extra_parameters[6]; real INaCa_multplicator_change = 1.0f - INaCa_multplicator; INaCa_multplicator = INaCa_multplicator + INaCa_multplicator_change*fibrosis; + + real Ikatp_multiplicator = extra_parameters[7]; + real Ikatp_multiplicator_change = 1.0f - Ikatp_multiplicator; + Ikatp_multiplicator = Ikatp_multiplicator + Ikatp_multiplicator_change*fibrosis; //real katp = 0.306; @@ -203,7 +207,7 @@ IbNa=GbNa*(svolt-Ena); IbCa=GbCa*(svolt-Eca); - IKatp = gkbaratp*(svolt-Ek); + IKatp = gkbaratp*(svolt-Ek) * Ikatp_multiplicator; //Determine total current diff --git a/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.c b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.c new file mode 100644 index 00000000..110cedc6 --- /dev/null +++ b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.c @@ -0,0 +1,165 @@ +#include +#include +#include "ten_tusscher_tt3_mixed_endo_mid_epi.h" + +GET_CELL_MODEL_DATA(init_cell_model_data) { + + assert(cell_model); + + if(get_initial_v) + cell_model->initial_v = INITIAL_V; + if(get_neq) + cell_model->number_of_ode_equations = NEQ; + +} + +SET_ODE_INITIAL_CONDITIONS_CPU(set_model_initial_conditions_cpu) { + + log_info("Using Mixed Ten & Tusscher 3 CPU model\n"); + + uint32_t num_cells = solver->original_num_cells; + + solver->sv = (real*)malloc(NEQ*num_cells*sizeof(real)); + + OMP(parallel for) + for(uint32_t i = 0; i < num_cells; i++) { + + real *sv = &solver->sv[i * NEQ]; + + sv[0] = -86.2f; // V; millivolt + sv[1] = 0.0f; //M + sv[2] = 0.75; //H + sv[3] = 0.75; //J + sv[4] = 0.0f; //Xr1 + sv[5] = 0.0f; //Xs + sv[6] = 1.0f; //S + sv[7] = 1.0f; //F + sv[8] = 1.0f; //F2 + sv[9] = 0.0; //D_INF + sv[10] = 0.0; //R_INF + sv[11] = 0.0; //Xr2_INF + } +} + +SOLVE_MODEL_ODES(solve_model_odes_cpu) { + + uint32_t sv_id; + real *transmurality = NULL; + real *fibrosis = NULL; + + size_t num_cells_to_solve = ode_solver->num_cells_to_solve; + uint32_t * cells_to_solve = ode_solver->cells_to_solve; + real *sv = ode_solver->sv; + real dt = ode_solver->min_dt; + uint32_t num_steps = ode_solver->num_steps; + + int num_extra_parameters = 8; + real extra_par[num_extra_parameters]; + + struct extra_data_for_tt3* extra_data_from_solver = (struct extra_data_for_tt3*)ode_solver->ode_extra_data; + bool deallocate = false; + + if(ode_solver->ode_extra_data) { + transmurality = extra_data_from_solver->transmurality; + fibrosis = extra_data_from_solver->fibrosis; + extra_par[0] = extra_data_from_solver->atpi; + extra_par[1] = extra_data_from_solver->Ko; + extra_par[2] = extra_data_from_solver->Ki; + extra_par[3] = extra_data_from_solver->Vm_modifier; + extra_par[4] = extra_data_from_solver->GNa_multiplicator; + extra_par[5] = extra_data_from_solver->GCaL_multiplicator; + extra_par[6] = extra_data_from_solver->INaCa_multiplicator; + extra_par[7] = extra_data_from_solver->Ikatp_multiplicator; + } + else { + // Default values for a healthy cell /////////// + extra_par[0] = 6.8f; + extra_par[1] = 5.4f; + extra_par[2] = 138.3f; + extra_par[3] = 0.0; + extra_par[4] = 1.0f; + extra_par[5] = 1.0f; + extra_par[6] = 1.0f; + extra_par[7] = 1.0f; + + fibrosis = (real*)malloc(num_cells_to_solve*sizeof(real)); + transmurality = (real*)malloc(num_cells_to_solve*sizeof(real)); + + // Default case: all cells ENDO and Healthy + for(uint64_t i = 0; i < num_cells_to_solve; i++) { + fibrosis[i] = 1.0; // Healthy + transmurality[i] = 0.0; // ENDO + } + + deallocate = true; + } + + int i; + + OMP(parallel for private(sv_id)) + for (i = 0; i < num_cells_to_solve; i++) { + if(cells_to_solve) + sv_id = cells_to_solve[i]; + else + sv_id = i; + + for (int j = 0; j < num_steps; ++j) { + solve_model_ode_cpu(dt, sv + (sv_id * NEQ), stim_currents[i], fibrosis[i], transmurality[i], extra_par); + } + } + + if(deallocate) { + free(fibrosis); + free(transmurality); + } +} + + +void solve_model_ode_cpu(real dt, real *sv, real stim_current, real fibrosis, real transmurality, real *extra_parameters) { + + assert(sv); + + real rY[NEQ], rDY[NEQ]; + + for(int i = 0; i < NEQ; i++) + rY[i] = sv[i]; + + RHS_cpu(rY, rDY, stim_current, dt, fibrosis, transmurality, extra_parameters); + + //THIS MODEL USES THE Rush Larsen Method TO SOLVE THE EDOS + sv[0] = dt*rDY[0] + rY[0]; + sv[1] = rDY[1]; + sv[2] = rDY[2]; + sv[3] = rDY[3]; + sv[4] = rDY[4]; + sv[5] = rDY[5]; + sv[6] = rDY[6]; + sv[7] = rDY[7]; + sv[8] = rDY[8]; + sv[9] = rDY[9]; + sv[10] = rDY[10]; + sv[11] = rDY[11]; +} + + +void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, real fibrosis, real transmurality, real const *extra_parameters) { + + //fibrosis = 0 means that the cell is fibrotic, 1 is not fibrotic. Anything between 0 and 1 means border zone + //transmurality = 0 means cell is endocardium, 1 is mid-myocardium and 2 is epicardium + + //THIS IS THE STATE VECTOR THAT WE NEED TO SAVE IN THE STEADY STATE + const real svolt = sv[0]; + const real sm = sv[1]; + const real sh = sv[2]; + const real sj = sv[3]; + const real sxr1 = sv[4]; + const real sxs = sv[5]; + const real ss = sv[6]; + const real sf = sv[7]; + const real sf2 = sv[8]; + const real D_INF = sv[9]; + const real R_INF = sv[10]; + const real Xr2_INF = sv[11]; + + #include "ten_tusscher_tt3_mixed_endo_mid_epi.common.c" +} diff --git a/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.common.c b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.common.c new file mode 100644 index 00000000..a5e1b569 --- /dev/null +++ b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.common.c @@ -0,0 +1,321 @@ + const real natp = 0.24; // K dependence of ATP-sensitive K current + const real nicholsarea = 0.00005; // Nichol's areas (cm^2) + const real hatp = 2; // Hill coefficient + + //Linear changing of atpi depending on the fibrosis and distance from the center of the scar (only for border zone cells) + real atpi = extra_parameters[0]; + real atpi_change = 6.8f - atpi; + atpi = atpi + atpi_change*fibrosis; + + //Extracellular potassium concentration was elevated + //from its default value of 5.4 mM to values between 6.0 and 8.0 mM + //Ref: A Comparison of Two Models of Human Ventricular Tissue: Simulated Ischemia and Re-entry + real Ko = extra_parameters[1]; + real Ko_change = 5.4f - Ko; + Ko = Ko + Ko_change*fibrosis; + + real Ki = extra_parameters[2]; + real Ki_change = 138.3 - Ki; + Ki = Ki + Ki_change*fibrosis; + + real Vm_modifier = extra_parameters[3]; + Vm_modifier = Vm_modifier - Vm_modifier*fibrosis; + + real GNa_multplicator = extra_parameters[4]; + real GNa_multplicator_change = 1.0f - GNa_multplicator; + GNa_multplicator = GNa_multplicator + GNa_multplicator_change*fibrosis; + + real GCaL_multplicator = extra_parameters[5]; + real GCaL_multplicator_change = 1.0f - GCaL_multplicator; + GCaL_multplicator = GCaL_multplicator + GCaL_multplicator_change*fibrosis; + + real INaCa_multplicator = extra_parameters[6]; + real INaCa_multplicator_change = 1.0f - INaCa_multplicator; + INaCa_multplicator = INaCa_multplicator + INaCa_multplicator_change*fibrosis; + + real Ikatp_multiplicator = extra_parameters[7]; + real Ikatp_multiplicator_change = 1.0f - Ikatp_multiplicator; + Ikatp_multiplicator = Ikatp_multiplicator + Ikatp_multiplicator_change*fibrosis; + + + //real katp = 0.306; + //Ref: A Comparison of Two Models of Human Ventricular Tissue: Simulated Ischaemia and Re-entry + //real katp = 0.306; + const real katp = -0.0942857142857*atpi + 0.683142857143; //Ref: A Comparison of Two Models of Human Ventricular Tissue: Simulated Ischaemia and Re-entry + + const real patp = 1/(1 + pow((atpi/katp),hatp)); + const real gkatp = 0.000195/nicholsarea; + const real gkbaratp = gkatp*patp*pow((Ko/5.4),natp); + + const real katp2= 1.4; + const real hatp2 = 2.6; + const real pcal = 1.0/(1.0 + pow((katp2/atpi),hatp2)); + + + const real Cao=2.0; + const real Nao=140.0; + const real Cai=0.00007; + const real Nai=7.67; + +//Constants + const real R=8314.472; + const real F=96485.3415; + const real T=310.0; + const real RTONF=(R*T)/F; + +//Parameters for currents +//Parameters for IKr + const real Gkr=0.101; +//Parameters for Iks + const real pKNa=0.03; + + real Gks = 0.257; + if (transmurality == EPI) { + Gks=0.257; + } + else if (transmurality == ENDO) { + Gks=0.392; + } + else if (transmurality == MID) { + Gks=0.098; + } +//Parameters for Ik1 + real GK1=5.405; +//Parameters for Ito + real Gto = 0.294; + if (transmurality == EPI) { + Gto=0.294; + } + else if (transmurality == ENDO) { + Gto=0.073; + } + else if (transmurality == MID) { + Gto=0.294; + } +//Parameters for INa + const real GNa=14.838*GNa_multplicator; //ACIDOSIS +//Parameters for IbNa + const real GbNa=0.00029; +//Parameters for INaK + const real KmK=1.0; + const real KmNa=40.0; + const real knak=2.724; +//Parameters for ICaL + const real GCaL=0.2786*pcal*GCaL_multplicator; //ACIDOSIS +//Parameters for IbCa + const real GbCa=0.000592; +//Parameters for INaCa + const real knaca=1000; + const real KmNai=87.5; + const real KmCa=1.38; + const real ksat=0.1; + const real n=0.35; +//Parameters for IpCa + const real GpCa=0.1238; + const real KpCa=0.0005; +//Parameters for IpK; + const real GpK=0.0293; + + const real Ek=RTONF*(log((Ko/Ki))); + const real Ena=RTONF*(log((Nao/Nai))); + const real Eks=RTONF*(log((Ko+pKNa*Nao)/(Ki+pKNa*Nai))); + const real Eca=0.5*RTONF*(log((Cao/Cai))); + real IKr; + real IKs; + real IK1; + real Ito; + real INa; + real IbNa; + real ICaL; + real IbCa; + real INaCa; + real IpCa; + real IpK; + real INaK; + real IKatp; + + real Ak1; + real Bk1; + real rec_iK1; + real rec_ipK; + real rec_iNaK; + real AM; + real BM; + real AH_1; + real BH_1; + real AH_2; + real BH_2; + real AJ_1; + real BJ_1; + real AJ_2; + real BJ_2; + real M_INF; + real H_INF; + real J_INF; + real TAU_M; + real TAU_H; + real TAU_J; + real axr1; + real bxr1; + real Xr1_INF; + real Xr2_INF_new; + real TAU_Xr1; + real Axs; + real Bxs; + real Xs_INF; + real TAU_Xs; + real R_INF_new; + real S_INF; + real TAU_S; + real Af; + real Bf; + real Cf; + real Af2; + real Bf2; + real Cf2; + real D_INF_new; + real TAU_F; + real F_INF; + real TAU_F2; + real F2_INF; + real sItot; + + + //Needed to compute currents + Ak1=0.1/(1.+exp(0.06*(svolt-Ek-200))); + Bk1=(3.*exp(0.0002*(svolt-Ek+100))+ + exp(0.1*(svolt-Ek-10)))/(1.+exp(-0.5*(svolt-Ek))); + rec_iK1=Ak1/(Ak1+Bk1); + rec_iNaK=(1./(1.+0.1245*exp(-0.1*svolt*F/(R*T))+0.0353*exp(-svolt*F/(R*T)))); + rec_ipK=1./(1.+exp((25-svolt)/5.98)); + + + //Compute currents + INa=GNa*sm*sm*sm*sh*sj*((svolt-Vm_modifier)-Ena); //ACIDOSIS + ICaL=GCaL*D_INF*sf*sf2*((svolt-Vm_modifier)-60); //ACIDOSIS + Ito=Gto*R_INF*ss*(svolt-Ek); + IKr=Gkr*sqrt(Ko/5.4)*sxr1*Xr2_INF*(svolt-Ek); + IKs=Gks*sxs*sxs*(svolt-Eks); + IK1=GK1*rec_iK1*(svolt-Ek); + INaCa=knaca*(1./(KmNai*KmNai*KmNai+Nao*Nao*Nao))*(1./(KmCa+Cao))* + (1./(1+ksat*exp((n-1)*svolt*F/(R*T))))* + (exp(n*svolt*F/(R*T))*Nai*Nai*Nai*Cao- + exp((n-1)*svolt*F/(R*T))*Nao*Nao*Nao*Cai*2.5); + + INaCa = INaCa*INaCa_multplicator; //ACIDOSIS + + INaK=knak*(Ko/(Ko+KmK))*(Nai/(Nai+KmNa))*rec_iNaK; + IpCa=GpCa*Cai/(KpCa+Cai); + IpK=GpK*rec_ipK*(svolt-Ek); + IbNa=GbNa*(svolt-Ena); + IbCa=GbCa*(svolt-Eca); + + IKatp = gkbaratp*(svolt-Ek) * Ikatp_multiplicator; + + + //Determine total current + (sItot) = IKr + + IKs + + IK1 + + Ito + + INa + + IbNa + + ICaL + + IbCa + + INaK + + INaCa + + IpCa + + IpK + + IKatp + + stim_current; + + //compute steady state values and time constants + AM=1./(1.+exp((-60.-svolt)/5.)); + BM=0.1/(1.+exp((svolt+35.)/5.))+0.10/(1.+exp((svolt-50.)/200.)); + TAU_M=AM*BM; + M_INF=1./((1.+exp((-56.86-svolt)/9.03))*(1.+exp((-56.86-svolt)/9.03))); + if (svolt>=-40.) + { + AH_1=0.; + BH_1=(0.77/(0.13*(1.+exp(-(svolt+10.66)/11.1)))); + TAU_H= 1.0/(AH_1+BH_1); + } + else + { + AH_2=(0.057*exp(-(svolt+80.)/6.8)); + BH_2=(2.7*exp(0.079*svolt)+(3.1e5)*exp(0.3485*svolt)); + TAU_H=1.0/(AH_2+BH_2); + } + H_INF=1./((1.+exp((svolt+71.55)/7.43))*(1.+exp((svolt+71.55)/7.43))); + if(svolt>=-40.) + { + AJ_1=0.; + BJ_1=(0.6*exp((0.057)*svolt)/(1.+exp(-0.1*(svolt+32.)))); + TAU_J= 1.0/(AJ_1+BJ_1); + } + else + { + AJ_2=(((-2.5428e4)*exp(0.2444*svolt)-(6.948e-6)* + exp(-0.04391*svolt))*(svolt+37.78)/ + (1.+exp(0.311*(svolt+79.23)))); + BJ_2=(0.02424*exp(-0.01052*svolt)/(1.+exp(-0.1378*(svolt+40.14)))); + TAU_J= 1.0/(AJ_2+BJ_2); + } + J_INF=H_INF; + + Xr1_INF=1./(1.+exp((-26.-svolt)/7.)); + axr1=450./(1.+exp((-45.-svolt)/10.)); + bxr1=6./(1.+exp((svolt-(-30.))/11.5)); + TAU_Xr1=axr1*bxr1; + Xr2_INF_new=1./(1.+exp((svolt-(-88.))/24.)); + + + Xs_INF=1./(1.+exp((-5.-svolt)/14.)); + Axs=(1400./(sqrt(1.+exp((5.-svolt)/6)))); + Bxs=(1./(1.+exp((svolt-35.)/15.))); + TAU_Xs=Axs*Bxs+80; + + if (transmurality == EPI) { + R_INF_new=1./(1.+exp((20-svolt)/6.)); + S_INF=1./(1.+exp((svolt+20)/5.)); + TAU_S=85.*exp(-(svolt+45.)*(svolt+45.)/320.)+5./(1.+exp((svolt-20.)/5.))+3.; + } + else if (transmurality == ENDO) { + R_INF_new=1./(1.+exp((20-svolt)/6.)); + S_INF=1./(1.+exp((svolt+28)/5.)); + TAU_S=1000.*exp(-(svolt+67)*(svolt+67)/1000.)+8.; + } + else if (transmurality == MID) { + R_INF_new=1./(1.+exp((20-svolt)/6.)); + S_INF=1./(1.+exp((svolt+20)/5.)); + TAU_S=85.*exp(-(svolt+45.)*(svolt+45.)/320.)+5./(1.+exp((svolt-20.)/5.))+3.; + } + + D_INF_new=1./(1.+exp((-8-svolt)/7.5)); + F_INF=1./(1.+exp((svolt+20)/7)); + Af=1102.5*exp(-(svolt+27)*(svolt+27)/225); + Bf=200./(1+exp((13-svolt)/10.)); + Cf=(180./(1+exp((svolt+30)/10)))+20; + TAU_F=Af+Bf+Cf; + F2_INF=0.67/(1.+exp((svolt+35)/7))+0.33; + Af2=600*exp(-(svolt+27)*(svolt+27)/170); + Bf2=7.75/(1.+exp((25-svolt)/10)); + Cf2=16/(1.+exp((svolt+30)/10)); + TAU_F2=Af2+Bf2+Cf2; + + //update voltage + rDY_[0] = -sItot; + + //Update gates + rDY_[1] = M_INF-(M_INF-sm)*exp(-dt/TAU_M); + rDY_[2] = H_INF-(H_INF-sh)*exp(-dt/TAU_H); + rDY_[3] = J_INF-(J_INF-sj)*exp(-dt/TAU_J); + rDY_[4] = Xr1_INF-(Xr1_INF-sxr1)*exp(-dt/TAU_Xr1); + rDY_[5] = Xs_INF-(Xs_INF-sxs)*exp(-dt/TAU_Xs); + rDY_[6]= S_INF-(S_INF-ss)*exp(-dt/TAU_S); + rDY_[7] =F_INF-(F_INF-sf)*exp(-dt/TAU_F); + rDY_[8] =F2_INF-(F2_INF-sf2)*exp(-dt/TAU_F2); + + rDY_[9] = D_INF_new; + rDY_[10] = R_INF_new; + rDY_[11] = Xr2_INF_new; \ No newline at end of file diff --git a/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.cu b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.cu new file mode 100644 index 00000000..5dff738b --- /dev/null +++ b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.cu @@ -0,0 +1,197 @@ +#include "../../gpu_utils/gpu_utils.h" +#include "../../monodomain/constants.h" +#include "ten_tusscher_tt3_mixed_endo_mid_epi.h" +#include + +extern "C" SET_ODE_INITIAL_CONDITIONS_GPU(set_model_initial_conditions_gpu) { + + log_info("Using Ten & Tusscher 3 GPU model\n"); + + uint32_t num_volumes = solver->original_num_cells; + + // execution configuration + const int GRID = (num_volumes + BLOCK_SIZE - 1)/BLOCK_SIZE; + + size_t size = num_volumes*sizeof(real); + + check_cuda_error(cudaMallocPitch((void **) &(solver->sv), &pitch_h, size, (size_t )NEQ)); + check_cuda_error(cudaMemcpyToSymbol(pitch, &pitch_h, sizeof(size_t))); + + kernel_set_model_inital_conditions <<>>(solver->sv, NULL, num_volumes); + + check_cuda_error( cudaPeekAtLastError() ); + cudaDeviceSynchronize(); + return pitch_h; +} + +extern "C" SOLVE_MODEL_ODES(solve_model_odes_gpu) { + + size_t num_cells_to_solve = ode_solver->num_cells_to_solve; + uint32_t * cells_to_solve = ode_solver->cells_to_solve; + real *sv = ode_solver->sv; + real dt = ode_solver->min_dt; + uint32_t num_steps = ode_solver->num_steps; + + // execution configuration + const int GRID = ((int)num_cells_to_solve + BLOCK_SIZE - 1)/BLOCK_SIZE; + + size_t stim_currents_size = sizeof(real)*num_cells_to_solve; + size_t cells_to_solve_size = sizeof(uint32_t)*num_cells_to_solve; + + real *stims_currents_device; + check_cuda_error(cudaMalloc((void **) &stims_currents_device, stim_currents_size)); + check_cuda_error(cudaMemcpy(stims_currents_device, stim_currents, stim_currents_size, cudaMemcpyHostToDevice)); + + uint32_t *cells_to_solve_device = NULL; + if(cells_to_solve != NULL) { + check_cuda_error(cudaMalloc((void **) &cells_to_solve_device, cells_to_solve_size)); + check_cuda_error(cudaMemcpy(cells_to_solve_device, cells_to_solve, cells_to_solve_size, cudaMemcpyHostToDevice)); + } + + real *fibrosis_device = NULL; + real *fibrosis = NULL; + real *transmurality_device = NULL; + real *transmurality = NULL; + int num_extra_parameters = 8; + real extra_par[num_extra_parameters]; + + size_t extra_parameters_size = num_extra_parameters * sizeof(real); + + real *extra_parameters_device; + real fibs_size = num_cells_to_solve*sizeof(real); + + struct extra_data_for_tt3* extra_data_from_cpu = (struct extra_data_for_tt3*)ode_solver->ode_extra_data; + + bool deallocate = false; + + if(ode_solver->ode_extra_data) { + fibrosis = extra_data_from_cpu->fibrosis; + transmurality = extra_data_from_cpu->transmurality; + extra_par[0] = extra_data_from_cpu->atpi; + extra_par[1] = extra_data_from_cpu->Ko; + extra_par[2] = extra_data_from_cpu->Ki; + extra_par[3] = extra_data_from_cpu->Vm_modifier; + extra_par[4] = extra_data_from_cpu->GNa_multiplicator; + extra_par[5] = extra_data_from_cpu->GCaL_multiplicator; + extra_par[6] = extra_data_from_cpu->INaCa_multiplicator; + extra_par[7] = extra_data_from_cpu->Ikatp_multiplicator; + } + else { + extra_par[0] = 6.8f; + extra_par[1] = 5.4f; + extra_par[2] = 138.3f; + extra_par[3] = 0.0; + extra_par[4] = 1.0f; + extra_par[5] = 1.0f; + extra_par[6] = 1.0f; + extra_par[7] = 1.0f; + + fibrosis = (real*)malloc(fibs_size); + transmurality = (real*)malloc(fibs_size); + + for(uint64_t i = 0; i < num_cells_to_solve; i++) { + fibrosis[i] = 1.0; // Healthy + transmurality[i] = 0.0; // ENDO + } + + deallocate = true; + } + + check_cuda_error(cudaMalloc((void **) &extra_parameters_device, extra_parameters_size)); + check_cuda_error(cudaMemcpy(extra_parameters_device, extra_par, extra_parameters_size, cudaMemcpyHostToDevice)); + + check_cuda_error(cudaMalloc((void **) &fibrosis_device, fibs_size)); + check_cuda_error(cudaMemcpy(fibrosis_device, fibrosis, fibs_size, cudaMemcpyHostToDevice)); + + check_cuda_error(cudaMalloc((void **) &transmurality_device, fibs_size)); + check_cuda_error(cudaMemcpy(transmurality_device, transmurality, fibs_size, cudaMemcpyHostToDevice)); + + solve_gpu<<>>(dt, sv, stims_currents_device, cells_to_solve_device, num_cells_to_solve, num_steps, fibrosis_device, transmurality_device, extra_parameters_device); + + check_cuda_error( cudaPeekAtLastError() ); + + check_cuda_error(cudaFree(stims_currents_device)); + check_cuda_error(cudaFree(fibrosis_device)); + check_cuda_error(cudaFree(transmurality_device)); + check_cuda_error(cudaFree(extra_parameters_device)); + + if(cells_to_solve_device) check_cuda_error(cudaFree(cells_to_solve_device)); + + if(deallocate) { + free(fibrosis); + free(transmurality); + } +} + +__global__ void kernel_set_model_inital_conditions(real *sv, real*IC, int num_volumes) { + + // Thread ID + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + + if(threadID < num_volumes) { + + *((real *) ((char *) sv + pitch * 0) + threadID) = INITIAL_V; // V; millivolt + *((real *) ((char *) sv + pitch * 1) + threadID) = 0.0f; //M + *((real *) ((char *) sv + pitch * 2) + threadID) = 0.75; //H + *((real *) ((char *) sv + pitch * 3) + threadID) = 0.75; //J + *((real *) ((char *) sv + pitch * 4) + threadID) = 0.0f; //Xr1 + *((real *) ((char *) sv + pitch * 5) + threadID) = 0.0f; //Xs + *((real *) ((char *) sv + pitch * 6) + threadID) = 1.0; //S + *((real *) ((char *) sv + pitch * 7) + threadID) = 1.0; //F + *((real *) ((char *) sv + pitch * 8) + threadID) = 1.0; //F2 + *((real *) ((char *) sv + pitch * 9) + threadID) = 0.0; //D_INF + *((real *) ((char *) sv + pitch * 10) + threadID) = 0.0; //R_INF + *((real *) ((char *) sv + pitch * 11) + threadID) = 0.0; //Xr2_INF + } +} + +// Solving the model for each cell in the tissue matrix ni x nj +__global__ void solve_gpu(real dt, real *sv, real* stim_currents, + uint32_t *cells_to_solve, uint32_t num_cells_to_solve, + int num_steps, real *fibrosis, real *transmurality, real *extra_parameters) { + + int threadID = blockDim.x * blockIdx.x + threadIdx.x; + int sv_id; + + // Each thread solves one cell model + if(threadID < num_cells_to_solve) { + if(cells_to_solve) + sv_id = cells_to_solve[threadID]; + else + sv_id = threadID; + + real rDY[NEQ]; + + for (int n = 0; n < num_steps; ++n) { + + RHS_gpu(sv, rDY, stim_currents[threadID], sv_id, dt, fibrosis[threadID], transmurality[threadID], extra_parameters); + + *((real*)((char*)sv) + sv_id) = dt*rDY[0] + *((real*)((char*)sv) + sv_id); + + for(int i = 1; i < NEQ; i++) { + *((real*)((char*)sv + pitch * i) + sv_id) = rDY[i]; + } + } + } +} + + +inline __device__ void RHS_gpu(real *sv_, real *rDY_, real stim_current, int threadID_, real dt, real fibrosis, real transmurality, real *extra_parameters) { + + //fibrosis = 0 means that the cell is fibrotic, 1 is not fibrotic. Anything between 0 and 1 means border zone + //transmurality = 0 means cell is endocardium, 1 is mid-myocardium and 2 is epicardium + const real svolt = *((real*)((char*)sv_ + pitch * 0) + threadID_); + const real sm = *((real*)((char*)sv_ + pitch * 1) + threadID_); + const real sh = *((real*)((char*)sv_ + pitch * 2) + threadID_); + const real sj = *((real*)((char*)sv_ + pitch * 3) + threadID_); + const real sxr1 = *((real*)((char*)sv_ + pitch * 4) + threadID_); + const real sxs = *((real*)((char*)sv_ + pitch * 5) + threadID_); + const real ss = *((real*)((char*)sv_ + pitch * 6) + threadID_); + const real sf = *((real*)((char*)sv_ + pitch * 7) + threadID_); + const real sf2 = *((real*)((char*)sv_ + pitch * 8) + threadID_); + const real D_INF = *((real*)((char*)sv_ + pitch * 9) + threadID_); + const real R_INF = *((real*)((char*)sv_ + pitch * 10) + threadID_); + const real Xr2_INF = *((real*)((char*)sv_ + pitch * 11) + threadID_); + + #include "ten_tusscher_tt3_mixed_endo_mid_epi.common.c" +} diff --git a/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.h b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.h new file mode 100644 index 00000000..324cb662 --- /dev/null +++ b/src/models_library/ten_tusscher/ten_tusscher_tt3_mixed_endo_mid_epi.h @@ -0,0 +1,34 @@ +#ifndef MONOALG3D_MODEL_TEN_TUSSCHER_3_MIXED_ENDO_MID_EPI_H +#define MONOALG3D_MODEL_TEN_TUSSCHER_3_MIXED_ENDO_MID_EPI_H + +#include "../model_common.h" +#include "../../extra_data_library/helper_functions.h" + +#define NEQ 12 +#define INITIAL_V (-86.2f) + +#define ENDO 0.0 +#define MID 1.0 +#define EPI 2.0 + +#ifdef __CUDACC__ + +#include "../../gpu_utils/gpu_utils.h" + +static __device__ size_t pitch; +static size_t pitch_h; + +__global__ void kernel_set_model_inital_conditions(real *sv, real* ICs, int num_volumes); + +__global__ void solve_gpu(real dt, real *sv, real *stim_currents, uint32_t *cells_to_solve, uint32_t num_cells_to_solve, + int num_steps, real *fibrosis, real *transmurality, real *extra_parameters); + +inline __device__ void RHS_gpu(real *sv_, real *rDY_, real stim_current, int threadID_, real dt, real fibrosis, real transmurality, + real *extra_parameters); + +#endif + +void RHS_cpu(const real *sv, real *rDY_, real stim_current, real dt, real fibrosis, real transmurality, real const *extra_parameters); +void solve_model_ode_cpu(real dt, real *sv, real stim_current, real fibrosis, real transmurality, real *extra_parameters); + +#endif // MONOALG3D_MODEL_TEN_TUSSCHER_3_MIXED_ENDO_MID_EPI_H