Hi everyone,
I am trying to run a simulation in OpenCARP using a spherical mesh with volume instead of the typical block. So far, I have followed these steps:
- I generated a spherical mesh with volume using MATLAB and exported the .elem, .pts, and .lon files.
- I loaded the mesh into OpenCARP and attempted to run the simulation.
- However, I encountered an error, and the simulation does not run correctly.
Here is the error message from the output:
NaN detected: gregion[0], element 330: 587 572 585 564
Volume: -nan, fib: (1, 0, 0), she: (0, 0, 0), G: (0.136108, 0.0175843, 0.0175843)
Below are the configuration files I am using:
- monodomain.par:
# ionic setup
num_imp_regions = 1
imp_region[0].im = LuoRudy91
# electrical setup
num_stim = 1
stimulus[0].name = "S1"
stimulus[0].stimtype = 0
stimulus[0].strength = 400.0 # Stimulus strength in microamperes
stimulus[0].duration = 2.0
stimulus[0].start = 20.0 # First stimulus starts at 20 ms
bidomain = 0
parab_solve = 1
- run.py:
#!/usr/bin/env python
import os
from datetime import date
from carputils import settings
from carputils import tools
from carputils import mesh
from carputils.carpio import txt
import numpy as np
def parser():
parser = tools.standard_parser()
group = parser.add_argument_group('experiment specific options')
group.add_argument('--duration',
type=float,
default=1000.,
help='Duration of simulation in [ms] (default: 1000.)')
group.add_argument('--S1-strength',
type=float,
default=400., # Stimulus strength in microamperes
help='Stimulus strength in [uA/cm^2] (default: 400.)')
group.add_argument('--S1-dur',
type=float,
default=2., # Stimulus duration in ms
help='Stimulus duration in [ms] (default: 2.)')
return parser
def jobID(args):
"""
Generate name of top-level output directory.
"""
today = date.today()
return '{}_monodomain_{}'.format(today.isoformat(), args.duration)
@tools.carpexample(parser, jobID)
def run(args, job):
# Generate mesh
meshname = 'sphere'
# Query for element labels
_, etags, _ = txt.read(meshname + '.elem')
etags = np.unique(etags)
IntraTags = etags[etags != 0].tolist()
ExtraTags = etags.tolist().copy()
# Calculate number of stimuli
num_stim = int((args.duration - 20) / 500) + 1 # Includes initial stimulus at 20 ms
# Create stimulus settings
stim = ['-num_stim', num_stim]
for i in range(num_stim):
stim += [
f'-stimulus[{i}].name', 'S1',
f'-stimulus[{i}].stimtype', 0,
f'-stimulus[{i}].strength', args.S1_strength,
f'-stimulus[{i}].duration', args.S1_dur,
f'-stimulus[{i}].start', 20 + i * 500 # Stimulus every 500 ms starting at 20 ms
]
# Generate CARP command
cmd = tools.carp_cmd(os.path.join(os.path.dirname(__file__), 'monodomain.par'))
cmd += ['-simID', job.ID,
'-meshname', meshname,
'-dt', 5,
'-tend', args.duration]
cmd += tools.gen_physics_opts(ExtraTags=ExtraTags, IntraTags=IntraTags)
cmd += stim
if args.visualize:
cmd += ['-gridout_i', 3, '-gridout_e', 3, '-spacedt', 0.1]
# Run simulation
job.carp(cmd)
# Visualization (optional)
if args.visualize and not settings.platform.BATCH:
geom = os.path.join(job.ID, os.path.basename(meshname) + '_i')
data = os.path.join(job.ID, 'vm.igb.gz')
view = os.path.join(os.path.dirname(__file__), 'view_vm.mshz')
job.meshalyzer(geom, data, view)
if __name__ == '__main__':
run()
I have checked the OpenCARP documentation and examples, but I haven't found any that use a spherical mesh instead of a block.
Questions:
- Has anyone successfully run a simulation on a volumetric spherical mesh?
- Is there an example in OpenCARP that uses a sphere instead of a block?
Any advice or references would be greatly appreciated.
Thanks in advance!
