Correct usage of adjustments files to apply ionic gradient in biventricular simulations

Question

Hi, I am doing a monodomain simulation on a biventricular mesh cut at the base using the Tomek model, I have already tagged 4 regions of the mesh based on the transmural coordinates and applied to them flags for the EPI, MCELL and ENDO ionic variant of the model. I want also to apply a gradient for the GKs_b variable of the model based on the apicobasal and transmural coordinates. 
To do so since the original Tomek in the limpet library does not expose any variable, I modified and recompiled the Tomek.model to expose the GKs_b var as nodal value and I created an adjustment file for mesh nodes based on ab coordinates. Then I applied the adj file with the ajustment parameters, the following is part of the carputils script with imp regions and phys region.

    cmd += [
        '-num_imp_regions', 3,
        '-imp_region[0].name', 'Endocardio',
        '-imp_region[0].num_IDs', 2,
        '-imp_region[0].ID[0]', 3,      #tag endocardio 3 
        '-imp_region[0].ID[1]', 4,      #tag fast endocardio 4
        '-imp_region[0].im', 'Tomek_editv3',
        '-imp_region[0].im_param', 'flags=ENDO',
       
        '-imp_region[1].name', 'Mid_Miocardio',
        '-imp_region[1].num_IDs', 1,
        '-imp_region[1].ID[0]', 2,  
        '-imp_region[1].im', 'Tomek_editv3',
        '-imp_region[1].im_param', 'flags=MCELL',

        '-imp_region[2].name', 'Epicardio',
        '-imp_region[2].num_IDs', 1,
        '-imp_region[2].ID[0]', 1,    
        '-imp_region[2].im', 'Tomek_editv3',
        '-imp_region[2].im_param', 'flags=EPI',
        '-num_adjustments', 1,   
        '-adjustment[0].variable', 'Tomek_editv3.GKs_b',
        '-adjustment[0].file', './sb3901_meshes/gks_tomek_gradient30abtm.adj',
        '-adjustment[0].dump', 1,
    ]
I added the flag dump, the simulation export also 2 dat files Gks_b_canonical.dat and Gks_perm.dat However If I open them on the input mesh in Meshalyzer, one shows only a portion of the mesh with the gradient values the rest is black while the other file shows various colored points distributed non uniformally on the mesh. 

I don't understand if I am making mistakes in the generation and application of adj file or it is a problem of the dump option.
At the folder at this link I uploaded the the file used for the simulation and the dat files output from the simulation. opencarp_simulation_files

Thanks for the help!

Answer 1

Hey,

The _perm.dat file is in PETSc parallel ordering by design, it's a debug artefact and is expected to look scrambled on the mesh. Only _canonical.dat is meant to be displayed.

Can you send the simulator stdout line Adjusted N values for Tomek_editv3.GKs_b. This tells us whether the lookup itself dropped most indices (N = 126323) or whether only the dump path is broken (N = 1429736).

Comments

Hi, the stdout log does not show that line, only assignment to tagged regions are present. I tried also to increase the output level to 10 but it does not help.
Here is the full stout log if it can help

Command Line Options
----------------------------------------------------------------------------------------------------------
   NP:        12
   BUILD:     CPU
   FLAVOR:    petsc
   DRY:       False
   VISUALIZE: False
----------------------------------------------------------------------------------------------------------
Solver Option Files
----------------------------------------------------------------------------------------------------------
   Parabolic: /usr/local/lib/python3.10/dist-packages/carputils/resources/petsc_options/ilu_cg_opts
   Elliptic:  /usr/local/lib/python3.10/dist-packages/carputils/resources/petsc_options/gamg_cg_opts
-----------------------------------------------------------------------------------------------------------------------------------
Paths
-----------------------------------------------------------------------------------------------------------------------------------
   OPENCARP_SRC:  /usr/local/lib
   MPIEXEC:       /usr/local/lib/opencarp/lib/petsc/bin/mpiexec
-----------------------------------------------------------------------------------------------------------------------------------
System Info
-----------------------------------------------------------------------------------------------------------------------------------
   Python:    3.10.12
   OS:        Ubuntu 22.04.5 LTS
   SETTINGS:  /.config/carputils/settings.yaml
   Email:     you@example.com
   Name:      Your Name
   ORCID:     0000-0000-0000-0000
============================================================
Leggendo i nodi da ./sb3901_rootnodes/sb3901_fine500um_candidate_root_nodes.vtx e i tempi da ./sb3901_rootnodes/sb3901_candidate_root_nodes_times.csv...
mkdir -p /shared/sb3901/test_mono_torord_sb3901_norm5

#--------------------------------------------------------------------------------------------------------
#                                  Launching openCARP Simulation ./test_mono_torord_sb3901_norm5
#--------------------------------------------------------------------------------------------------------

mpiexec -n 12 /usr/local/bin/openCARP \
  -parab_options_file /usr/local/lib/python3.10/dist-packages/carputils/resources/petsc_options/ilu_cg_opts \
  -ellip_options_file /usr/local/lib/python3.10/dist-packages/carputils/resources/petsc_options/gamg_cg_opts \
  -simID ./test_mono_torord_sb3901_norm5 \
  -meshname ./sb3901_meshes/sb3901_500um_taggedv2_fastendo \
  -tend 600.0 \
  -spacedt 1.0 \
  -timedt 1.0 \
  -dt 10.0 \
  -bidomain 0 \
  -parab_solve 1 \
  -mass_lumping 1 \
  -output_level 5 \
  -phie_rec_ptf ./sb3901_meshes/sb3901_electrodesum \
  -phie_recovery_file sb3901_phie_recovery_norm2_prep_fast \
  -num_LATs 1 \
  -lats[0].ID ACT \
  -lats[0].measurand 0 \
  -lats[0].all 0 \
  -lats[0].threshold -10.0 \
  -lats[0].mode 0 \
  -num_external_imp 1 \
  -external_imp[0] /shared/tomek-model/Tomek_editv3.so \
  -num_gregions 2 \
  -gregion[0].name Miocardio \
  -gregion[0].num_IDs 3 \
  -gregion[0].ID[0] 1 \
  -gregion[0].ID[1] 2 \
  -gregion[0].ID[2] 3 \
  -gregion[0].g_il 0.2993 \
  -gregion[0].g_it 0.1114 \
  -gregion[0].g_in 0.0601 \
  -gregion[1].name FastEndo \
  -gregion[1].num_IDs 1 \
  -gregion[1].ID[0] 4 \
  -gregion[1].g_il 1.1972 \
  -gregion[1].g_it 0.4456 \
  -gregion[1].g_in 0.2404 \
  -num_imp_regions 3 \
  -imp_region[0].name Endocardio \
  -imp_region[0].num_IDs 2 \
  -imp_region[0].ID[0] 3 \
  -imp_region[0].ID[1] 4 \
  -imp_region[0].im Tomek_editv3 \
  -imp_region[0].im_param "flags=ENDO" \
  -imp_region[1].name Mid_Miocardio \
  -imp_region[1].num_IDs 1 \
  -imp_region[1].ID[0] 2 \
  -imp_region[1].im Tomek_editv3 \
  -imp_region[1].im_param "flags=MCELL" \
  -imp_region[2].name Epicardio \
  -imp_region[2].num_IDs 1 \
  -imp_region[2].ID[0] 1 \
  -imp_region[2].im Tomek_editv3 \
  -imp_region[2].im_param "flags=EPI" \
  -num_adjustments 1 \
  -adjustment[0].variable Tomek_editv3.GKs_b \
  -adjustment[0].file ./sb3901_meshes/gks_tomek_gradient30abtm.adj \
  -adjustment[0].dump 1 \
  -num_phys_regions 1 \
  -phys_region[0].ptype 0 \
  -phys_region[0].name Intracellular domain \
  -phys_region[0].num_IDs 4 \
  -phys_region[0].ID[0] 1 \
  -phys_region[0].ID[1] 2 \
  -phys_region[0].ID[2] 3 \
  -phys_region[0].ID[3] 4 \
  -num_stim 13 \
  -stim[0].name RN1 \
  -stim[0].elec.vtx_file ./test_mono_torord_sb3901_norm5/sb3901_rootnodes_vtx/RN1.vtx \
  -stim[0].ptcl.duration 4.0 \
  -stim[0].ptcl.npls 1 \
  -stim[0].ptcl.start 26.3145 \
  -stim[0].pulse.strength 53 \
  -stim[0].crct.type 0 \
  -stim[1].name RN2 \
  -stim[1].elec.vtx_file ./test_mono_torord_sb3901_norm5/sb3901_rootnodes_vtx/RN2.vtx \
  -stim[1].ptcl.duration 4.0 \
  -stim[1].ptcl.npls 1 \
  -stim[1].ptcl.start 35.982 \
  -stim[1].pulse.strength 53 \
  -stim[1].crct.type 0 \
  -stim[2].name RN3 \
  -stim[2].elec.vtx_file ./test_mono_torord_sb3901_norm5/sb3901_rootnodes_vtx/RN3.vtx \
  -stim[2].ptcl.duration 4.0 \
  -stim[2].ptcl.npls 1 \
  -stim[2].ptcl.start 44.987 \
  -stim[2].pulse.strength 53 \
  -stim[2].crct.type 0 \
.......

*** GIT tag:            4dd37e6
*** GIT hash:           4dd37e6fc547565cf18a86a1931bca490faf90d8
*** GIT repo:           https://gitlab-ci-token:glcbt-64_rsFoNNCFNEGDU5R7xhGV@git.opencarp.org/openCARP/openCARP.git
*** dependency commits:

L2 : Output directory exists: ./test_mono_torord_sb3901_norm5

    *** Processing meshes ***

Reading reference mesh: ./sb3901_meshes/sb3901_500um_taggedv2_fastendo.*
Done in 35.292099 sec.
Done in 0.002841 sec.
Processing submeshes

Submesh name: Intracellular domain
Extraction done in 1.617777 sec.

Submesh name: Intracellular domain
Using KDtree partitioner ..
Partitioning done in 4.773577 sec.
Redistributing done in 0.123178 sec.
Canonical numbering done in 0.171225 sec.
PETSc numbering done in 0.925716 sec.
===== Parallel mesh statistics =====
#pid    #nodes  #elems  #interf #alg
0       135082  698892  8759    126323
1       132283  698888  14327   123260
2       131987  698887  13924   123256
3       125684  698888  9905    118774
4       123628  667121  15817   117725
5       127668  667118  18068   117725
6       122721  667122  12185   117117
7       124120  667119  15260   117117
8       121580  667119  11400   117112
9       121109  667121  12413   117112
10      124005  667120  17754   117109
11      121032  667114  10913   117106

Multiplicities :
2: 155942
3: 4743
4: 40
All done in 45.818863 sec.

    *** Initializing physics ***

Initializing Electrics ..

Setting up ionic models and plugins
-----------------------------------

Assigning IMPS to tagged regions:
        Ionic model: Tomek_editv3 to tag region(s) [3], [4].
        Ionic model: Tomek_editv3 to tag region(s) [2].
        Ionic model: Tomek_editv3 to tag region(s) [1].
Ionic model: Tomek_editv3

Ionic model: Tomek_editv3

Ionic model: Tomek_editv3

Stimulus 0: Selecting vertices from file ./test_mono_torord_sb3901_norm5/sb3901_rootnodes_vtx/RN1.vtx
L3 : Warning: openCARP requires input vtx indices to be of from the input mesh, i.e. of type "extra".
Stimulus 1: Selecting vertices from file ./test_mono_torord_sb3901_norm5/sb3901_rootnodes_vtx/RN2.vtx
L3 : Warning: openCARP requires input vtx indices to be of from the input mesh, i.e. of type "extra".
......
Computed parabolic stiffness matrix in 4.494 seconds.
Computed parabolic mass matrix in 1.321 seconds.
Initializing parabolic solver in 0.00059 seconds.

---

Ok, this helps a lot already. It seems that there is a bug with the dump option. Your log shows rank 0 has #alg = 126323, which exactly matches the non-zero count in your dump files. It looks like all ranks write to rank 0 global indices and the last writes stick. I will open an issue for that.

The "Adjusted N values for Tomek_editv3.GKs_b" line is still important as it confirms that the simulation still runs fine. It actually should be in the electrics.log file, not in stdout output as I said before. If you find it there with N == 1,429,736, your simulation is fine and the dump output is misleading due to this potential bug.

---

Ok, I checked the electrics.log file and I actually found the line saying "Adjusted 1429736 values for Tomek_editv3.GKs_b " so it is confirmed that it is a bug. Thanks for the help!

---

Another problem always related to MPI ranks that I encountered during the simulations is that I can't go over 12 ranks (--np 12 ), because the simulation randomly crash with the error "MPI_Abort(compute, 59) - process 12"
I am running the simulations using the official opencarp docker container on an Ubuntu machine with 24 cores and 64 gb of ram. It could be a problem related to my system or the container or a problem with the mesh partitioning?