diff --git a/.github/workflows/check-pep8.yml b/.github/workflows/check-pep8.yml index e89bbbbae..0bd0d44c9 100644 --- a/.github/workflows/check-pep8.yml +++ b/.github/workflows/check-pep8.yml @@ -15,7 +15,7 @@ jobs: - uses: actions/checkout@v2 - name: autopep8 id: autopep8 - uses: peter-evans/autopep8@v1 + uses: peter-evans/autopep8@v2 with: args: --recursive --diff --aggressive --aggressive --exit-code --ignore E402 --max-line-length 120 . - name: Fail if autopep8 made changes diff --git a/muscle-tendon-complex/.gitignore b/muscle-tendon-complex/.gitignore new file mode 100644 index 000000000..03904b82a --- /dev/null +++ b/muscle-tendon-complex/.gitignore @@ -0,0 +1,13 @@ +**/__pycache__ +**/lib +**/logs +**/out +**/*.log +**/*.txt +**/precice-profiling +**/build_release +**/muscle-solver +**/tendon-solver +**/.sconf_temp +**/.scons* +muscle-opendihu/src diff --git a/muscle-tendon-complex/README.md b/muscle-tendon-complex/README.md new file mode 100644 index 000000000..685027e15 --- /dev/null +++ b/muscle-tendon-complex/README.md @@ -0,0 +1,125 @@ +--- +title: Muscle-tendon complex +permalink: tutorials-muscle-tendon-complex.html +keywords: multi-coupling, OpenDiHu, skeletal muscle +summary: In this case, a skeletal muscle (biceps) and three tendons are coupled together using a fully-implicit multi-coupling scheme. +--- + +{% note %} +Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/muscle-tendon-complex). Read how in the [tutorials introduction](https://www.precice.org/tutorials.html). +{% endnote %} + +## Case Setup + +In the following tutorial, we model the contraction of a muscle (the biceps). The biceps is attached to the bones by three tendons (one at the bottom and two at the top). We enforce an activation in the muscle which results in its contraction. The tendons move as a result of the muscle contraction. In this case, a muscle and three tendons are coupled together using a fully-implicit multi-coupling scheme. The case setup is shown in the following figure: + +![Setup](images/tutorials-muscle-tendon-complex-setup.png) + +The muscle participant (in red) is connected to three tendons. The muscle sends traction values to the tendons, which send displacement and velocity values back to the muscle. The end of each tendon which is not attached to the muscle is fixed by a dirichlet boundary condition (in reality, it would be fixed to the bones). + +The muscle and tendon meshes are obtained from patient imaging. The interfaces of the tendons and the muscle do not perfectly match, which is a quite common issue due to the limitations of imaging methods and postprocessing tools. Nonetheless, preCICE coupling methods are robust and can handle meshes that do not perfectly match. + +## Why multi-coupling? + +This is a case with four participants: the muscle and each tendon. In preCICE, there are two options to [couple more than two participants](https://www.precice.org/configuration-coupling-multi.html). The first option is a composition of bi-coupling schemes, in which we must specify the exchange of data in a participant-to-participant manner, limited to primarily explicit coupling schemes. However, such a composition is not suited for combining multiple strong interactions [1]. Thus, in this case, we use the second option, fully-implicit multi-coupling. For another multi-coupling tutorial, you can refer to the [multiple perpendicular flaps tutorial](http://precice.org/tutorials-multiple-perpendicular-flaps.html). + +We can set this in our `precice-config.xml`: + +```xml + + + + + +``` + +The participant that has the control is the one that it is connected to all other participants. This is why we have chosen the muscle participant for this task. + +## About the solvers + +We use solvers based on [OpenDiHu](https://github.com/opendihu/opendihu) for all participants. + +**The muscle solver** consists of a multi-physcis multi-scale solver itself. It combines two OpenDiHu solvers in one: the *FastMonodomainSolver* and the *MuscleContractionSolver*. The two solvers are coupled using the OpenDiHu coupling tool for weak coupling. + +- The [FastMonodomainSolver](https://opendihu.readthedocs.io/en/latest/settings/fast_monodomain_solver.html) models the electrochemical processes that take place in the muscle fibers. Motor neurons fire electrical signals that propagate from the neuromuscular junctions. i.e. the center of the muscle fibers, to the extremes of the muscle fibers. The propagation of the electrical signal triggers chemical reactions which lead to the contraction of sarcomeres, the smallest contraction unit in the muscle. +The solver solves the so called "monodomain equation" independently for each fiber [2]. The modonomain equation has a reaction term (small time scale) and a diffusion term (large time scale) and is solved using Strang splitting. The sarcomeres, i.e., the reaction term, are modelled using a variant of the Shorten model, specified by the CellML file `opendihu/examples/electrophysiology/input/2020_06_03_hodgkin-huxley_shorten_ocallaghan_davidson_soboleva_2007.cellml`. The firing of the motor neurons is modelled by the `opendihu/examples/electrophysiology/input/MU_firing_times_always.txt` file. + +- The [MuscleContractionSolver](https://opendihu.readthedocs.io/en/latest/settings/muscle_contraction_solver.html) models the mechanics of the muscle. It consists of a dynamic FEM solver that models an hyperelastic active material. The active component is calculated from the activation parameter $\gamma$, which ranges from 0 (no activation) to 1 (maximum activation) and is calculated in the *FastMonodomainSolver*. The material parameters are chosen as in [Heidlauf et al.](https://link.springer.com/article/10.1007/s10237-016-0772-7) + +**The tendon solver** is a dynamic FEM mechanical solver. It models an hyperelastic passive material. The material parameters are chosen as in [Carniel et al.](https://pubmed.ncbi.nlm.nih.gov/28238424/) + +## Running the Simulation + +1. Preparation: + - Install OpenDiHu + + In the OpenDiHu website you can find detailed [installation instructions](https://opendihu.readthedocs.io/en/latest/user/installation.html). + We recommend to download the code from the [GitHub repository](https://github.com/opendihu/opendihu) and to run `make release_without_tests` in the parent directory. + + {% note %} + OpenDiHu automatically downloads dependencies and installs them in the `opendihu/dependencies/` folder. You can avoid that by setting e.g., `PRECICE_DOWNLOAD = False` in the [user-variables.scons.py](https://github.com/opendihu/opendihu/blob/develop/user-variables.scons.py) before building OpenDiHu. + {% endnote %} + + - Download input files for OpenDiHu + + OpenDiHu requires input files hosted in [Zenodo](https://zenodo.org/records/4705982) which include CellML files (containing model equations) and mesh files. Downloading these files is necessary to simulate muscles and/or tendons with OpenDiHu. You can [directly download the necessary files](https://zenodo.org/record/4705982/files/input.tgz?download=1). + + - Setup `$OPENDIHU_HOME` and `$OPENDIHU_INPUT_DIR` in your `.bashrc` file + + ```bash + export OPENDIHU_HOME=/path/to/opendihu + export OPENDIHU_INPUT_DIR=/path/to/input + ``` + + - Compile muscle and tendon solvers + + ```bash + cd solver-opendihu + ./build.sh + ``` + +2. Starting the simulation: + + We are going to run each solver in a different terminal. It is important that first we navigate to the respective case directory, so that all solvers start from directories with same parent directory (see `exchange-directory` in the `precice-config.xml`). + To start the `Muscle` participant, run: + + ```bash + cd muscle-opendihu + ./run.sh + ``` + + To start the `Tendon-Bottom` participant, run: + + ```bash + cd tendon-bottom-opendihu + ./run.sh + ``` + + To start the `Tendon-Top-A` participant, run: + + ```bash + cd tendon-top-A-opendihu + ./run.sh + ``` + + Finally, to start the `Tendon-Top-B` participant, run: + + ```bash + cd tendon-top-B-opendihu + ./run.sh + ``` + +## Postprocessing... TODO + +After the simulation has finished, you can visualize your results using e.g. ParaView. + +## References TODO + + +[1] H. Bungartz, F. Linder, M. Mehl, B. Uekermann. A plug-and-play coupling approach for parallel multi-field simulations. *Comput Mech* **55**, 1119-1129 (2015). https://doi.org/10.1007/s00466-014-1113-2 + +[2] T. Heidlauf, O. Röhrle. A multiscale chemo-electro-mechanical skeletal muscle model to analyze muscle contraction and force generation for different muscle fiber arrangements. *Front. Physiology* **5** (2014). https://doi.org/10.3389/fphys.2014.00498 + +{% disclaimer %} +This offering is not approved or endorsed by OpenCFD Limited, producer and distributor of the OpenFOAM software via www.openfoam.com, and owner of the OPENFOAM® and OpenCFD® trade marks. +{% enddisclaimer %} diff --git a/muscle-tendon-complex/clean-tutorial.sh b/muscle-tendon-complex/clean-tutorial.sh new file mode 120000 index 000000000..4713f5092 --- /dev/null +++ b/muscle-tendon-complex/clean-tutorial.sh @@ -0,0 +1 @@ +../tools/clean-tutorial-base.sh \ No newline at end of file diff --git a/muscle-tendon-complex/images/tutorials-muscle-tendon-complex-setup.png b/muscle-tendon-complex/images/tutorials-muscle-tendon-complex-setup.png new file mode 100644 index 000000000..d744e1f24 Binary files /dev/null and b/muscle-tendon-complex/images/tutorials-muscle-tendon-complex-setup.png differ diff --git a/muscle-tendon-complex/muscle-opendihu/clean.sh b/muscle-tendon-complex/muscle-opendihu/clean.sh new file mode 100755 index 000000000..56411170b --- /dev/null +++ b/muscle-tendon-complex/muscle-opendihu/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_opendihu . diff --git a/muscle-tendon-complex/muscle-opendihu/helper.py b/muscle-tendon-complex/muscle-opendihu/helper.py new file mode 100644 index 000000000..c01343ca1 --- /dev/null +++ b/muscle-tendon-complex/muscle-opendihu/helper.py @@ -0,0 +1,606 @@ +# Multiple 1D fibers (monodomain) with 3D contraction, biceps geometry +# This is a helper script that sets a lot of the internal variables which are all defined in variables.py +# +# if variables.fiber_file=cuboid.bin, it uses a small cuboid test example + +import numpy as np +import pickle +import sys +import os +import struct +import argparse +# sys.path.insert(0, '..') +import variables # file variables.py +from create_partitioned_meshes_for_settings import * # file create_partitioned_meshes_for_settings + +# parse arguments +rank_no = (int)(sys.argv[-2]) +n_ranks = (int)(sys.argv[-1]) + +# generate cuboid fiber file +if "cuboid.bin" in variables.fiber_file: + + if variables.n_fibers_y is None: + variables.n_fibers_x = 4 + variables.n_fibers_y = variables.n_fibers_x + variables.n_points_whole_fiber = 20 + + size_x = variables.n_fibers_x * 0.1 + size_y = variables.n_fibers_y * 0.1 + size_z = variables.n_points_whole_fiber / 100. + + if rank_no == 0: + print( + "create cuboid.bin with size [{},{},{}], n points [{},{},{}]".format( + size_x, + size_y, + size_z, + variables.n_fibers_x, + variables.n_fibers_y, + variables.n_points_whole_fiber)) + + # write header + with open(variables.fiber_file, "wb") as outfile: + + # write header + header_str = "opendihu self-generated cuboid " + outfile.write(struct.pack('32s', bytes(header_str, 'utf-8'))) # 32 bytes + outfile.write(struct.pack('i', 40)) # header length + outfile.write(struct.pack('i', variables.n_fibers_x * variables.n_fibers_y)) # n_fibers + outfile.write(struct.pack('i', variables.n_points_whole_fiber)) # variables.n_points_whole_fiber + outfile.write(struct.pack('i', 0)) # nBoundaryPointsXNew + outfile.write(struct.pack('i', 0)) # nBoundaryPointsZNew + outfile.write(struct.pack('i', 0)) # nFineGridFibers_ + outfile.write(struct.pack('i', 1)) # nRanks + outfile.write(struct.pack('i', 1)) # nRanksZ + outfile.write(struct.pack('i', 0)) # nFibersPerRank + outfile.write(struct.pack('i', 0)) # date + + # loop over points + for y in range(variables.n_fibers_y): + for x in range(variables.n_fibers_x): + for z in range(variables.n_points_whole_fiber): + point = [x * (float)(size_x) / (variables.n_fibers_x), y * (float)(size_y) / + (variables.n_fibers_y), z * (float)(size_z) / (variables.n_points_whole_fiber)] + outfile.write(struct.pack('3d', point[0], point[1], point[2])) # data point + +# output diffusion solver type +if rank_no == 0: + print("diffusion solver type: {}".format(variables.diffusion_solver_type)) + +variables.load_fiber_data = False # load all local node positions from fiber_file, in order to infer partitioning for fat_layer mesh + +# create the partitioning using the script in create_partitioned_meshes_for_settings.py +result = create_partitioned_meshes_for_settings( + variables.n_subdomains_x, + variables.n_subdomains_y, + variables.n_subdomains_z, + variables.fiber_file, + variables.load_fiber_data, + variables.sampling_stride_x, + variables.sampling_stride_y, + variables.sampling_stride_z, + variables.generate_linear_3d_mesh, + variables.generate_quadratic_3d_mesh) +[variables.meshes, + variables.own_subdomain_coordinate_x, + variables.own_subdomain_coordinate_y, + variables.own_subdomain_coordinate_z, + variables.n_fibers_x, + variables.n_fibers_y, + variables.n_points_whole_fiber] = result + +variables.n_subdomains_xy = variables.n_subdomains_x * variables.n_subdomains_y +variables.n_fibers_total = variables.n_fibers_x * variables.n_fibers_y + +# create mappings between meshes +# variables.mappings_between_meshes = {"MeshFiber_{}".format(i) : "3Dmesh" for i in range(variables.n_fibers_total)} +variables.mappings_between_meshes = {"MeshFiber_{}".format( + i): {"name": "3Dmesh", "xiTolerance": 1e-3} for i in range(variables.n_fibers_total)} + +# a higher tolerance includes more fiber dofs that may be almost out of the 3D mesh +variables.mappings_between_meshes = { + "MeshFiber_{}".format(i): { + "name": "3Dmesh_quadratic", + "xiTolerance": variables.mapping_tolerance, + "enableWarnings": False, + "compositeUseOnlyInitializedMappings": False, + "fixUnmappedDofs": True, + "defaultValue": 0, + } for i in range(variables.n_fibers_total) +} +# set output writer +variables.output_writer_fibers = [] +variables.output_writer_elasticity = [] +variables.output_writer_emg = [] +variables.output_writer_0D_states = [] + +subfolder = "" +if variables.paraview_output: + if variables.adios_output: + subfolder = "paraview/" + variables.output_writer_emg.append({"format": "Paraview", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D_emg), + "filename": "out/" + subfolder + variables.scenario_name + "/hd_emg", + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"}) + variables.output_writer_elasticity.append({"format": "Paraview", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D), + "filename": "out/" + subfolder + variables.scenario_name + "/elasticity", + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"}) + variables.output_writer_fibers.append({"format": "Paraview", + "outputInterval": int(1. / variables.dt_splitting * variables.output_timestep_fibers), + "filename": "out/" + subfolder + variables.scenario_name + "/fibers", + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"}) + if variables.states_output: + variables.output_writer_0D_states.append({"format": "Paraview", + "outputInterval": 1, + "filename": "out/" + subfolder + variables.scenario_name + "/0D_states", + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"}) + +if variables.adios_output: + if variables.paraview_output: + subfolder = "adios/" + variables.output_writer_emg.append({"format": "MegaMol", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D_emg), + "filename": "out/" + subfolder + variables.scenario_name + "/hd_emg", + "useFrontBackBuffer": False, + "combineNInstances": 1, + "fileNumbering": "incremental"}) + variables.output_writer_elasticity.append({"format": "MegaMol", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D), + "filename": "out/" + subfolder + variables.scenario_name + "/elasticity", + "useFrontBackBuffer": False, + "fileNumbering": "incremental"}) + variables.output_writer_fibers.append({"format": "MegaMol", + "outputInterval": int(1. / variables.dt_splitting * variables.output_timestep_fibers), + "filename": "out/" + subfolder + variables.scenario_name + "/fibers", + "combineNInstances": variables.n_subdomains_xy, + "useFrontBackBuffer": False, + "fileNumbering": "incremental"}) + # variables.output_writer_fibers.append({"format": "MegaMol", + # "outputInterval": + # int(1./variables.dt_splitting*variables.output_timestep_fibers), + # "filename": "out/" + variables.scenario_name + "/fibers", + # "combineNInstances": 1, "useFrontBackBuffer": False, "fileNumbering": + # "incremental"} + +if variables.python_output: + if variables.adios_output: + subfolder = "python/" + variables.output_writer_emg.append({"format": "PythonFile", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D_emg), + "filename": "out/" + subfolder + variables.scenario_name + "/hd_emg", + "binary": True, + "fileNumbering": "incremental"}) + variables.output_writer_elasticity.append({"format": "PythonFile", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D), + "filename": "out/" + subfolder + variables.scenario_name + "/elasticity", + "binary": True, + "fileNumbering": "incremental"}) + variables.output_writer_fibers.append({"format": "PythonFile", + "outputInterval": int(1. / variables.dt_splitting * variables.output_timestep_fibers), + "filename": "out/" + subfolder + variables.scenario_name + "/fibers", + "binary": True, + "fileNumbering": "incremental"}) + +if variables.exfile_output: + if variables.adios_output: + subfolder = "exfile/" + variables.output_writer_emg.append({"format": "Exfile", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D_emg), + "filename": "out/" + subfolder + variables.scenario_name + "/hd_emg", + "fileNumbering": "incremental"}) + variables.output_writer_elasticity.append({"format": "Exfile", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D), + "filename": "out/" + subfolder + variables.scenario_name + "/elasticity", + "fileNumbering": "incremental"}) + variables.output_writer_fibers.append({"format": "Exfile", + "outputInterval": int(1. / variables.dt_splitting * variables.output_timestep_fibers), + "filename": "out/" + subfolder + variables.scenario_name + "/fibers", + "fileNumbering": "incremental"}) + +# set variable mappings for cellml model +if "hodgkin-huxley_shorten_ocallaghan_davidson_soboleva_2007.cellml" in variables.cellml_file: # hodgkin huxley membrane with fatigue from shorten + # parameters: I_stim, fiber stretch λ, fiber contraction velocity \dot{λ} + variables.mappings = { + ("parameter", 0): ("constant", "membrane/i_Stim"), # parameter 0 is I_stim + ("parameter", 1): ("constant", "razumova/l_hs"), # parameter 1 is fiber stretch λ + ("parameter", 2): ("constant", "razumova/velocity"), # parameter 2 is fiber contraction velocity \dot{λ} + ("connectorSlot", "vm"): ("state", "membrane/V"), # expose state Vm to the operator splitting + ("connectorSlot", "stress"): ("algebraic", "razumova/activestress"), # expose algebraic γ to the operator splitting + ("connectorSlot", "alpha"): ("algebraic", "razumova/activation"), # expose algebraic α to the operator splitting + + # connect output "lamda" of mechanics solver to parameter 1 (l_hs) + ("connectorSlot", "lambda"): "razumova/l_hs", + # connect output "ldot" of mechanics solver to parameter 2 (rel_velo) + ("connectorSlot", "ldot"): "razumova/velocity", + } + # Aliev_Panfilov/I_HH = I_stim, Razumova/l_hs = λ, Razumova/rel_velo = \dot{λ} + variables.parameters_initial_values = [0, 1, 0] + variables.nodal_stimulation_current = 40. # not used + # to which value of Vm the stimulated node should be set (option "valueForStimulatedPoint" of FastMonodomainSolver) + variables.vm_value_stimulated = 40. + # disable computation of force-length relation in opendihu, as it is carried out in CellML model + variables.enable_force_length_relation = False + # scaling factor to convert dimensionless contraction velocity to + # shortening velocity, velocity = factor*\dot{lambda} + variables.lambda_dot_scaling_factor = 7.815e-05 + +elif "hodgkin_huxley" in variables.cellml_file: + # parameters: I_stim + variables.mappings = { + ("parameter", 0): ("constant", "membrane/i_Stim"), # parameter 0 is constant 2 = I_stim + ("connectorSlot", "vm"): "membrane/V", # expose state 0 = Vm to the operator splitting + } + variables.parameters_initial_values = [0.0] # initial value for stimulation current + variables.nodal_stimulation_current = 40. # not used + # to which value of Vm the stimulated node should be set (option "valueForStimulatedPoint" of FastMonodomainSolver) + variables.vm_value_stimulated = 20. + +elif "shorten" in variables.cellml_file: + # parameters: stimulation current I_stim, fiber stretch λ + variables.mappings = { + ("parameter", 0): ("algebraic", "wal_environment/I_HH"), # parameter is algebraic 32 + # parameter is constant 65, fiber stretch λ, this indicates how much the + # fiber has stretched, 1 means no extension + ("parameter", 1): ("constant", "razumova/L_x"), + ("connectorSlot", "vm"): "wal_environment/vS", # expose state 0 = Vm to the operator splitting + } + variables.parameters_initial_values = [0.0, 1.0] # stimulation current I_stim, fiber stretch λ + variables.nodal_stimulation_current = 1200. # not used + # to which value of Vm the stimulated node should be set (option "valueForStimulatedPoint" of FastMonodomainSolver) + variables.vm_value_stimulated = 40. + +elif "slow_TK_2014" in variables.cellml_file: # this is (3a, "MultiPhysStrain", old tomo mechanics) in OpenCMISS + # parameters: I_stim, fiber stretch λ + variables.mappings = { + ("parameter", 0): ("constant", "wal_environment/I_HH"), # parameter 0 is constant 54 = I_stim + ("parameter", 1): ("constant", "razumova/L_S"), # parameter 1 is constant 67 = fiber stretch λ + ("connectorSlot", "vm"): "wal_environment/vS", # expose state 0 = Vm to the operator splitting + # expose algebraic 12 = γ to the operator splitting + ("connectorSlot", "stress"): "razumova/stress", + } + # wal_environment/I_HH = I_stim, razumova/L_S = λ + variables.parameters_initial_values = [0.0, 1.0] + variables.nodal_stimulation_current = 40. # not used + # to which value of Vm the stimulated node should be set (option "valueForStimulatedPoint" of FastMonodomainSolver) + variables.vm_value_stimulated = 40. + +# this is (3, "MultiPhysStrain", numerically more stable) in OpenCMISS, this only computes A1,A2,x1,x2 not the stress +elif "Aliev_Panfilov_Razumova_2016_08_22" in variables.cellml_file: + # parameters: I_stim, fiber stretch λ, fiber contraction velocity \dot{λ} + variables.mappings = { + ("parameter", 0): ("constant", "Aliev_Panfilov/I_HH"), # parameter 0 is constant 0 = I_stim + ("parameter", 1): ("constant", "Razumova/l_hs"), # parameter 1 is constant 8 = fiber stretch λ + # parameter 2 is constant 9 = fiber contraction velocity \dot{λ} + ("parameter", 2): ("constant", "Razumova/velo"), + ("connectorSlot", "vm"): "Aliev_Panfilov/V_m", # expose state 0 = Vm to the operator splitting + # expose algebraic 0 = γ to the operator splitting + ("connectorSlot", "stress"): "Razumova/sigma", + } + # Aliev_Panfilov/I_HH = I_stim, Razumova/l_hs = λ, Razumova/velo = \dot{λ} + variables.parameters_initial_values = [0, 1, 0] + variables.nodal_stimulation_current = 40. # not used + # to which value of Vm the stimulated node should be set (option "valueForStimulatedPoint" of FastMonodomainSolver) + variables.vm_value_stimulated = 40. + +elif "Aliev_Panfilov_Razumova_Titin" in variables.cellml_file: # this is (4, "Titin") in OpenCMISS + # parameters: I_stim, fiber stretch λ, fiber contraction velocity \dot{λ} + variables.mappings = { + ("parameter", 0): ("constant", "Aliev_Panfilov/I_HH"), # parameter 0 is constant 0 = I_stim + ("parameter", 1): ("constant", "Razumova/l_hs"), # parameter 1 is constant 11 = fiber stretch λ + # parameter 2 is constant 12 = fiber contraction velocity \dot{λ} + ("parameter", 2): ("constant", "Razumova/rel_velo"), + ("connectorSlot", "vm"): ("state", "Aliev_Panfilov/V_m"), # expose state 0 = Vm to the operator splitting + # expose algebraic 4 = γ to the operator splitting + ("connectorSlot", "stress"): ("algebraic", "Razumova/ActiveStress"), + # expose algebraic 5 = α to the operator splitting + ("connectorSlot", "alpha"): ("algebraic", "Razumova/Activation"), + } + # Aliev_Panfilov/I_HH = I_stim, Razumova/l_hs = λ, Razumova/rel_velo = \dot{λ} + variables.parameters_initial_values = [0, 1, 0] + variables.nodal_stimulation_current = 40. # not used + # to which value of Vm the stimulated node should be set (option "valueForStimulatedPoint" of FastMonodomainSolver) + variables.vm_value_stimulated = 40. + + +# callback functions +# -------------------------- +def get_motor_unit_no(fiber_no): + return int(variables.fiber_distribution[fiber_no % len(variables.fiber_distribution)] - 1) + + +def get_diffusion_prefactor(fiber_no, mu_no): + diffusion_prefactor = variables.get_conductivity( + fiber_no, mu_no) / (variables.get_am(fiber_no, mu_no) * variables.get_cm(fiber_no, mu_no)) + # print("diffusion_prefactor: {}/({}*{}) = {}".format(variables.get_conductivity(fiber_no, mu_no), variables.get_am(fiber_no, mu_no), variables.get_cm(fiber_no, mu_no), diffusion_prefactor)) + return diffusion_prefactor + + +def fiber_gets_stimulated(fiber_no, frequency, current_time): + """ + determine if fiber fiber_no gets stimulated at simulation time current_time + """ + + # determine motor unit + alpha = 1.0 # 0.8 + mu_no = (int)(get_motor_unit_no(fiber_no) * alpha) + + # determine if fiber fires now + index = int(np.round(current_time * frequency)) + n_firing_times = np.size(variables.firing_times, 0) + + # if variables.firing_times[index % n_firing_times, mu_no] == 1: + # print("{}: fiber {} is mu {}, t = {}, row: {}, stimulated: {} {}".format(rank_no, fiber_no, mu_no, current_time, (index % n_firing_times), variables.firing_times[index % n_firing_times, mu_no], "true" if variables.firing_times[index % n_firing_times, mu_no] == 1 else "false")) + + return variables.firing_times[index % n_firing_times, mu_no] == 1 + + +def set_parameters(n_nodes_global, time_step_no, current_time, parameters, dof_nos_global, fiber_no): + + # determine if fiber gets stimulated at the current time + is_fiber_gets_stimulated = fiber_gets_stimulated(fiber_no, variables.stimulation_frequency, current_time) + + # determine nodes to stimulate (center node, left and right neighbour) + innervation_zone_width_n_nodes = variables.innervation_zone_width * 100 # 100 nodes per cm + # + np.random.randint(-innervation_zone_width_n_nodes/2,innervation_zone_width_n_nodes/2+1) + innervation_node_global = int(n_nodes_global / 2) + nodes_to_stimulate_global = [innervation_node_global] + if innervation_node_global > 0: + nodes_to_stimulate_global.insert(0, innervation_node_global - 1) + if innervation_node_global < n_nodes_global - 1: + nodes_to_stimulate_global.append(innervation_node_global + 1) + + # stimulation value + if is_fiber_gets_stimulated: + stimulation_current = variables.nodal_stimulation_current + else: + stimulation_current = 0. + + first_dof_global = dof_nos_global[0] + last_dof_global = dof_nos_global[-1] + + for node_no_global in nodes_to_stimulate_global: + if first_dof_global <= node_no_global <= last_dof_global: + # get local no for global no (1D) + dof_no_local = node_no_global - first_dof_global + parameters[dof_no_local] = stimulation_current + +# callback function that can set parameters, i.e. stimulation current + + +def set_specific_parameters(n_nodes_global, time_step_no, current_time, parameters, fiber_no): + + # determine if fiber gets stimulated at the current time + is_fiber_gets_stimulated = fiber_gets_stimulated(fiber_no, variables.stimulation_frequency, current_time) + + # determine nodes to stimulate (center node, left and right neighbour) + innervation_zone_width_n_nodes = variables.innervation_zone_width * 100 # 100 nodes per cm + # + np.random.randint(-innervation_zone_width_n_nodes/2,innervation_zone_width_n_nodes/2+1) + innervation_node_global = int(n_nodes_global / 2) + nodes_to_stimulate_global = [innervation_node_global] + + for k in range(10): + if innervation_node_global - k >= 0: + nodes_to_stimulate_global.insert(0, innervation_node_global - k) + if innervation_node_global + k <= n_nodes_global - 1: + nodes_to_stimulate_global.append(innervation_node_global + k) + + # stimulation value + if is_fiber_gets_stimulated: + stimulation_current = 40. + else: + stimulation_current = 0. + + for node_no_global in nodes_to_stimulate_global: + parameters[(node_no_global, 0)] = stimulation_current # key: ((x,y,z),nodal_dof_index) + +# callback function that can set states, i.e. prescribed values for stimulation + + +def set_specific_states(n_nodes_global, time_step_no, current_time, states, fiber_no): + + # print("call set_specific_states at time {}".format(current_time)) + + # determine if fiber gets stimulated at the current time + is_fiber_gets_stimulated = fiber_gets_stimulated(fiber_no, variables.stimulation_frequency, current_time) + + if is_fiber_gets_stimulated: + # determine nodes to stimulate (center node, left and right neighbour) + innervation_zone_width_n_nodes = variables.innervation_zone_width * 100 # 100 nodes per cm + # + np.random.randint(-innervation_zone_width_n_nodes/2,innervation_zone_width_n_nodes/2+1) + innervation_node_global = int(n_nodes_global / 2) + nodes_to_stimulate_global = [innervation_node_global] + if innervation_node_global > 0: + nodes_to_stimulate_global.insert(0, innervation_node_global - 1) + if innervation_node_global < n_nodes_global - 1: + nodes_to_stimulate_global.append(innervation_node_global + 1) + # if rank_no == 0: + # print("t: {}, stimulate fiber {} at nodes {}".format(current_time, fiber_no, nodes_to_stimulate_global)) + + for node_no_global in nodes_to_stimulate_global: + states[(node_no_global, 0, 0)] = variables.vm_value_stimulated # key: ((x,y,z),nodal_dof_index,state_no) + +# callback function for artifical stress values, instead of monodomain + + +def set_stress_values(n_dofs_global, n_nodes_global_per_coordinate_direction, + time_step_no, current_time, values, global_natural_dofs, fiber_no): + # n_dofs_global: (int) global number of dofs in the mesh where to set the values + # n_nodes_global_per_coordinate_direction (list of ints) [mx, my, mz] number of global nodes in each coordinate direction. + # For composite meshes, the values are only for the first submesh, for other meshes sum(...) equals n_dofs_global + # time_step_no: (int) current time step number + # current_time: (float) the current simulation time + # values: (list of floats) all current local values of the field variable, if there are multiple components, they are stored in struct-of-array memory layout + # i.e. [point0_component0, point0_component1, ... pointN_component0, point1_component0, point1_component1, ...] + # After the call, these values will be assigned to the field variable. + # global_natural_dofs (list of ints) for every local dof no. the dof no. in global natural ordering + # additional_argument: The value of the option "additionalArgument", can be any Python object. + + # loop over nodes in fiber + for local_dof_no in range(len(values)): + # get the global no. of the current dof + global_dof_no = global_natural_dofs[local_dof_no] + + n_nodes_per_fiber = n_nodes_global_per_coordinate_direction[0] + + k = global_dof_no + N = n_nodes_per_fiber + + if k > N / 2: + k = N / 2 - k + else: + k = k - N / 2 + + values[local_dof_no] = 0.1 * np.sin((current_time / 100 + 0.2 * k / N + 0.1 * + fiber_no / variables.n_fibers_total) * 2 * np.pi) ** 2 + + +# load MU distribution and firing times +variables.fiber_distribution = np.genfromtxt(variables.fiber_distribution_file, delimiter=" ") +variables.firing_times = np.genfromtxt(variables.firing_times_file) + +# for debugging output show when the first 20 fibers will fire +if rank_no == 0 and not variables.disable_firing_output: + print("Debugging output about fiber firing: Taking input from file \"{}\"".format(variables.firing_times_file)) + import timeit + t_start = timeit.default_timer() + + first_stimulation_info = [] + + n_firing_times = np.size(variables.firing_times, 0) + for fiber_no_index in range(variables.n_fibers_total): + if fiber_no_index % 100 == 0: + t_algebraic = timeit.default_timer() + if t_algebraic - t_start > 100: + print("Note: break after {}/{} fibers ({:.0f}%) because it already took {:.3f}s".format(fiber_no_index, + variables.n_fibers_total, 100.0 * fiber_no_index / (variables.n_fibers_total - 1.), t_algebraic - t_start)) + break + + first_stimulation = None + for current_time in np.linspace(0, 1. / variables.stimulation_frequency * n_firing_times, n_firing_times): + if fiber_gets_stimulated(fiber_no_index, variables.stimulation_frequency, current_time): + first_stimulation = current_time + break + mu_no = get_motor_unit_no(fiber_no_index) + first_stimulation_info.append([fiber_no_index, mu_no, first_stimulation]) + + first_stimulation_info.sort( + key=lambda x: 1e6 + + 1e-6 * + x[1] + + 1e-12 * + x[0] if x[2] is None else x[2] + + 1e-6 * + x[1] + + 1e-12 * + x[0]) + + print("First stimulation times") + print(" Time MU fibers") + n_stimulated_mus = 0 + n_not_stimulated_mus = 0 + stimulated_fibers = [] + last_time = 0 + last_mu_no = first_stimulation_info[0][1] + for stimulation_info in first_stimulation_info: + mu_no = stimulation_info[1] + fiber_no = stimulation_info[0] + if mu_no == last_mu_no: + stimulated_fibers.append(fiber_no) + else: + if last_time is not None: + if len(stimulated_fibers) > 10: + print("{:8.2f} {:3} {} (only showing first 10, {} total)".format( + last_time, last_mu_no, str(stimulated_fibers[0:10]), len(stimulated_fibers))) + else: + print("{:8.2f} {:3} {}".format(last_time, last_mu_no, str(stimulated_fibers))) + n_stimulated_mus += 1 + else: + if len(stimulated_fibers) > 10: + print(" never stimulated: MU {:3}, fibers {} (only showing first 10, {} total)".format( + last_mu_no, str(stimulated_fibers[0:10]), len(stimulated_fibers))) + else: + print(" never stimulated: MU {:3}, fibers {}".format(last_mu_no, str(stimulated_fibers))) + n_not_stimulated_mus += 1 + stimulated_fibers = [fiber_no] + + last_time = stimulation_info[2] + last_mu_no = mu_no + + print("stimulated MUs: {}, not stimulated MUs: {}".format(n_stimulated_mus, n_not_stimulated_mus)) + + t_end = timeit.default_timer() + print("duration of assembling this list: {:.3f} s\n".format(t_end - t_start)) + +# compute partitioning +if rank_no == 0: + if n_ranks != variables.n_subdomains_x * variables.n_subdomains_y * variables.n_subdomains_z: + print( + "\n\nError! Number of ranks {} does not match given partitioning {} x {} x {} = {}.\n\n".format( + n_ranks, + variables.n_subdomains_x, + variables.n_subdomains_y, + variables.n_subdomains_z, + variables.n_subdomains_x * + variables.n_subdomains_y * + variables.n_subdomains_z)) + quit() + +# n_fibers_per_subdomain_* is already set + +#################################### +# set Dirichlet BC for the flow problem + +n_points_3D_mesh_linear_global_x = sum([n_sampled_points_in_subdomain_x(subdomain_coordinate_x) + for subdomain_coordinate_x in range(variables.n_subdomains_x)]) +n_points_3D_mesh_linear_global_y = sum([n_sampled_points_in_subdomain_y(subdomain_coordinate_y) + for subdomain_coordinate_y in range(variables.n_subdomains_y)]) +n_points_3D_mesh_linear_global_z = sum([n_sampled_points_in_subdomain_z(subdomain_coordinate_z) + for subdomain_coordinate_z in range(variables.n_subdomains_z)]) +n_points_3D_mesh_linear_global = n_points_3D_mesh_linear_global_x * \ + n_points_3D_mesh_linear_global_y * n_points_3D_mesh_linear_global_z + +n_points_3D_mesh_quadratic_global_x = 2 * n_points_3D_mesh_linear_global_x - 1 +n_points_3D_mesh_quadratic_global_y = 2 * n_points_3D_mesh_linear_global_y - 1 +n_points_3D_mesh_quadratic_global_z = 2 * n_points_3D_mesh_linear_global_z - 1 + +# set boundary conditions for the elasticity +[mx, my, mz] = variables.meshes["3Dmesh_quadratic"]["nPointsGlobal"] +[nx, ny, nz] = variables.meshes["3Dmesh_quadratic"]["nElements"] + +variables.fiber_mesh_names = [mesh_name for mesh_name in variables.meshes.keys() if "MeshFiber" in mesh_name] + +# set Dirichlet BC at top nodes for linear elasticity problem, fix muscle at top +variables.elasticity_dirichlet_bc = {} +if False: + for j in range(my): + for i in range(mx): + variables.elasticity_dirichlet_bc[(mz - 1) * mx * my + j * mx + i] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + +# fix muscle at bottom +if False: + k = 0 + for j in range(my): + for i in range(mx): + variables.elasticity_dirichlet_bc[k * mx * my + j * mx + i] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + +# Neumann BC at top nodes, traction upwards +# k = nz-1 +# variables.elasticity_neumann_bc = [{"element": k*nx*ny + j*nx + i, "constantVector": [0.0,0.0,10.0], "face": "2+"} for j in range(ny) for i in range(nx)] +variables.elasticity_neumann_bc = [] + +# with open("mesh","w") as f: +# f.write(str(variables.meshes["3Dmesh_quadratic"])) diff --git a/muscle-tendon-complex/muscle-opendihu/run.sh b/muscle-tendon-complex/muscle-opendihu/run.sh new file mode 100755 index 000000000..d09eff9c4 --- /dev/null +++ b/muscle-tendon-complex/muscle-opendihu/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +./muscle-solver settings-muscle.py variables.py diff --git a/muscle-tendon-complex/muscle-opendihu/settings-muscle.py b/muscle-tendon-complex/muscle-opendihu/settings-muscle.py new file mode 100644 index 000000000..e5675cb94 --- /dev/null +++ b/muscle-tendon-complex/muscle-opendihu/settings-muscle.py @@ -0,0 +1,882 @@ +# Multiple 1D fibers (monodomain) with 3D dynamic mooney rivlin with active contraction term, on biceps geometry + +import sys +import os +import timeit +import argparse +import importlib +import distutils.util + +# set title of terminal +title = "muscle" +print('\33]0;{}\a'.format(title), end='', flush=True) + +# parse rank arguments +rank_no = (int)(sys.argv[-2]) +n_ranks = (int)(sys.argv[-1]) + +# add variables subfolder to python path where the variables script is located +script_path = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, script_path) +sys.path.insert(0, os.path.join(script_path, 'variables')) + +import variables # file variables.py, defined default values for all parameters, you can set the parameters there +# file create_partitioned_meshes_for_settings with helper functions about own subdomain +from create_partitioned_meshes_for_settings import * + +# if first argument contains "*.py", it is a custom variable definition file, load these values +if ".py" in sys.argv[0]: + variables_path_and_filename = sys.argv[0] + variables_path, variables_filename = os.path.split(variables_path_and_filename) # get path and filename + # add the directory of the variables file to python path + sys.path.insert(0, os.path.join(script_path, variables_path)) + # remove the ".py" extension to get the name of the module + variables_module, _ = os.path.splitext(variables_filename) + + if rank_no == 0: + print("Loading variables from \"{}\".".format(variables_path_and_filename)) + + custom_variables = importlib.import_module(variables_module, + package=variables_filename) # import variables module + variables.__dict__.update(custom_variables.__dict__) + sys.argv = sys.argv[1:] # remove first argument, which now has already been parsed +else: + if rank_no == 0: + print("Error: no variables file was specified, e.g:\n ./biceps_contraction ../settings_biceps_contraction.py ramp.py") + exit(0) + +# -------------- begin user parameters ---------------- +variables.output_timestep_3D = 50 # [ms] output timestep of mechanics +variables.output_timestep_fibers = 50 # [ms] output timestep of fibers +# -------------- end user parameters ---------------- + +# define command line arguments + + +def mbool(x): return bool(distutils.util.strtobool(x)) # function to parse bool arguments + + +parser = argparse.ArgumentParser(description='muscle') +parser.add_argument( + '--scenario_name', + help='The name to identify this run in the log.', + default=variables.scenario_name) +parser.add_argument('--n_subdomains', nargs=3, help='Number of subdomains in x,y,z direction.', type=int) +parser.add_argument( + '--n_subdomains_x', + '-x', + help='Number of subdomains in x direction.', + type=int, + default=variables.n_subdomains_x) +parser.add_argument( + '--n_subdomains_y', + '-y', + help='Number of subdomains in y direction.', + type=int, + default=variables.n_subdomains_y) +parser.add_argument( + '--n_subdomains_z', + '-z', + help='Number of subdomains in z direction.', + type=int, + default=variables.n_subdomains_z) +parser.add_argument( + '--diffusion_solver_type', + help='The solver for the diffusion.', + default=variables.diffusion_solver_type, + choices=[ + "gmres", + "cg", + "lu", + "gamg", + "richardson", + "chebyshev", + "cholesky", + "jacobi", + "sor", + "preonly"]) +parser.add_argument( + '--diffusion_preconditioner_type', + help='The preconditioner for the diffusion.', + default=variables.diffusion_preconditioner_type, + choices=[ + "jacobi", + "sor", + "lu", + "ilu", + "gamg", + "none"]) +parser.add_argument( + '--potential_flow_solver_type', + help='The solver for the potential flow (non-spd matrix).', + default=variables.potential_flow_solver_type, + choices=[ + "gmres", + "cg", + "lu", + "gamg", + "richardson", + "chebyshev", + "cholesky", + "jacobi", + "sor", + "preonly"]) +parser.add_argument( + '--potential_flow_preconditioner_type', + help='The preconditioner for the potential flow.', + default=variables.potential_flow_preconditioner_type, + choices=[ + "jacobi", + "sor", + "lu", + "ilu", + "gamg", + "none"]) +parser.add_argument( + '--paraview_output', + help='Enable the paraview output writer.', + default=variables.paraview_output, + action='store_true') +parser.add_argument( + '--adios_output', + help='Enable the MegaMol/ADIOS output writer.', + default=variables.adios_output, + action='store_true') +parser.add_argument( + '--fiber_file', + help='The filename of the file that contains the fiber data.', + default=variables.fiber_file) +parser.add_argument( + '--fiber_distribution_file', + help='The filename of the file that contains the MU firing times.', + default=variables.fiber_distribution_file) +parser.add_argument( + '--firing_times_file', + help='The filename of the file that contains the cellml model.', + default=variables.firing_times_file) +parser.add_argument( + '--end_time', + '--tend', + '-t', + help='The end simulation time.', + type=float, + default=variables.end_time) +parser.add_argument( + '--output_timestep', + help='The timestep for writing outputs.', + type=float, + default=variables.output_timestep) +parser.add_argument('--dt_0D', help='The timestep for the 0D model.', type=float, default=variables.dt_0D) +parser.add_argument('--dt_1D', help='The timestep for the 1D model.', type=float, default=variables.dt_1D) +parser.add_argument( + '--dt_splitting', + help='The timestep for the splitting.', + type=float, + default=variables.dt_splitting) +parser.add_argument( + '--dt_3D', + help='The timestep for the 3D model, i.e. dynamic solid mechanics.', + type=float, + default=variables.dt_3D) +parser.add_argument( + '--disable_firing_output', + help='Disables the initial list of fiber firings.', + default=variables.disable_firing_output, + action='store_true') +parser.add_argument( + '--enable_coupling', + help='Enables the precice coupling.', + type=mbool, + default=variables.enable_coupling) +parser.add_argument('--v', help='Enable full verbosity in c++ code') +parser.add_argument('-v', help='Enable verbosity level in c++ code', action="store_true") +parser.add_argument('-vmodule', help='Enable verbosity level for given file in c++ code') +parser.add_argument('-pause', help='Stop at parallel debugging barrier', action="store_true") + +# parse command line arguments and assign values to variables module +args, other_args = parser.parse_known_args(args=sys.argv[:-2], namespace=variables) +if len(other_args) != 0 and rank_no == 0: + print("Warning: These arguments were not parsed by the settings python file\n " + "\n ".join(other_args), file=sys.stderr) + + +variables.n_subdomains = variables.n_subdomains_x * variables.n_subdomains_y * variables.n_subdomains_z + +# automatically initialize partitioning if it has not been set +if n_ranks != variables.n_subdomains: + + # create all possible partitionings to the given number of ranks + optimal_value = n_ranks**(1 / 3) + possible_partitionings = [] + for i in range(1, n_ranks + 1): + for j in range(1, n_ranks + 1): + if i * j <= n_ranks and n_ranks % (i * j) == 0: + k = (int)(n_ranks / (i * j)) + performance = (k - optimal_value)**2 + (j - optimal_value)**2 + 1.1 * (i - optimal_value)**2 + possible_partitionings.append([i, j, k, performance]) + + # if no possible partitioning was found + if len(possible_partitionings) == 0: + if rank_no == 0: + print("\n\n\033[0;31mError! Number of ranks {} does not match given partitioning {} x {} x {} = {} and no automatic partitioning could be done.\n\n\033[0m".format( + n_ranks, variables.n_subdomains_x, variables.n_subdomains_y, variables.n_subdomains_z, variables.n_subdomains_x * variables.n_subdomains_y * variables.n_subdomains_z)) + quit() + + # select the partitioning with the lowest value of performance which is the best + lowest_performance = possible_partitionings[0][3] + 1 + for i in range(len(possible_partitionings)): + if possible_partitionings[i][3] < lowest_performance: + lowest_performance = possible_partitionings[i][3] + variables.n_subdomains_x = possible_partitionings[i][0] + variables.n_subdomains_y = possible_partitionings[i][1] + variables.n_subdomains_z = possible_partitionings[i][2] + +# output information of run +if rank_no == 0: + print("scenario_name: {}, n_subdomains: {} {} {}, n_ranks: {}, end_time: {}".format(variables.scenario_name, + variables.n_subdomains_x, variables.n_subdomains_y, variables.n_subdomains_z, n_ranks, variables.end_time)) + print("dt_0D: {:0.0e}, diffusion_solver_type: {}".format( + variables.dt_0D, variables.diffusion_solver_type)) + print("dt_1D: {:0.0e}, potential_flow_solver_type: {}".format( + variables.dt_1D, variables.potential_flow_solver_type)) + print("dt_splitting: {:0.0e}, emg_solver_type: {}, emg_initial_guess_nonzero: {}".format( + variables.dt_splitting, variables.emg_solver_type, variables.emg_initial_guess_nonzero)) + print("dt_3D: {:0.0e}, paraview_output: {}".format(variables.dt_3D, variables.paraview_output)) + print("output_timestep: {:0.0e} stimulation_frequency: {} 1/ms = {} Hz".format(variables.output_timestep, + variables.stimulation_frequency, variables.stimulation_frequency * 1e3)) + print("fiber_file: {}".format(variables.fiber_file)) + print("cellml_file: {}".format(variables.cellml_file)) + print("fiber_distribution_file: {}".format(variables.fiber_distribution_file)) + print("firing_times_file: {}".format(variables.firing_times_file)) + print("********************************************************************************") + + print("prefactor: sigma_eff/(Am*Cm) = {} = {} / ({}*{})".format(variables.Conductivity / + (variables.Am * variables.Cm), variables.Conductivity, variables.Am, variables.Cm)) + + # start timer to measure duration of parsing of this script + t_start_script = timeit.default_timer() + +# initialize all helper variables +from helper import * + +variables.n_subdomains_xy = variables.n_subdomains_x * variables.n_subdomains_y +variables.n_fibers_total = variables.n_fibers_x * variables.n_fibers_y + +if False: + for subdomain_coordinate_y in range(variables.n_subdomains_y): + for subdomain_coordinate_x in range(variables.n_subdomains_x): + + print( + "subdomain (x{},y{}) ranks: {} n fibers in subdomain: x{},y{}".format( + subdomain_coordinate_x, + subdomain_coordinate_y, + list( + range( + subdomain_coordinate_y * + variables.n_subdomains_x + + subdomain_coordinate_x, + n_ranks, + variables.n_subdomains_x * + variables.n_subdomains_y)), + n_fibers_in_subdomain_x(subdomain_coordinate_x), + n_fibers_in_subdomain_y(subdomain_coordinate_y))) + + for fiber_in_subdomain_coordinate_y in range(n_fibers_in_subdomain_y(subdomain_coordinate_y)): + for fiber_in_subdomain_coordinate_x in range(n_fibers_in_subdomain_x(subdomain_coordinate_x)): + print( + "({},{}) n instances: {}".format( + fiber_in_subdomain_coordinate_x, + fiber_in_subdomain_coordinate_y, + n_fibers_in_subdomain_x(subdomain_coordinate_x) * + n_fibers_in_subdomain_y(subdomain_coordinate_y))) + + +# define the config dict +config = { + "scenarioName": variables.scenario_name, + "logFormat": "csv", + # output file of a diagram that shows data connection between solvers + "solverStructureDiagramFile": "out/muscle_solver_structure.txt", + "mappingsBetweenMeshesLogFile": "out/muscle_mappings_between_meshes.txt", # log file of when mappings between meshes occur + "Meshes": variables.meshes, + "MappingsBetweenMeshes": variables.mappings_between_meshes, + "Solvers": { + "diffusionTermSolver": { # solver for the implicit timestepping scheme of the diffusion time step + "maxIterations": 1e4, + "relativeTolerance": 1e-10, + "absoluteTolerance": 1e-10, # 1e-10 absolute tolerance of the residual + "solverType": variables.diffusion_solver_type, + "preconditionerType": variables.diffusion_preconditioner_type, + "dumpFilename": "", # "out/dump_" + "dumpFormat": "matlab", + }, + "potentialFlowSolver": { # solver for the initial potential flow, that is needed to estimate fiber directions for the bidomain equation + "relativeTolerance": 1e-10, + "absoluteTolerance": 1e-10, # 1e-10 absolute tolerance of the residual + "maxIterations": 1e4, + "solverType": variables.potential_flow_solver_type, + "preconditionerType": variables.potential_flow_preconditioner_type, + "dumpFilename": "", + "dumpFormat": "matlab", + }, + "mechanicsSolver": { # solver for the dynamic mechanics problem + "relativeTolerance": 1e-10, # 1e-10 relative tolerance of the linear solver + "absoluteTolerance": 1e-10, # 1e-10 absolute tolerance of the residual of the linear solver + "solverType": "preonly", # type of the linear solver: cg groppcg pipecg pipecgrr cgne nash stcg gltr richardson chebyshev gmres tcqmr fcg pipefcg bcgs ibcgs fbcgs fbcgsr bcgsl cgs tfqmr cr pipecr lsqr preonly qcg bicg fgmres pipefgmres minres symmlq lgmres lcd gcr pipegcr pgmres dgmres tsirm cgls + "preconditionerType": "lu", # type of the preconditioner + "maxIterations": 1e4, # maximum number of iterations in the linear solver + "snesMaxFunctionEvaluations": 1e8, # maximum number of function iterations + "snesMaxIterations": 140, # maximum number of iterations in the nonlinear solver + "snesRelativeTolerance": 1e-5, # relative tolerance of the nonlinear solver + "snesAbsoluteTolerance": 1e-5, # absolute tolerance of the nonlinear solver + "snesLineSearchType": "l2", # type of linesearch, possible values: "bt" "nleqerr" "basic" "l2" "cp" "ncglinear" + "snesRebuildJacobianFrequency": 3, # how often the jacobian should be recomputed, -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time the Jacobian is built etc. -2 means rebuild at next chance but then never again + "dumpFilename": "", # dump system matrix and right hand side after every solve + "dumpFormat": "matlab", # default, ascii, matlab + } + }, + "PreciceAdapter": { # precice adapter for muscle + "timeStepOutputInterval": 100, # interval in which to display current timestep and time in console + "timestepWidth": 1, # coupling time step width, must match the value in the precice config + # if the precice coupling is enabled, if not, it simply calls the nested solver, for debugging + "couplingEnabled": variables.enable_coupling, + "preciceConfigFilename": variables.precice_config_file, # the preCICE configuration file + # name of the own precice participant, has to match the name given in the precice xml config file + "preciceParticipantName": "Muscle", + "scalingFactor": 1, # a factor to scale the exchanged data, prior to communication + # if the output writers should be called only after a time window of + # precice is complete, this means the timestep has converged + "outputOnlyConvergedTimeSteps": True, + "preciceMeshes": [ # the precice meshes get created as the top or bottom surface of the main geometry mesh of the nested solver + { + "preciceMeshName": "Muscle-Bottom-Mesh", # precice name of the 2D coupling mesh + "face": "2-", # face of the 3D mesh where the 2D mesh is located, "2-" = bottom, "2+" = top + }, + { + "preciceMeshName": "Muscle-Top-A-Mesh", # precice name of the 2D coupling mesh + "face": "2+", # face of the 3D mesh where the 2D mesh is located, "2-" = bottom, "2+" = top + }, + { + "preciceMeshName": "Muscle-Top-B-Mesh", # precice name of the 2D coupling mesh + "face": "2+", # face of the 3D mesh where the 2D mesh is located, "2-" = bottom, "2+" = top + }, + ], + "preciceData": [ + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "read-displacements-velocities", + # name of the precice coupling surface mesh, as given in the precice xml settings file + "preciceMeshName": "Muscle-Bottom-Mesh", + # name of the displacements "data", i.e. field variable, as given in the precice xml settings file + "displacementsName": "Displacement", + # name of the velocity "data", i.e. field variable, as given in the precice xml settings file + "velocitiesName": "Velocity", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "read-displacements-velocities", + # name of the precice coupling surface mesh, as given in the precice xml settings file + "preciceMeshName": "Muscle-Top-A-Mesh", + # name of the displacements "data", i.e. field variable, as given in the precice xml settings file + "displacementsName": "Displacement", + # name of the velocity "data", i.e. field variable, as given in the precice xml settings file + "velocitiesName": "Velocity", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "read-displacements-velocities", + # name of the precice coupling surface mesh, as given in the precice xml settings file + "preciceMeshName": "Muscle-Top-B-Mesh", + # name of the displacements "data", i.e. field variable, as given in the precice xml settings file + "displacementsName": "Displacement", + # name of the velocity "data", i.e. field variable, as given in the precice xml settings file + "velocitiesName": "Velocity", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "write-traction", + # name of the precice coupling surface mesh, as given in the precice xml settings + "preciceMeshName": "Muscle-Bottom-Mesh", + # name of the traction "data", i.e. field variable, as given in the precice xml settings file + "tractionName": "Traction", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "write-traction", + # name of the precice coupling surface mesh, as given in the precice xml settings + "preciceMeshName": "Muscle-Top-A-Mesh", + # name of the traction "data", i.e. field variable, as given in the precice xml settings file + "tractionName": "Traction", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "write-traction", + # name of the precice coupling surface mesh, as given in the precice xml settings + "preciceMeshName": "Muscle-Top-B-Mesh", + # name of the traction "data", i.e. field variable, as given in the precice xml settings file + "tractionName": "Traction", + }, + ], + + "Coupling": { + "timeStepWidth": variables.dt_3D, # 1e-1 + "logTimeStepWidthAsKey": "dt_3D", + "durationLogKey": "duration_total", + "timeStepOutputInterval": 1, + "endTime": variables.end_time, + # transfer gamma to MuscleContractionSolver, the receiving slots are λ, λdot, γ + "connectedSlotsTerm1To2": {1: 2}, + "connectedSlotsTerm2To1": None, # transfer nothing back + "Term1": { # monodomain, fibers + "MultipleInstances": { + "logKey": "duration_subdomains_xy", + "ranksAllComputedInstances": list(range(n_ranks)), + "nInstances": variables.n_subdomains_xy, + "instances": + [{ + "ranks": list(range(subdomain_coordinate_y * variables.n_subdomains_x + subdomain_coordinate_x, n_ranks, variables.n_subdomains_x * variables.n_subdomains_y)), + + # this is for the actual model with fibers + "StrangSplitting": { + # "numberTimeSteps": 1, + "timeStepWidth": variables.dt_splitting, # 1e-1 + "logTimeStepWidthAsKey": "dt_splitting", + "durationLogKey": "duration_monodomain", + "timeStepOutputInterval": 100, + "endTime": variables.dt_splitting, + # transfer slot 0 = state Vm from Term1 (CellML) to Term2 (Diffusion) + "connectedSlotsTerm1To2": [0, 1, 2], + "connectedSlotsTerm2To1": [0, None, 2], # transfer the same back, this avoids data copy + + "Term1": { # CellML, i.e. reaction term of Monodomain equation + "MultipleInstances": { + "logKey": "duration_subdomains_z", + "nInstances": n_fibers_in_subdomain_x(subdomain_coordinate_x) * n_fibers_in_subdomain_y(subdomain_coordinate_y), + "instances": + [{ + # these rank nos are local nos to the outer instance of MultipleInstances, + # i.e. from 0 to number of ranks in z direction + "ranks": list(range(variables.n_subdomains_z)), + "Heun": { + "timeStepWidth": variables.dt_0D, # timestep width of 0D problem + # key under which the time step width will be written to the log file + "logTimeStepWidthAsKey": "dt_0D", + "durationLogKey": "duration_0D", # log key of duration for this solver + "timeStepOutputInterval": 1e4, # how often to print the current timestep + "initialValues": [], # no initial values are specified + "dirichletBoundaryConditions": {}, # no Dirichlet boundary conditions are specified + # filename for a vtp file that contains the Dirichlet boundary condition + # nodes and their values, set to None to disable + "dirichletOutputFilename": None, + # the boundary conditions and initial values would be given as global + # numbers + "inputMeshIsGlobal": True, + "checkForNanInf": True, # abort execution if the solution contains nan or inf values + "nAdditionalFieldVariables": 0, # number of additional field variables + "additionalSlotNames": "", + + "CellML": { + "modelFilename": variables.cellml_file, # input C++ source file or cellml XML file + # "statesInitialValues": [], # if given, the initial values for the the states of one instance + "statesInitialValues": variables.states_initial_values, # initial values for new_slow_TK + # if the equilibrium values of the states should be computed before the + # simulation starts + "initializeStatesToEquilibrium": False, + # if initializeStatesToEquilibrium is enable, the timestep width to use to + # solve the equilibrium equation + "initializeStatesToEquilibriumTimestepWidth": 1e-4, + + # optimization parameters + # "vc", "simd", "openmp" type of generated optimizated source file + "optimizationType": "vc", + # if optimizationType is "vc", whether the exponential function exp(x) + # should be approximate by (1+x/n)^n with n=1024 + "approximateExponentialFunction": True, + # compiler flags used to compile the optimized model code + "compilerFlags": "-fPIC -O3 -march=native -Wno-deprecated-declarations -shared ", + # if optimizationType is "openmp", the maximum number of threads to use. + # Default value 0 means no restriction. + "maximumNumberOfThreads": 0, + + # stimulation callbacks + # "libraryFilename": "cellml_simd_lib.so", # compiled library + # "setSpecificParametersFunction": set_specific_parameters, # callback function that sets parameters like stimulation current + # "setSpecificParametersCallInterval": int(1./variables.stimulation_frequency/variables.dt_0D), # set_specific_parameters should be called every 0.1, 5e-5 * 1e3 = 5e-2 = 0.05 + # callback function that sets states like Vm, activation can be + # implemented by using this method and directly setting Vm values, or by + # using setParameters/setSpecificParameters + "setSpecificStatesFunction": set_specific_states, + # "setSpecificStatesCallInterval": 2*int(1./variables.stimulation_frequency/variables.dt_0D), # set_specific_states should be called variables.stimulation_frequency times per ms, the factor 2 is needed because every Heun step includes two calls to rhs + "setSpecificStatesCallInterval": 0, # 0 means disabled + # set_specific_states should be called variables.stimulation_frequency + # times per ms + "setSpecificStatesCallFrequency": variables.get_specific_states_call_frequency(fiber_no, motor_unit_no), + # random value to add or substract to setSpecificStatesCallFrequency every + # stimulation, this is to add random jitter to the frequency + "setSpecificStatesFrequencyJitter": variables.get_specific_states_frequency_jitter(fiber_no, motor_unit_no), + # [ms] simulation time span for which the setSpecificStates callback will be called after a call was triggered + "setSpecificStatesRepeatAfterFirstCall": 0.01, + # [ms] first time when to call setSpecificStates + "setSpecificStatesCallEnableBegin": variables.get_specific_states_call_enable_begin(fiber_no, motor_unit_no), + # last argument that will be passed to the callback functions + # set_specific_states, set_specific_parameters, etc. + "additionalArgument": fiber_no, + + # parameters to the cellml model + # mappings between parameters and algebraics/constants and between + # outputConnectorSlots and states, algebraics or parameters, they are + # defined in helper.py + "mappings": variables.mappings, + # [0.0, 1.0], # initial values for the parameters: I_Stim, l_hs + "parametersInitialValues": variables.parameters_initial_values, + + # reference to the fiber mesh + "meshName": "MeshFiber_{}".format(fiber_no), + # a file that will contain the times of stimulations + "stimulationLogFilename": "out/stimulation.log", + }, + "OutputWriter": [ + {"format": "Paraview", + "outputInterval": 1, + "filename": "out/" + variables.scenario_name + "/0D_states({},{})".format(fiber_in_subdomain_coordinate_x, + fiber_in_subdomain_coordinate_y), + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"} + ] if variables.states_output else [] + + }, + } for fiber_in_subdomain_coordinate_y in range(n_fibers_in_subdomain_y(subdomain_coordinate_y)) \ + for fiber_in_subdomain_coordinate_x in range(n_fibers_in_subdomain_x(subdomain_coordinate_x)) \ + for fiber_no in [get_fiber_no(subdomain_coordinate_x, subdomain_coordinate_y, fiber_in_subdomain_coordinate_x, fiber_in_subdomain_coordinate_y)] \ + for motor_unit_no in [get_motor_unit_no(fiber_no)]], + } + }, + "Term2": { # Diffusion + "MultipleInstances": { + "nInstances": n_fibers_in_subdomain_x(subdomain_coordinate_x) * n_fibers_in_subdomain_y(subdomain_coordinate_y), + "instances": + [{ + # these rank nos are local nos to the outer instance of MultipleInstances, + # i.e. from 0 to number of ranks in z direction + "ranks": list(range(variables.n_subdomains_z)), + "CrankNicolson": { + "initialValues": [], # no initial values are given + # "numberTimeSteps": 1, + "timeStepWidth": variables.dt_1D, # timestep width for the diffusion problem + "timeStepWidthRelativeTolerance": 1e-10, + # key under which the time step width will be written to the log file + "logTimeStepWidthAsKey": "dt_1D", + "durationLogKey": "duration_1D", # log key of duration for this solver + "timeStepOutputInterval": 1e4, # how often to print the current timestep + # old Dirichlet BC that are not used in FastMonodomainSolver: {0: + # -75.0036, -1: -75.0036}, + "dirichletBoundaryConditions": {}, + # filename for a vtp file that contains the Dirichlet boundary condition + # nodes and their values, set to None to disable + "dirichletOutputFilename": None, + "inputMeshIsGlobal": True, # initial values would be given as global numbers + "solverName": "diffusionTermSolver", # reference to the linear solver + # number of additional field variables that will be written to the output + # file, here for stress + "nAdditionalFieldVariables": 2, + "additionalSlotNames": ["stress", "activation"], + "checkForNanInf": True, # abort execution if the solution contains nan or inf values + + "FiniteElementMethod": { + "inputMeshIsGlobal": True, + "meshName": "MeshFiber_{}".format(fiber_no), + "solverName": "diffusionTermSolver", + # resolves to Conductivity / (Am * Cm) + "prefactor": get_diffusion_prefactor(fiber_no, motor_unit_no), + "slotName": None, + }, + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": int(1./variables.dt_1D*variables.output_timestep), "filename": "out/fiber_"+str(fiber_no), "binary": True, "fixedFormat": False, "combineFiles": True}, + # {"format": "Paraview", "outputInterval": 1./variables.dt_1D*variables.output_timestep, "filename": "out/fiber_"+str(i)+"_txt", "binary": False, "fixedFormat": False}, + # {"format": "ExFile", "filename": "out/fiber_"+str(i), "outputInterval": 1./variables.dt_1D*variables.output_timestep, "sphereSize": "0.02*0.02*0.02"}, + # {"format": "PythonFile", "filename": "out/fiber_"+str(i), "outputInterval": 1./variables.dt_1D*variables.output_timestep, "binary":True, "onlyNodalValues":True}, + ] + }, + } for fiber_in_subdomain_coordinate_y in range(n_fibers_in_subdomain_y(subdomain_coordinate_y)) \ + for fiber_in_subdomain_coordinate_x in range(n_fibers_in_subdomain_x(subdomain_coordinate_x)) \ + for fiber_no in [get_fiber_no(subdomain_coordinate_x, subdomain_coordinate_y, fiber_in_subdomain_coordinate_x, fiber_in_subdomain_coordinate_y)] \ + for motor_unit_no in [get_motor_unit_no(fiber_no)]], + "OutputWriter": [ + {"format": "Paraview", + "outputInterval": int(1 / variables.dt_splitting * variables.output_timestep_fibers), + "filename": "out/fibers", + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"} + ], + }, + }, + }, + + # this is for biceps_contraction_no_cell, i.e. PrescribedValues instead of fibers + "GodunovSplitting": { # this splitting scheme is only needed to replicate the solver structure as with the fibers + "timeStepWidth": variables.dt_3D, + "logTimeStepWidthAsKey": "dt_splitting", + "durationLogKey": "duration_prescribed_values", + "timeStepOutputInterval": 100, + "endTime": variables.dt_3D, + "connectedSlotsTerm1To2": [], + "connectedSlotsTerm2To1": [], # transfer the same back, this avoids data copy + + "Term1": { + "MultipleInstances": { + "logKey": "duration_subdomains_z", + "nInstances": n_fibers_in_subdomain_x(subdomain_coordinate_x) * n_fibers_in_subdomain_y(subdomain_coordinate_y), + "instances": + [{ + # these rank nos are local nos to the outer instance of MultipleInstances, + # i.e. from 0 to number of ranks in z direction + "ranks": list(range(variables.n_subdomains_z)), + "PrescribedValues": { + # reference to the fiber mesh + "meshName": "MeshFiber_{}".format(fiber_no), + "numberTimeSteps": 1, # number of timesteps to call the callback functions subsequently, this is usually 1 for prescribed values, because it is enough to set the reaction term only once per time step + "timeStepOutputInterval": 20, # if the time step should be written to console, a value > 10 produces no output + "slotNames": [], # names of the data connector slots + + # a list of field variables that will get values assigned in every + # timestep, by the provided callback function + "fieldVariables1": [ + {"name": "Vm", "callback": None}, + {"name": "stress", "callback": set_stress_values}, + ], + "fieldVariables2": [], + # a custom argument to the fieldVariables callback functions, this will be + # passed on as the last argument + "additionalArgument": fiber_no, + + "OutputWriter": [ + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_fibers), + "filename": "out/prescribed_fibers", + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"} + ] + }, + } for fiber_in_subdomain_coordinate_y in range(n_fibers_in_subdomain_y(subdomain_coordinate_y)) \ + for fiber_in_subdomain_coordinate_x in range(n_fibers_in_subdomain_x(subdomain_coordinate_x)) \ + for fiber_no in [get_fiber_no(subdomain_coordinate_x, subdomain_coordinate_y, fiber_in_subdomain_coordinate_x, fiber_in_subdomain_coordinate_y)] \ + for motor_unit_no in [get_motor_unit_no(fiber_no)]], + + # "OutputWriter" : variables.output_writer_fibers, + "OutputWriter": [ + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_fibers), + "filename": "out/fibers", + "binary": True, + "fixedFormat": False, + "combineFiles": True, + "fileNumbering": "incremental"} + ] + } + }, + + # term2 is unused, it is needed to be similar to the actual fiber solver structure + "Term2": {} + } + + } if (subdomain_coordinate_x, subdomain_coordinate_y) == (variables.own_subdomain_coordinate_x, variables.own_subdomain_coordinate_y) else None + for subdomain_coordinate_y in range(variables.n_subdomains_y) + for subdomain_coordinate_x in range(variables.n_subdomains_x)] + }, + # for FastMonodomainSolver, e.g. MU_fibre_distribution_3780.txt + "fiberDistributionFile": variables.fiber_distribution_file, + "firingTimesFile": variables.firing_times_file, # for FastMonodomainSolver, e.g. MU_firing_times_real.txt + # only compute fibers after they have been stimulated for the first time + "onlyComputeIfHasBeenStimulated": True, + # optimization where states that are close to their equilibrium will not be computed again + "disableComputationWhenStatesAreCloseToEquilibrium": True, + "valueForStimulatedPoint": variables.vm_value_stimulated, # to which value of Vm the stimulated node should be set + # range where the neuromuscular junction is located around the center, + # relative to fiber length. The actual position is draws randomly from the + # interval [0.5-s/2, 0.5+s/2) with s being this option. 0 means sharply at + # the center, 0.1 means located approximately at the center, but it can + # vary 10% in total between all fibers. + "neuromuscularJunctionRelativeSize": 0.1, + }, + "Term2": { # solid mechanics + "MuscleContractionSolver": { + "numberTimeSteps": 1, # only use 1 timestep per interval + "timeStepOutputInterval": 100, # do not output time steps + "Pmax": variables.pmax, # maximum PK2 active stress + # if the factor f_l(λ_f) modeling the force-length relation (as in + # Heidlauf2013) should be multiplied. Set to false if this relation is + # already considered in the CellML model. + "enableForceLengthRelation": variables.enable_force_length_relation, + # scaling factor for the output of the lambda dot slot, i.e. the + # contraction velocity. Use this to scale the unit-less quantity to, e.g., + # micrometers per millisecond for the subcellular model. + "lambdaDotScalingFactor": variables.lambda_dot_scaling_factor, + "slotNames": ["lambda", "ldot", "gamma", "T"], # names of the data connector slots + "OutputWriter": [ + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D), + "filename": "out/muscle_3D", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + ], + "mapGeometryToMeshes": [], # the mesh names of the meshes that will get the geometry transferred + "dynamic": True, # if the dynamic solid mechanics solver should be used, else it computes the quasi-static problem + + # the actual solid mechanics solver, this is either + # "DynamicHyperelasticitySolver" or "HyperelasticitySolver", depending on + # the value of "dynamic" + "DynamicHyperelasticitySolver": { + "timeStepWidth": variables.dt_3D, # time step width + "durationLogKey": "nonlinear", # key to find duration of this solver in the log file + "timeStepOutputInterval": 1, # how often the current time step should be printed to console + + "materialParameters": variables.material_parameters, # material parameters of the Mooney-Rivlin material + "density": variables.rho, # density of the material + # scaling factor for displacements, only set to sth. other than 1 only to + # increase visual appearance for very small displacements + "displacementsScalingFactor": 1.0, + # log file where residual norm values of the nonlinear solver will be written + "residualNormLogFilename": "muscle_log_residual_norm.txt", + # whether to use the analytically computed jacobian matrix in the nonlinear solver (fast) + "useAnalyticJacobian": True, + # whether to use the numerically computed jacobian matrix in the nonlinear + # solver (slow), only works with non-nested matrices, if both numeric and + # analytic are enable, it uses the analytic for the preconditioner and the + # numeric as normal jacobian + "useNumericJacobian": False, + + # whether to have extra output of matlab vectors, x,r, jacobian matrix (very slow) + "dumpDenseMatlabVariables": False, + # if useAnalyticJacobian,useNumericJacobian and dumpDenseMatlabVariables + # all all three true, the analytic and numeric jacobian matrices will get + # compared to see if there are programming errors for the analytic + # jacobian + + # mesh + "inputMeshIsGlobal": True, # the mesh is given locally + "meshName": "3Dmesh_quadratic", # name of the 3D mesh, it is defined under "Meshes" at the beginning of this config + # fiber meshes that will be used to determine the fiber direction, for + # multidomain there are no fibers so this would be empty list + "fiberMeshNames": variables.fiber_mesh_names, + # "fiberDirection": [0,0,1], # if fiberMeshNames is empty, directly set the constant fiber direction, in element coordinate system + + # solving + # name of the nonlinear solver configuration, it is defined under + # "Solvers" at the beginning of this config + "solverName": "mechanicsSolver", + # "loadFactors": [0.25, 0.66, 1.0], # load factors for every timestep + "loadFactors": [], # no load factors, solve problem directly + # a threshold for the load factor, when to abort the solve of the current + # time step. The load factors are adjusted automatically if the nonlinear + # solver diverged. If the progression between two subsequent load factors + # gets smaller than this value, the solution is aborted. + "loadFactorGiveUpThreshold": 4e-2, + "nNonlinearSolveCalls": 1, # how often the nonlinear solve should be repeated + + # boundary and initial conditions + # the initial Dirichlet boundary conditions that define values for + # displacements u and velocity v + "dirichletBoundaryConditions": variables.elasticity_dirichlet_bc, + # Neumann boundary conditions that define traction forces on surfaces of elements + "neumannBoundaryConditions": variables.elasticity_neumann_bc, + # if the given Neumann boundary condition values under + # "neumannBoundaryConditions" are total forces instead of surface loads + # and therefore should be scaled by the surface area of all elements where + # Neumann BC are applied + "divideNeumannBoundaryConditionValuesByTotalArea": True, + # function that updates the dirichlet BCs while the simulation is running + "updateDirichletBoundaryConditionsFunction": None, + # every which step the update function should be called, 1 means every time step + "updateDirichletBoundaryConditionsFunctionCallInterval": 1, + + # the initial values for the displacements, vector of values for every + # node [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesDisplacements": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # the initial values for the velocities, vector of values for every node + # [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesVelocities": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # if the initial values for the dynamic nonlinear problem should be + # computed by extrapolating the previous displacements and velocities + "extrapolateInitialGuess": True, + "constantBodyForce": variables.constant_body_force, # a constant force that acts on the whole body, e.g. for gravity + + # filename for a vtp file that contains the Dirichlet boundary condition + # nodes and their values, set to None to disable + "dirichletOutputFilename": "out/muscle_dirichlet_boundary_conditions", + # filename of a log file that will contain the total (bearing) forces and + # moments at the top and bottom of the volume + "totalForceLogFilename": "out/muscle_force.csv", + "totalForceLogOutputInterval": 10, # output interval when to write the totalForceLog file + # global element nos of the bottom elements used to compute the total + # forces in the log file totalForceLogFilename + "totalForceBottomElementNosGlobal": [j * nx + i for j in range(ny) for i in range(nx)], + # global element nos of the top elements used to compute the total forces + # in the log file totalForceTopElementsGlobal + "totalForceTopElementNosGlobal": [(nz - 1) * ny * nx + j * nx + i for j in range(ny) for i in range(nx)], + + + # define which file formats should be written + # 1. main output writer that writes output files using the quadratic + # elements function space. Writes displacements, velocities and PK2 + # stresses. + "OutputWriter": [ + + # Paraview files + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/u", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + + # Python callback function "postprocess" + # {"format": "PythonCallback", "outputInterval": 1, "callback": postprocess, "onlyNodalValues":True, "filename": ""}, + ], + # 2. additional output writer that writes also the hydrostatic pressure + "pressure": { # output files for pressure function space (linear elements), contains pressure values, as well as displacements and velocities + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/p", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, + # 3. additional output writer that writes virtual work terms + "dynamic": { # output of the dynamic solver, has additional virtual work values + "OutputWriter": [ # output files for displacements function space (quadratic elements) + # {"format": "Paraview", "outputInterval": int(output_interval/dt), "filename": "out/dynamic", "binary": False, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_3D * variables.output_timestep_3D), + "filename": "out/muscle_virtual_work", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + ], + }, + # 4. output writer for debugging, outputs files after each load increment, + # the geometry is not changed but u and v are written + "LoadIncrements": { + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/load_increments", "binary": False, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, + } + } + } + } + } +} + + +# stop timer and calculate how long parsing lasted +if rank_no == 0: + t_stop_script = timeit.default_timer() + print("Python config parsed in {:.1f}s.".format(t_stop_script - t_start_script)) diff --git a/muscle-tendon-complex/muscle-opendihu/variables/variables.py b/muscle-tendon-complex/muscle-opendihu/variables/variables.py new file mode 100644 index 000000000..b75aebcf6 --- /dev/null +++ b/muscle-tendon-complex/muscle-opendihu/variables/variables.py @@ -0,0 +1,170 @@ +# material parameters +# -------------------- + +c1 = 3.176e-10 # [N/cm^2] +c2 = 1.813 # [N/cm^2] +b = 1.075e-2 # [N/cm^2] anisotropy parameter +d = 9.1733 # [-] anisotropy parameter +material_parameters = [c1, c2, b, d] # material parameters +pmax = 7.3 # maximum stress [N/cm^2] +Conductivity = 3.828 # sigma, conductivity [mS/cm] +Am = 500.0 # surface area to volume ratio [cm^-1] +Cm = 0.58 # membrane capacitance [uF/cm^2] +# not used [cm], this will later be used to specify a variance of positions of the innervation point at the fibers +innervation_zone_width = 0. +rho = 10 + +# solvers +# ------- +diffusion_solver_type = "cg" # solver and preconditioner for the diffusion part of the Monodomain equation +diffusion_preconditioner_type = "none" # preconditioner +# solver and preconditioner for an initial Laplace flow on the domain, from which fiber directions are determined +potential_flow_solver_type = "gmres" +potential_flow_preconditioner_type = "none" # preconditioner +# solver and preconditioner for the 3D static Bidomain equation that solves the intra-muscular EMG signal +emg_solver_type = "cg" +emg_preconditioner_type = "none" # preconditioner +emg_initial_guess_nonzero = False # < If the initial guess for the emg linear system should be set to the previous solution + +# timing parameters +# ----------------- +end_time = 1000.0 # [ms] end time of the simulation +# [ms^-1] sampling frequency of stimuli in firing_times_file, in stimulations per ms, number before 1e-3 factor is in Hertz. +stimulation_frequency = 100 * 1e-3 +dt_0D = 1e-3 # [ms] timestep width of ODEs +dt_1D = 1.5e-3 # [ms] timestep width of diffusion +dt_splitting = 3e-3 # [ms] overall timestep width of strang splitting +# [ms] time step width of coupling, when 3D should be performed, also sampling time of monopolar EMG +dt_3D = 1e0 +output_timestep = 1e0 # [ms] timestep for output files +output_timestep_3D_emg = 1e0 # [ms] timestep for output files +output_timestep_3D = 1e0 # [ms] timestep for output files +output_timestep_fibers = 1e0 # [ms] timestep for output files +activation_start_time = 0 # [ms] time when to start checking for stimulation + +# input files +# ----------- + +import os +input_dir = os.environ.get('OPENDIHU_INPUT_DIR') +fiber_file = input_dir + "/left_biceps_brachii_9x9fibers.bin" +fat_mesh_file = fiber_file + "_fat.bin" +# use setSpecificStatesCallEnableBegin and setSpecificStatesCallFrequency +firing_times_file = input_dir + "/MU_firing_times_always.txt" +fiber_distribution_file = input_dir + "/MU_fibre_distribution_10MUs.txt" +cellml_file = input_dir + "/2020_06_03_hodgkin-huxley_shorten_ocallaghan_davidson_soboleva_2007.cellml" +firing_times_file = input_dir + "/MU_firing_times_real.txt" +precice_config_file = "../precice-config.xml" +# If the fiber geometry data should be loaded completely in the python +# script. If True, this reads the binary file and assigns the node +# positions in the config. If False, the C++ code will read the binary +# file and only extract the local node positions. This is more performant +# for highly parallel runs. +load_fiber_data = False +debug_output = False # verbose output in this python script, for debugging the domain decomposition +disable_firing_output = True # Disables the initial list of fiber firings on the console to save some console space +paraview_output = False # If the paraview output writer should be enabled +adios_output = False # If the MegaMol/ADIOS output writer should be enabled +python_output = False # If the Python output writer should be enabled +exfile_output = False # If the Exfile output writer should be enabled + +# partitioning +# ------------ +# this has to match the total number of processes +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 + +# stride for sampling the 3D elements from the fiber data +# here any number is possible +sampling_stride_x = 2 +sampling_stride_y = 2 +sampling_stride_z = 50 + +mapping_tolerance = 0.1 + +# scenario name for log file +scenario_name = "" + +# functions, here, Am, Cm and Conductivity are constant for all fibers and MU's +# These functions can be redefined differently in a custom variables script + + +def get_am(fiber_no, mu_no): + return Am + + +def get_cm(fiber_no, mu_no): + return Cm + + +def get_conductivity(fiber_no, mu_no): + return Conductivity + + +def get_specific_states_call_frequency(fiber_no, mu_no): + return stimulation_frequency + + +def get_specific_states_frequency_jitter(fiber_no, mu_no): + return [0] + + +def get_specific_states_call_enable_begin(fiber_no, mu_no): + return activation_start_time + + +# further internal variables that will be set by the helper.py script and used in the config in settings_fibers_emg.py +n_fibers_total = None +n_subdomains_xy = None +own_subdomain_coordinate_x = None +own_subdomain_coordinate_y = None +own_subdomain_coordinate_z = None +n_fibers_x = None +n_fibers_y = None +n_points_whole_fiber = None +n_points_3D_mesh_global_x = None +n_points_3D_mesh_global_y = None +n_points_3D_mesh_global_z = None +output_writer_fibers = None +output_writer_emg = None +output_writer_0D_states = None +states_output = False +parameters_used_as_algebraic = None +parameters_used_as_constant = None +parameters_initial_values = None +output_algebraic_index = None +output_state_index = None +nodal_stimulation_current = None +fiber_file_handle = None +fibers = None +fiber_distribution = None +firing_times = None +n_fibers_per_subdomain_x = None +n_fibers_per_subdomain_y = None +n_points_per_subdomain_z = None +z_point_index_start = None +z_point_index_end = None +meshes = None +potential_flow_dirichlet_bc = None +elasticity_dirichlet_bc = None +elasticity_neumann_bc = None +fibers_on_own_rank = None +n_fiber_nodes_on_subdomain = None +fiber_start_node_no = None +generate_linear_3d_mesh = False +generate_quadratic_3d_mesh = True +nx = None +ny = None +nz = None +constant_body_force = None +bottom_traction = None +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 +states_initial_values = [] +enable_coupling = True +enable_force_length_relation = True +lambda_dot_scaling_factor = 1 +mappings = None +vm_value_stimulated = None diff --git a/muscle-tendon-complex/precice-config.xml b/muscle-tendon-complex/precice-config.xml new file mode 100644 index 000000000..2c641bbb6 --- /dev/null +++ b/muscle-tendon-complex/precice-config.xml @@ -0,0 +1,181 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/muscle-tendon-complex/solver-opendihu/SConscript b/muscle-tendon-complex/solver-opendihu/SConscript new file mode 100644 index 000000000..e70812192 --- /dev/null +++ b/muscle-tendon-complex/solver-opendihu/SConscript @@ -0,0 +1,8 @@ +# This script declares to SCons how to compile the example. +# It has to be called from a SConstruct file. +# The 'env' object is passed from there and contains further specification like directory and debug/release flags. + +Import('env') + +env.Program(target = 'muscle-solver', source = "src/muscle-solver.cpp") +env.Program(target = 'tendon-solver', source = "src/tendon-solver.cpp") diff --git a/muscle-tendon-complex/solver-opendihu/SConstruct b/muscle-tendon-complex/solver-opendihu/SConstruct new file mode 100644 index 000000000..78c806c09 --- /dev/null +++ b/muscle-tendon-complex/solver-opendihu/SConstruct @@ -0,0 +1,10 @@ +import os + +# get the directory where opendihu is installed (the top level directory of opendihu) +opendihu_home = os.environ.get('OPENDIHU_HOME') +# set path where the "SConscript" file is located (set to current path) +path_where_to_call_sconscript = Dir('.').srcnode().abspath + +# call general SConstruct that will configure everything and then call SConscript at the given path +SConscript(os.path.join(opendihu_home,'SConstructGeneral'), + exports={"path": path_where_to_call_sconscript}) \ No newline at end of file diff --git a/muscle-tendon-complex/solver-opendihu/build.sh b/muscle-tendon-complex/solver-opendihu/build.sh new file mode 100755 index 000000000..6eaf7927e --- /dev/null +++ b/muscle-tendon-complex/solver-opendihu/build.sh @@ -0,0 +1,14 @@ +#!/bin/bash +set -e -u + +if [ -n "$OPENDIHU_HOME" ] +then + "$OPENDIHU_HOME/scripts/shortcuts/sr.sh" && "$OPENDIHU_HOME/scripts/shortcuts/mkorn.sh" + # copy executables to partipant folders + cp build_release/muscle-solver ../muscle-opendihu/ + cp build_release/tendon-solver ../tendon-bottom-opendihu/ + cp build_release/tendon-solver ../tendon-top-A-opendihu/ + cp build_release/tendon-solver ../tendon-top-B-opendihu/ +else + echo "OPENDIHU_HOME is not defined" +fi diff --git a/muscle-tendon-complex/solver-opendihu/clean.sh b/muscle-tendon-complex/solver-opendihu/clean.sh new file mode 100755 index 000000000..ebfabac17 --- /dev/null +++ b/muscle-tendon-complex/solver-opendihu/clean.sh @@ -0,0 +1,12 @@ +#!/bin/bash +set -e -u + +rm -rf ./*.log +rm -rf build_release +rm -rf precice-profiling + +# remove executables from partipant folders +rm -f ../muscle-opendihu/muscle-solver +rm -f ../tendon-bottom-opendihu/tendon-solver +rm -f ../tendon-top-A-opendihu/tendon-solver +rm -f ../tendon-top-B-opendihu/tendon-solver diff --git a/muscle-tendon-complex/solver-opendihu/src/muscle-solver.cpp b/muscle-tendon-complex/solver-opendihu/src/muscle-solver.cpp new file mode 100644 index 000000000..27eca51f7 --- /dev/null +++ b/muscle-tendon-complex/solver-opendihu/src/muscle-solver.cpp @@ -0,0 +1,36 @@ +#include +#include +#include +#include +#include + +int main(int argc, char *argv[]) +{ + + DihuContext settings(argc, argv); + + Control::PreciceAdapter, + BasisFunction::LagrangeOfOrder<1>>>>>, + Control::MultipleInstances< + TimeSteppingScheme::CrankNicolson< + SpatialDiscretization::FiniteElementMethod< + Mesh::StructuredDeformableOfDimension<1>, + BasisFunction::LagrangeOfOrder<1>, + Quadrature::Gauss<2>, + Equation::Dynamic::IsotropicDiffusion>>>>>>, + MuscleContractionSolver<>>> + problem(settings); + + problem.run(); + + return EXIT_SUCCESS; +} diff --git a/muscle-tendon-complex/solver-opendihu/src/tendon-solver.cpp b/muscle-tendon-complex/solver-opendihu/src/tendon-solver.cpp new file mode 100644 index 000000000..8a729edfa --- /dev/null +++ b/muscle-tendon-complex/solver-opendihu/src/tendon-solver.cpp @@ -0,0 +1,19 @@ +#include +#include +#include +#include +#include + +int main(int argc, char *argv[]) +{ + + DihuContext settings(argc, argv); + + Control::PreciceAdapter> + problem(settings); + + problem.run(); + + return EXIT_SUCCESS; +} diff --git a/muscle-tendon-complex/tendon-bottom-opendihu/clean.sh b/muscle-tendon-complex/tendon-bottom-opendihu/clean.sh new file mode 100755 index 000000000..56411170b --- /dev/null +++ b/muscle-tendon-complex/tendon-bottom-opendihu/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_opendihu . diff --git a/muscle-tendon-complex/tendon-bottom-opendihu/run.sh b/muscle-tendon-complex/tendon-bottom-opendihu/run.sh new file mode 100755 index 000000000..4de786bf9 --- /dev/null +++ b/muscle-tendon-complex/tendon-bottom-opendihu/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +./tendon-solver settings-tendon-bottom.py diff --git a/muscle-tendon-complex/tendon-bottom-opendihu/settings-tendon-bottom.py b/muscle-tendon-complex/tendon-bottom-opendihu/settings-tendon-bottom.py new file mode 100644 index 000000000..c27771d6d --- /dev/null +++ b/muscle-tendon-complex/tendon-bottom-opendihu/settings-tendon-bottom.py @@ -0,0 +1,406 @@ +# Transversely-isotropic Mooney Rivlin on a tendon geometry +# Note, this is not possible to be run in parallel because the fibers +# cannot be initialized without MultipleInstances class. +import sys +import os +import numpy as np +import argparse +import sys +from stl import mesh +import json + +# set title of terminal +title = "tendon-bottom" +print('\33]0;{}\a'.format(title), end='', flush=True) + +# add variables subfolder to python path where the variables script is located +script_path = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, script_path) +sys.path.insert(0, os.path.join(script_path, 'variables')) + +import variables +from create_partitioned_meshes_for_settings import * + +# update material parameters +if (variables.tendon_material == "nonLinear"): + c = 9.98 # [N/cm^2=kPa] + ca = 14.92 # [-] + ct = 14.7 # [-] + cat = 9.64 # [-] + ctt = 11.24 # [-] + mu = 3.76 # [N/cm^2=kPa] + k1 = 42.217e3 # [N/cm^2=kPa] + k2 = 411.360e3 # [N/cm^2=kPa] + + variables.material_parameters = [c, ca, ct, cat, ctt, mu, k1, k2] + +if (variables.tendon_material == "SaintVenantKirchoff"): + # material parameters for Saint Venant-Kirchhoff material + # https://www.researchgate.net/publication/230248067_Bulk_Modulus + + youngs_modulus = 7e4 # [N/cm^2 = 10kPa] + shear_modulus = 3e4 + + lambd = shear_modulus * (youngs_modulus - 2 * shear_modulus) / \ + (3 * shear_modulus - youngs_modulus) # Lamé parameter lambda + mu = shear_modulus # Lamé parameter mu or G (shear modulus) + + variables.material_parameters = [lambd, mu] + + +# If the fiber geometry data should be loaded completely in the python +# script. If True, this reads the binary file and assigns the node +# positions in the config. If False, the C++ code will read the binary +# file and only extract the local node positions. This is more performant +# for highly parallel runs. +load_fiber_data = False + +# parse arguments +rank_no = (int)(sys.argv[-2]) +n_ranks = (int)(sys.argv[-1]) + +# define command line arguments +parser = argparse.ArgumentParser(description='tendon') +parser.add_argument('--n_subdomains', nargs=3, help='Number of subdomains in x,y,z direction.', type=int) +parser.add_argument( + '--n_subdomains_x', + '-x', + help='Number of subdomains in x direction.', + type=int, + default=variables.n_subdomains_x) +parser.add_argument( + '--n_subdomains_y', + '-y', + help='Number of subdomains in y direction.', + type=int, + default=variables.n_subdomains_y) +parser.add_argument( + '--n_subdomains_z', + '-z', + help='Number of subdomains in z direction.', + type=int, + default=variables.n_subdomains_z) +parser.add_argument( + '--fiber_file', + help='The filename of the file that contains the fiber data.', + default=variables.fiber_file) +parser.add_argument('-vmodule', help='ignore') + +# parse command line arguments and assign values to variables module +args, other_args = parser.parse_known_args(args=sys.argv[:-2], namespace=variables) +if len(other_args) != 0 and rank_no == 0: + print("Warning: These arguments were not parsed by the settings python file\n " + "\n ".join(other_args), file=sys.stderr) + +# partitioning +# ------------ +# this has to match the total number of processes +if variables.n_subdomains is not None: + variables.n_subdomains_x = variables.n_subdomains[0] + variables.n_subdomains_y = variables.n_subdomains[1] + variables.n_subdomains_z = variables.n_subdomains[2] + +# compute partitioning +if rank_no == 0: + if n_ranks != variables.n_subdomains_x * variables.n_subdomains_y * variables.n_subdomains_z: + print( + "\n\nError! Number of ranks {} does not match given partitioning {} x {} x {} = {}.\n\n".format( + n_ranks, + variables.n_subdomains_x, + variables.n_subdomains_y, + variables.n_subdomains_z, + variables.n_subdomains_x * + variables.n_subdomains_y * + variables.n_subdomains_z)) + sys.exit(-1) + +# stride for sampling the 3D elements from the fiber data +# here any number is possible +sampling_stride_x = 1 +sampling_stride_y = 1 +sampling_stride_z = 2 + +# create the partitioning using the script in create_partitioned_meshes_for_settings.py +result = create_partitioned_meshes_for_settings( + variables.n_subdomains_x, variables.n_subdomains_y, variables.n_subdomains_z, + variables.fiber_file, load_fiber_data, + sampling_stride_x, sampling_stride_y, sampling_stride_z, True, True) + +[variables.meshes, + variables.own_subdomain_coordinate_x, + variables.own_subdomain_coordinate_y, + variables.own_subdomain_coordinate_z, + variables.n_fibers_x, + variables.n_fibers_y, + variables.n_points_whole_fiber] = result + +n_points_3D_mesh_linear_global_x = sum([n_sampled_points_in_subdomain_x(subdomain_coordinate_x) + for subdomain_coordinate_x in range(variables.n_subdomains_x)]) +n_points_3D_mesh_linear_global_y = sum([n_sampled_points_in_subdomain_y(subdomain_coordinate_y) + for subdomain_coordinate_y in range(variables.n_subdomains_y)]) +n_points_3D_mesh_linear_global_z = sum([n_sampled_points_in_subdomain_z(subdomain_coordinate_z) + for subdomain_coordinate_z in range(variables.n_subdomains_z)]) +n_points_3D_mesh_linear_global = n_points_3D_mesh_linear_global_x * \ + n_points_3D_mesh_linear_global_y * n_points_3D_mesh_linear_global_z +nx = n_points_3D_mesh_linear_global_x - 1 +ny = n_points_3D_mesh_linear_global_y - 1 +nz = n_points_3D_mesh_linear_global_z - 1 + +node_positions = variables.meshes["3Dmesh_quadratic"]["nodePositions"] + +# boundary conditions (for quadratic elements) +# -------------------------------------------- +[mx, my, mz] = variables.meshes["3Dmesh_quadratic"]["nPointsGlobal"] +[nx, ny, nz] = variables.meshes["3Dmesh_quadratic"]["nElements"] + +# set Dirichlet BC, fix x and y coordinates of the end of tendon that is attached to the bone +variables.elasticity_dirichlet_bc = {} +k = 0 + +# fix z value on the whole x-y-plane +for j in range(my): + for i in range(mx): + variables.elasticity_dirichlet_bc[k * mx * my + j * mx + i] = [0.0, 0.0, None, None, None, None] + +# set Neumann BC, set traction at the end of the tendon that is attached to the bone +k = 0 +# start with 0 BC +variables.elasticity_neumann_bc = [{"element": k * nx * ny + j * nx + i, + "constantVector": [0, 0, 0], "face": "2-"} for j in range(ny) for i in range(nx)] + + +def update_neumann_bc(t): + + # set new Neumann boundary conditions + k = 0 + factor = min(1, t / 100) # for t ∈ [0,100] from 0 to 1 + elasticity_neumann_bc = [{ + "element": k * nx * ny + j * nx + i, + "constantVector": [0, 0, -variables.force * factor], # force pointing to bottom + "face": "2-", + "isInReferenceConfiguration": True + } for j in range(ny) for i in range(nx)] + + config = { + "inputMeshIsGlobal": True, + # if the given Neumann boundary condition values under + # "neumannBoundaryConditions" are total forces instead of surface loads + # and therefore should be scaled by the surface area of all elements where + # Neumann BC are applied + "divideNeumannBoundaryConditionValuesByTotalArea": True, + "neumannBoundaryConditions": elasticity_neumann_bc, + } + print("prescribed pulling force to bottom: {}".format(variables.force * factor)) + return config + + +config_hyperelasticity = { # for both "HyperelasticitySolver" and "DynamicHyperelasticitySolver" + "timeStepWidth": variables.dt_elasticity, # time step width + "endTime": variables.end_time, # end time of the simulation time span + "durationLogKey": "duration_mechanics", # key to find duration of this solver in the log file + "timeStepOutputInterval": 1, # how often the current time step should be printed to console + + "materialParameters": variables.material_parameters, # material parameters of the Mooney-Rivlin material + "density": variables.rho, # density of the material + # scaling factor for displacements, only set to sth. other than 1 only to + # increase visual appearance for very small displacements + "displacementsScalingFactor": 1.0, + # log file where residual norm values of the nonlinear solver will be written + "residualNormLogFilename": "out/tendon_bottom_log_residual_norm.txt", + # whether to use the analytically computed jacobian matrix in the nonlinear solver (fast) + "useAnalyticJacobian": True, + # whether to use the numerically computed jacobian matrix in the nonlinear + # solver (slow), only works with non-nested matrices, if both numeric and + # analytic are enable, it uses the analytic for the preconditioner and the + # numeric as normal jacobian + "useNumericJacobian": False, + + # whether to have extra output of matlab vectors, x,r, jacobian matrix (very slow) + "dumpDenseMatlabVariables": False, + # if useAnalyticJacobian,useNumericJacobian and dumpDenseMatlabVariables + # all all three true, the analytic and numeric jacobian matrices will get + # compared to see if there are programming errors for the analytic + # jacobian + + # mesh + "meshName": "3Dmesh_quadratic", # mesh with quadratic Lagrange ansatz functions + # boundary conditions are specified in global numberings, whereas the mesh is given in local numberings + "inputMeshIsGlobal": True, + + "fiberMeshNames": [], # fiber meshes that will be used to determine the fiber direction + # "fiberDirection": [0,0,1], # if fiberMeshNames is empty, directly set the constant fiber direction, in element coordinate system + # if fiberMeshNames and fiberDirections are empty, directly set the + # constant fiber direction, in element coordinate system + "fiberDirectionInElement": [0, 0, 1], + + # nonlinear solver + "relativeTolerance": 1e-10, # 1e-10 relative tolerance of the linear solver + "absoluteTolerance": 1e-10, # 1e-10 absolute tolerance of the residual of the linear solver + "solverType": "preonly", # type of the linear solver: cg groppcg pipecg pipecgrr cgne nash stcg gltr richardson chebyshev gmres tcqmr fcg pipefcg bcgs ibcgs fbcgs fbcgsr bcgsl cgs tfqmr cr pipecr lsqr preonly qcg bicg fgmres pipefgmres minres symmlq lgmres lcd gcr pipegcr pgmres dgmres tsirm cgls + "preconditionerType": "lu", # type of the preconditioner + "maxIterations": 1e4, # maximum number of iterations in the linear solver + "snesMaxFunctionEvaluations": 1e8, # maximum number of function iterations + "snesMaxIterations": 240, # maximum number of iterations in the nonlinear solver + "snesRelativeTolerance": 1e-2, # relative tolerance of the nonlinear solver + # type of linesearch, possible values: "bt" "nleqerr" "basic" "l2" "cp" "ncglinear" + "snesLineSearchType": "l2", + "snesAbsoluteTolerance": 1e-5, # absolute tolerance of the nonlinear solver + # how often the jacobian should be recomputed, -1 indicates NEVER rebuild, + # 1 means rebuild every time the Jacobian is computed within a single + # nonlinear solve, 2 means every second time the Jacobian is built etc. -2 + # means rebuild at next chance but then never again + "snesRebuildJacobianFrequency": 5, + + # "dumpFilename": "out/r{}/m".format(sys.argv[-1]), # dump system matrix and right hand side after every solve + "dumpFilename": "", # dump disabled + "dumpFormat": "matlab", # default, ascii, matlab + + # "loadFactors": [0.1, 0.2, 0.35, 0.5, 1.0], # load factors for every timestep + # "loadFactors": [0.5, 1.0], # load factors for every timestep + "loadFactors": [], # no load factors, solve problem directly + # a threshold for the load factor, when to abort the solve of the current + # time step. The load factors are adjusted automatically if the nonlinear + # solver diverged. If the load factors get too small, it aborts the solve. + "loadFactorGiveUpThreshold": 1e-3, + "nNonlinearSolveCalls": 1, # how often the nonlinear solve should be called + + # boundary and initial conditions + # the initial Dirichlet boundary conditions that define values for displacements u and velocity v + "dirichletBoundaryConditions": variables.elasticity_dirichlet_bc, + # Neumann boundary conditions that define traction forces on surfaces of elements + "neumannBoundaryConditions": variables.elasticity_neumann_bc, + # if the given Neumann boundary condition values under + # "neumannBoundaryConditions" are total forces instead of surface loads + # and therefore should be scaled by the surface area of all elements where + # Neumann BC are applied + "divideNeumannBoundaryConditionValuesByTotalArea": True, + # update_dirichlet_bc, # function that updates the dirichlet BCs while the simulation is running + "updateDirichletBoundaryConditionsFunction": None, + # stide every which step the update function should be called, 1 means every time step + "updateDirichletBoundaryConditionsFunctionCallInterval": 1, + # a callback function to periodically update the Neumann boundary conditions + "updateNeumannBoundaryConditionsFunction": update_neumann_bc, + # every which step the update function should be called, 1 means every time step + "updateNeumannBoundaryConditionsFunctionCallInterval": 1, + + # the initial values for the displacements, vector of values for every node [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesDisplacements": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # the initial values for the velocities, vector of values for every node [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesVelocities": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # if the initial values for the dynamic nonlinear problem should be + # computed by extrapolating the previous displacements and velocities + "extrapolateInitialGuess": True, + "constantBodyForce": variables.constant_body_force, # a constant force that acts on the whole body, e.g. for gravity + + # filename for a vtp file that contains the Dirichlet boundary condition + # nodes and their values, set to None to disable + "dirichletOutputFilename": "out/tendon_bottom_dirichlet_boundary_conditions", + # filename of a log file that will contain the total (bearing) forces and + # moments at the top and bottom of the volume + "totalForceLogFilename": "out/tendon_bottom_force.csv", + "totalForceLogOutputInterval": 10, # output interval when to write the totalForceLog file + # global element nos of the bottom elements used to compute the total forces in the log file totalForceLogFilename + "totalForceBottomElementNosGlobal": [j * nx + i for j in range(ny) for i in range(nx)], + # global element nos of the top elements used to compute the total forces + # in the log file totalForceTopElementsGlobal + "totalForceTopElementNosGlobal": [(nz - 1) * ny * nx + j * nx + i for j in range(ny) for i in range(nx)], + # callback_total_force, # callback function that gets the total force at bottom and top of the domain + "totalForceFunction": None, + "totalForceFunctionCallInterval": 1, # how often the "totalForceFunction" is called + + # define which file formats should be written + # 1. main output writer that writes output files using the quadratic + # elements function space. Writes displacements, velocities and PK2 + # stresses. + "OutputWriter": [ + + # Paraview files + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_elasticity * variables.output_timestep_3D), + "filename": "out/tendon_bottom", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + ], + # 2. additional output writer that writes also the hydrostatic pressure + "pressure": { # output files for pressure function space (linear elements), contains pressure values, as well as displacements and velocities + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/p", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, + # 3. additional output writer that writes virtual work terms + "dynamic": { # output of the dynamic solver, has additional virtual work values + "OutputWriter": [ # output files for displacements function space (quadratic elements) + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_elasticity * variables.output_timestep_3D), + "filename": "out/tendon_bottom_virtual_work", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/virtual_work", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ], + }, + # 4. output writer for debugging, outputs files after each load increment, + # the geometry is not changed but u and v are written + "LoadIncrements": { + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/load_increments", "binary": False, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, +} + +config = { + "scenarioName": variables.scenario_name, # scenario name to identify the simulation runs in the log file + "logFormat": "csv", # "csv" or "json", format of the lines in the log file, csv gives smaller files + # output file of a diagram that shows data connection between solvers + "solverStructureDiagramFile": "out/tendon_bottom_solver_structure.txt", + "mappingsBetweenMeshesLogFile": "out/tendon_bottom_mappings_between_meshes_log.txt", # log file for mappings + "Meshes": variables.meshes, + + "PreciceAdapter": { # precice adapter for bottom tendon + "timeStepOutputInterval": 100, # interval in which to display current timestep and time in console + "timestepWidth": 1, # coupling time step width, must match the value in the precice config + # if the precice coupling is enabled, if not, it simply calls the nested solver, for debugging + "couplingEnabled": True, + "preciceConfigFilename": variables.precice_config_file, + # name of the own precice participant, has to match the name given in the precice xml config file + "preciceParticipantName": "Tendon-Bottom", + "scalingFactor": 1, # a factor to scale the exchanged data, prior to communication + # if the output writers should be called only after a time window of + # precice is complete, this means the timestep has converged + "outputOnlyConvergedTimeSteps": True, + "preciceMeshes": [ # the precice meshes get created as the top or bottom surface of the main geometry mesh of the nested solver + { + "preciceMeshName": "Tendon-Bottom-Mesh", # precice name of the 2D coupling mesh + "face": "2+", # face of the 3D mesh where the 2D mesh is located, "2-" = bottom, "2+" = top + } + ], + "preciceData": [ + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "write-displacements-velocities", + # name of the precice coupling surface mesh, as given in the precice xml settings file + "preciceMeshName": "Tendon-Bottom-Mesh", + # name of the displacements "data", i.e. field variable, as given in the precice xml settings file + "displacementsName": "Displacement", + # name of the velocity "data", i.e. field variable, as given in the precice xml settings file + "velocitiesName": "Velocity", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "read-traction", + # name of the precice coupling surface mesh, as given in the precice xml settings + "preciceMeshName": "Tendon-Bottom-Mesh", + # name of the traction "data", i.e. field variable, as given in the precice xml settings file + "tractionName": "Traction", + } + ], + "HyperelasticitySolver": config_hyperelasticity, + "DynamicHyperelasticitySolver": config_hyperelasticity, + } +} diff --git a/muscle-tendon-complex/tendon-bottom-opendihu/variables/variables.py b/muscle-tendon-complex/tendon-bottom-opendihu/variables/variables.py new file mode 100644 index 000000000..91eb72a1e --- /dev/null +++ b/muscle-tendon-complex/tendon-bottom-opendihu/variables/variables.py @@ -0,0 +1,122 @@ +scenario_name = "tendon-bottom" + +# time parameters +# --------------- +dt_elasticity = 1 # [ms] time step width for elasticity +end_time = 20000 # [ms] simulation time +output_timestep_3D = 50 # [ms] output timestep + +# setup +# ----- +constant_body_force = (0, 0, -9.81e-4) # [cm/ms^2], gravity constant for the body force +force = 100.0 # [N] pulling force to the bottom + + +# input files +# ----------- + +import os +input_dir = os.environ.get('OPENDIHU_INPUT_DIR') +fiber_file = input_dir + "/left_biceps_brachii_tendon1.bin" +cellml_file = input_dir + "/2020_06_03_hodgkin-huxley_shorten_ocallaghan_davidson_soboleva_2007.cellml" +precice_config_file = "../precice-config.xml" +# If the fiber geometry data should be loaded completely in the python +# script. If True, this reads the binary file and assigns the node +# positions in the config. If False, the C++ code will read the binary +# file and only extract the local node positions. This is more performant +# for highly parallel runs. +load_fiber_data = False +debug_output = False # verbose output in this python script, for debugging the domain decomposition +disable_firing_output = True # Disables the initial list of fiber firings on the console to save some console space +paraview_output = False # If the paraview output writer should be enabled +adios_output = False # If the MegaMol/ADIOS output writer should be enabled +python_output = False # If the Python output writer should be enabled +exfile_output = False # If the Exfile output writer should be enabled# material parameters + +# material parameters +# -------------------- +tendon_material = "nonLinear" +rho = 10 + +# solvers +# ------- +diffusion_solver_type = "cg" # solver and preconditioner for the diffusion part of the Monodomain equation +diffusion_preconditioner_type = "none" # preconditioner +# solver and preconditioner for an initial Laplace flow on the domain, from which fiber directions are determined +potential_flow_solver_type = "gmres" +potential_flow_preconditioner_type = "none" # preconditioner +# solver and preconditioner for the 3D static Bidomain equation that solves the intra-muscular EMG signal +emg_solver_type = "cg" +emg_preconditioner_type = "none" # preconditioner +emg_initial_guess_nonzero = False # < If the initial guess for the emg linear system should be set to the previous solution + +# partitioning +# ------------ +# this has to match the total number of processes +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 + +# stride for sampling the 3D elements from the fiber data +# here any number is possible +sampling_stride_x = 2 +sampling_stride_y = 2 +sampling_stride_z = 50 + +mapping_tolerance = 0.1 + + +# further internal variables that will be set by the helper.py script and used in the config in settings_fibers_emg.py +n_fibers_total = None +n_subdomains_xy = None +own_subdomain_coordinate_x = None +own_subdomain_coordinate_y = None +own_subdomain_coordinate_z = None +n_fibers_x = None +n_fibers_y = None +n_points_whole_fiber = None +n_points_3D_mesh_global_x = None +n_points_3D_mesh_global_y = None +n_points_3D_mesh_global_z = None +output_writer_fibers = None +output_writer_emg = None +output_writer_0D_states = None +states_output = False +parameters_used_as_algebraic = None +parameters_used_as_constant = None +parameters_initial_values = None +output_algebraic_index = None +output_state_index = None +nodal_stimulation_current = None +fiber_file_handle = None +fibers = None +fiber_distribution = None +firing_times = None +n_fibers_per_subdomain_x = None +n_fibers_per_subdomain_y = None +n_points_per_subdomain_z = None +z_point_index_start = None +z_point_index_end = None +meshes = None +potential_flow_dirichlet_bc = None +elasticity_dirichlet_bc = None +elasticity_neumann_bc = None +fibers_on_own_rank = None +n_fiber_nodes_on_subdomain = None +fiber_start_node_no = None +generate_linear_3d_mesh = False +generate_quadratic_3d_mesh = True +nx = None +ny = None +nz = None +constant_body_force = None +bottom_traction = None +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 +states_initial_values = [] +enable_coupling = True +enable_force_length_relation = True +lambda_dot_scaling_factor = 1 +mappings = None +vm_value_stimulated = None diff --git a/muscle-tendon-complex/tendon-top-A-opendihu/clean.sh b/muscle-tendon-complex/tendon-top-A-opendihu/clean.sh new file mode 100755 index 000000000..56411170b --- /dev/null +++ b/muscle-tendon-complex/tendon-top-A-opendihu/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_opendihu . diff --git a/muscle-tendon-complex/tendon-top-A-opendihu/run.sh b/muscle-tendon-complex/tendon-top-A-opendihu/run.sh new file mode 100755 index 000000000..c7b70e79a --- /dev/null +++ b/muscle-tendon-complex/tendon-top-A-opendihu/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +./tendon-solver settings-tendon-top-A.py diff --git a/muscle-tendon-complex/tendon-top-A-opendihu/settings-tendon-top-A.py b/muscle-tendon-complex/tendon-top-A-opendihu/settings-tendon-top-A.py new file mode 100644 index 000000000..68e17404a --- /dev/null +++ b/muscle-tendon-complex/tendon-top-A-opendihu/settings-tendon-top-A.py @@ -0,0 +1,330 @@ +# Transversely-isotropic Mooney Rivlin on a tendon geometry +# Note, this is not possible to be run in parallel because the fibers +# cannot be initialized without MultipleInstances class. +import sys +import os +import numpy as np +import sys + +# set title of terminal +title = "tendon-top-a" +print('\33]0;{}\a'.format(title), end='', flush=True) + +# add variables subfolder to python path where the variables script is located +script_path = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, script_path) +sys.path.insert(0, os.path.join(script_path, 'variables')) + +import variables +from create_partitioned_meshes_for_settings import * + +# update material parameters +if (variables.tendon_material == "nonLinear"): + c = 9.98 # [N/cm^2=kPa] + ca = 14.92 # [-] + ct = 14.7 # [-] + cat = 9.64 # [-] + ctt = 11.24 # [-] + mu = 3.76 # [N/cm^2=kPa] + k1 = 42.217e3 # [N/cm^2=kPa] + k2 = 411.360e3 # [N/cm^2=kPa] + + variables.material_parameters = [c, ca, ct, cat, ctt, mu, k1, k2] + +if (variables.tendon_material == "SaintVenantKirchoff"): + # material parameters for Saint Venant-Kirchhoff material + # https://www.researchgate.net/publication/230248067_Bulk_Modulus + + youngs_modulus = 7e4 # [N/cm^2 = 10kPa] + shear_modulus = 3e4 + + lambd = shear_modulus * (youngs_modulus - 2 * shear_modulus) / \ + (3 * shear_modulus - youngs_modulus) # Lamé parameter lambda + mu = shear_modulus # Lamé parameter mu or G (shear modulus) + + variables.material_parameters = [lambd, mu] + +# If the fiber geometry data should be loaded completely in the python +# script. If True, this reads the binary file and assigns the node +# positions in the config. If False, the C++ code will read the binary +# file and only extract the local node positions. This is more performant +# for highly parallel runs. +load_fiber_data = False + +# parse arguments +rank_no = (int)(sys.argv[-2]) +n_ranks = (int)(sys.argv[-1]) + +# compute partitioning +if rank_no == 0: + if n_ranks != variables.n_subdomains_x * variables.n_subdomains_y * variables.n_subdomains_z: + print( + "\n\nError! Number of ranks {} does not match given partitioning {} x {} x {} = {}.\n\n".format( + n_ranks, + variables.n_subdomains_x, + variables.n_subdomains_y, + variables.n_subdomains_z, + variables.n_subdomains_x * + variables.n_subdomains_y * + variables.n_subdomains_z)) + sys.exit(-1) + +# stride for sampling the 3D elements from the fiber data +# here any number is possible +sampling_stride_x = 1 +sampling_stride_y = 1 +sampling_stride_z = 2 + +# create the partitioning using the script in create_partitioned_meshes_for_settings.py +result = create_partitioned_meshes_for_settings( + variables.n_subdomains_x, variables.n_subdomains_y, variables.n_subdomains_z, + variables.fiber_file, load_fiber_data, + sampling_stride_x, sampling_stride_y, sampling_stride_z, True, True) + +[variables.meshes, + variables.own_subdomain_coordinate_x, + variables.own_subdomain_coordinate_y, + variables.own_subdomain_coordinate_z, + variables.n_fibers_x, + variables.n_fibers_y, + variables.n_points_whole_fiber] = result + +n_points_3D_mesh_linear_global_x = sum([n_sampled_points_in_subdomain_x(subdomain_coordinate_x) + for subdomain_coordinate_x in range(variables.n_subdomains_x)]) +n_points_3D_mesh_linear_global_y = sum([n_sampled_points_in_subdomain_y(subdomain_coordinate_y) + for subdomain_coordinate_y in range(variables.n_subdomains_y)]) +n_points_3D_mesh_linear_global_z = sum([n_sampled_points_in_subdomain_z(subdomain_coordinate_z) + for subdomain_coordinate_z in range(variables.n_subdomains_z)]) +n_points_3D_mesh_linear_global = n_points_3D_mesh_linear_global_x * \ + n_points_3D_mesh_linear_global_y * n_points_3D_mesh_linear_global_z +nx = n_points_3D_mesh_linear_global_x - 1 +ny = n_points_3D_mesh_linear_global_y - 1 +nz = n_points_3D_mesh_linear_global_z - 1 + +node_positions = variables.meshes["3Dmesh_quadratic"]["nodePositions"] + +# boundary conditions (for quadratic elements) +# -------------------------------------------- +[mx, my, mz] = variables.meshes["3Dmesh_quadratic"]["nPointsGlobal"] +[nx, ny, nz] = variables.meshes["3Dmesh_quadratic"]["nElements"] + +# set Dirichlet BC, fix top end of tendon that is attached to the bone +variables.elasticity_dirichlet_bc = {} +k = mz - 1 + +# fix the whole x-y plane +for j in range(my): + for i in range(mx): + variables.elasticity_dirichlet_bc[k * mx * my + j * mx + i] = [0.0, 0.0, 0.0, None, None, None] + +# set no Neumann BC +variables.elasticity_neumann_bc = [] + + +config_hyperelasticity = { # for both "HyperelasticitySolver" and "DynamicHyperelasticitySolver" + "timeStepWidth": variables.dt_elasticity, # time step width + "endTime": variables.end_time, # end time of the simulation time span + "durationLogKey": "duration_mechanics", # key to find duration of this solver in the log file + "timeStepOutputInterval": 1, # how often the current time step should be printed to console + + "materialParameters": variables.material_parameters, # material parameters of the Mooney-Rivlin material + "density": variables.rho, # density of the material + # scaling factor for displacements, only set to sth. other than 1 only to + # increase visual appearance for very small displacements + "displacementsScalingFactor": 1.0, + # log file where residual norm values of the nonlinear solver will be written + "residualNormLogFilename": "out/tendon_top_a_log_residual_norm.txt", + # whether to use the analytically computed jacobian matrix in the nonlinear solver (fast) + "useAnalyticJacobian": True, + # whether to use the numerically computed jacobian matrix in the nonlinear + # solver (slow), only works with non-nested matrices, if both numeric and + # analytic are enable, it uses the analytic for the preconditioner and the + # numeric as normal jacobian + "useNumericJacobian": False, + + # whether to have extra output of matlab vectors, x,r, jacobian matrix (very slow) + "dumpDenseMatlabVariables": False, + # if useAnalyticJacobian,useNumericJacobian and dumpDenseMatlabVariables + # all all three true, the analytic and numeric jacobian matrices will get + # compared to see if there are programming errors for the analytic + # jacobian + + # mesh + "meshName": "3Dmesh_quadratic", # mesh with quadratic Lagrange ansatz functions + # boundary conditions are specified in global numberings, whereas the mesh is given in local numberings + "inputMeshIsGlobal": True, + + "fiberMeshNames": [], # fiber meshes that will be used to determine the fiber direction + # "fiberDirection": [0,0,1], # if fiberMeshNames is empty, directly set the constant fiber direction, in element coordinate system + # if fiberMeshNames and fiberDirections are empty, directly set the + # constant fiber direction, in element coordinate system + "fiberDirectionInElement": [0, 0, 1], + + # nonlinear solver + "relativeTolerance": 1e-10, # 1e-10 relative tolerance of the linear solver + "absoluteTolerance": 1e-10, # 1e-10 absolute tolerance of the residual of the linear solver + "solverType": "preonly", # type of the linear solver: cg groppcg pipecg pipecgrr cgne nash stcg gltr richardson chebyshev gmres tcqmr fcg pipefcg bcgs ibcgs fbcgs fbcgsr bcgsl cgs tfqmr cr pipecr lsqr preonly qcg bicg fgmres pipefgmres minres symmlq lgmres lcd gcr pipegcr pgmres dgmres tsirm cgls + "preconditionerType": "lu", # type of the preconditioner + "maxIterations": 1e4, # maximum number of iterations in the linear solver + "snesMaxFunctionEvaluations": 1e8, # maximum number of function iterations + "snesMaxIterations": 240, # maximum number of iterations in the nonlinear solver + "snesRelativeTolerance": 1e-2, # relative tolerance of the nonlinear solver + # type of linesearch, possible values: "bt" "nleqerr" "basic" "l2" "cp" "ncglinear" + "snesLineSearchType": "l2", + "snesAbsoluteTolerance": 1e-5, # absolute tolerance of the nonlinear solver + # how often the jacobian should be recomputed, -1 indicates NEVER rebuild, + # 1 means rebuild every time the Jacobian is computed within a single + # nonlinear solve, 2 means every second time the Jacobian is built etc. -2 + # means rebuild at next chance but then never again + "snesRebuildJacobianFrequency": 5, + + # "dumpFilename": "out/r{}/m".format(sys.argv[-1]), # dump system matrix and right hand side after every solve + "dumpFilename": "", # dump disabled + "dumpFormat": "matlab", # default, ascii, matlab + + # "loadFactors": [0.1, 0.2, 0.35, 0.5, 1.0], # load factors for every timestep + # "loadFactors": [0.5, 1.0], # load factors for every timestep + "loadFactors": [], # no load factors, solve problem directly + # a threshold for the load factor, when to abort the solve of the current + # time step. The load factors are adjusted automatically if the nonlinear + # solver diverged. If the load factors get too small, it aborts the solve. + "loadFactorGiveUpThreshold": 1e-3, + "nNonlinearSolveCalls": 1, # how often the nonlinear solve should be called + + # boundary and initial conditions + # the initial Dirichlet boundary conditions that define values for displacements u and velocity v + "dirichletBoundaryConditions": variables.elasticity_dirichlet_bc, + # Neumann boundary conditions that define traction forces on surfaces of elements + "neumannBoundaryConditions": variables.elasticity_neumann_bc, + # if the given Neumann boundary condition values under + # "neumannBoundaryConditions" are total forces instead of surface loads + # and therefore should be scaled by the surface area of all elements where + # Neumann BC are applied + "divideNeumannBoundaryConditionValuesByTotalArea": False, + # function that updates the dirichlet BCs while the simulation is running + "updateDirichletBoundaryConditionsFunction": None, + # every which step the update function should be called, 1 means every time step + "updateDirichletBoundaryConditionsFunctionCallInterval": 1, + + # the initial values for the displacements, vector of values for every node [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesDisplacements": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # the initial values for the velocities, vector of values for every node [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesVelocities": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # if the initial values for the dynamic nonlinear problem should be + # computed by extrapolating the previous displacements and velocities + "extrapolateInitialGuess": False, + "constantBodyForce": variables.constant_body_force, # a constant force that acts on the whole body, e.g. for gravity + + # filename for a vtp file that contains the Dirichlet boundary condition + # nodes and their values, set to None to disable + "dirichletOutputFilename": "out/tendon_top_a_dirichlet_boundary_conditions_tendon_top_a", + # filename of a log file that will contain the total (bearing) forces and + # moments at the top and bottom of the volume + "totalForceLogFilename": "out/tendon_top_a_force.csv", + "totalForceLogOutputInterval": 10, # output interval when to write the totalForceLog file + # global element nos of the bottom elements used to compute the total forces in the log file totalForceLogFilename + "totalForceBottomElementNosGlobal": [j * nx + i for j in range(ny) for i in range(nx)], + # global element nos of the top elements used to compute the total forces + # in the log file totalForceTopElementsGlobal + "totalForceTopElementNosGlobal": [(nz - 1) * ny * nx + j * nx + i for j in range(ny) for i in range(nx)], + + # define which file formats should be written + # 1. main output writer that writes output files using the quadratic + # elements function space. Writes displacements, velocities and PK2 + # stresses. + "OutputWriter": [ + + # Paraview files + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_elasticity * variables.output_timestep_3D), + "filename": "out/tendon_top_a", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + + # Python callback function "postprocess" + # {"format": "PythonCallback", "outputInterval": 1, "callback": postprocess, "onlyNodalValues":True, "filename": ""}, + ], + # 2. additional output writer that writes also the hydrostatic pressure + "pressure": { # output files for pressure function space (linear elements), contains pressure values, as well as displacements and velocities + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/p", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, + # 3. additional output writer that writes virtual work terms + "dynamic": { # output of the dynamic solver, has additional virtual work values + "OutputWriter": [ # output files for displacements function space (quadratic elements) + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_elasticity * variables.output_timestep_3D), + "filename": "out/tendon_top_a_virtual_work", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/virtual_work", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ], + }, + # 4. output writer for debugging, outputs files after each load increment, + # the geometry is not changed but u and v are written + "LoadIncrements": { + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/load_increments", "binary": False, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, +} + +config = { + "scenarioName": variables.scenario_name, # scenario name to identify the simulation runs in the log file + "logFormat": "csv", # "csv" or "json", format of the lines in the log file, csv gives smaller files + # output file of a diagram that shows data connection between solvers + "solverStructureDiagramFile": "out/tendon_top_a_solver_structure.txt", + "mappingsBetweenMeshesLogFile": "out/tendon_top_a_mappings_between_meshes_log.txt", # log file for mappings + "Meshes": variables.meshes, + + "PreciceAdapter": { # precice adapter for bottom tendon + "timeStepOutputInterval": 100, # interval in which to display current timestep and time in console + "timestepWidth": 1, # coupling time step width, must match the value in the precice config + # if the precice coupling is enabled, if not, it simply calls the nested solver, for debugging + "couplingEnabled": True, + "preciceConfigFilename": variables.precice_config_file, # the preCICE configuration file + # name of the own precice participant, has to match the name given in the precice xml config file + "preciceParticipantName": "Tendon-Top-A", + "scalingFactor": 1, # a factor to scale the exchanged data, prior to communication + # if the output writers should be called only after a time window of + # precice is complete, this means the timestep has converged + "outputOnlyConvergedTimeSteps": True, + "preciceMeshes": [ # the precice meshes get created as the top or bottom surface of the main geometry mesh of the nested solver + { + "preciceMeshName": "Tendon-Top-A-Mesh", # precice name of the 2D coupling mesh + "face": "2-", # face of the 3D mesh where the 2D mesh is located, "2-" = bottom, "2+" = top + } + ], + "preciceData": [ + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "write-displacements-velocities", + # name of the precice coupling surface mesh, as given in the precice xml settings file + "preciceMeshName": "Tendon-Top-A-Mesh", + # name of the displacements "data", i.e. field variable, as given in the precice xml settings file + "displacementsName": "Displacement", + # name of the velocity "data", i.e. field variable, as given in the precice xml settings file + "velocitiesName": "Velocity", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "read-traction", + # name of the precice coupling surface mesh, as given in the precice xml settings + "preciceMeshName": "Tendon-Top-A-Mesh", + # name of the traction "data", i.e. field variable, as given in the precice xml settings file + "tractionName": "Traction", + } + ], + "HyperelasticitySolver": config_hyperelasticity, + "DynamicHyperelasticitySolver": config_hyperelasticity, + } +} diff --git a/muscle-tendon-complex/tendon-top-A-opendihu/variables/variables.py b/muscle-tendon-complex/tendon-top-A-opendihu/variables/variables.py new file mode 100644 index 000000000..433abb485 --- /dev/null +++ b/muscle-tendon-complex/tendon-top-A-opendihu/variables/variables.py @@ -0,0 +1,122 @@ +scenario_name = "tendon-top-A" + +# time parameters +# --------------- +dt_elasticity = 1 # [ms] time step width for elasticity +end_time = 20000 # [ms] simulation time +output_timestep_3D = 50 # [ms] output timestep + +# setup +# ----- +constant_body_force = (0, 0, -9.81e-4) # [cm/ms^2], gravity constant for the body force +force = 100.0 # [N] pulling force to the bottom + + +# input files +# ----------- + +import os +input_dir = os.environ.get('OPENDIHU_INPUT_DIR') +fiber_file = input_dir + "/left_biceps_brachii_tendon2a.bin" +cellml_file = input_dir + "/2020_06_03_hodgkin-huxley_shorten_ocallaghan_davidson_soboleva_2007.cellml" +precice_config_file = "../precice-config.xml" +# If the fiber geometry data should be loaded completely in the python +# script. If True, this reads the binary file and assigns the node +# positions in the config. If False, the C++ code will read the binary +# file and only extract the local node positions. This is more performant +# for highly parallel runs. +load_fiber_data = False +debug_output = False # verbose output in this python script, for debugging the domain decomposition +disable_firing_output = True # Disables the initial list of fiber firings on the console to save some console space +paraview_output = False # If the paraview output writer should be enabled +adios_output = False # If the MegaMol/ADIOS output writer should be enabled +python_output = False # If the Python output writer should be enabled +exfile_output = False # If the Exfile output writer should be enabled# material parameters + +# material parameters +# -------------------- +tendon_material = "nonLinear" +rho = 10 + +# solvers +# ------- +diffusion_solver_type = "cg" # solver and preconditioner for the diffusion part of the Monodomain equation +diffusion_preconditioner_type = "none" # preconditioner +# solver and preconditioner for an initial Laplace flow on the domain, from which fiber directions are determined +potential_flow_solver_type = "gmres" +potential_flow_preconditioner_type = "none" # preconditioner +# solver and preconditioner for the 3D static Bidomain equation that solves the intra-muscular EMG signal +emg_solver_type = "cg" +emg_preconditioner_type = "none" # preconditioner +emg_initial_guess_nonzero = False # < If the initial guess for the emg linear system should be set to the previous solution + +# partitioning +# ------------ +# this has to match the total number of processes +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 + +# stride for sampling the 3D elements from the fiber data +# here any number is possible +sampling_stride_x = 2 +sampling_stride_y = 2 +sampling_stride_z = 50 + +mapping_tolerance = 0.1 + + +# further internal variables that will be set by the helper.py script and used in the config in settings_fibers_emg.py +n_fibers_total = None +n_subdomains_xy = None +own_subdomain_coordinate_x = None +own_subdomain_coordinate_y = None +own_subdomain_coordinate_z = None +n_fibers_x = None +n_fibers_y = None +n_points_whole_fiber = None +n_points_3D_mesh_global_x = None +n_points_3D_mesh_global_y = None +n_points_3D_mesh_global_z = None +output_writer_fibers = None +output_writer_emg = None +output_writer_0D_states = None +states_output = False +parameters_used_as_algebraic = None +parameters_used_as_constant = None +parameters_initial_values = None +output_algebraic_index = None +output_state_index = None +nodal_stimulation_current = None +fiber_file_handle = None +fibers = None +fiber_distribution = None +firing_times = None +n_fibers_per_subdomain_x = None +n_fibers_per_subdomain_y = None +n_points_per_subdomain_z = None +z_point_index_start = None +z_point_index_end = None +meshes = None +potential_flow_dirichlet_bc = None +elasticity_dirichlet_bc = None +elasticity_neumann_bc = None +fibers_on_own_rank = None +n_fiber_nodes_on_subdomain = None +fiber_start_node_no = None +generate_linear_3d_mesh = False +generate_quadratic_3d_mesh = True +nx = None +ny = None +nz = None +constant_body_force = None +bottom_traction = None +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 +states_initial_values = [] +enable_coupling = True +enable_force_length_relation = True +lambda_dot_scaling_factor = 1 +mappings = None +vm_value_stimulated = None diff --git a/muscle-tendon-complex/tendon-top-B-opendihu/clean.sh b/muscle-tendon-complex/tendon-top-B-opendihu/clean.sh new file mode 100755 index 000000000..56411170b --- /dev/null +++ b/muscle-tendon-complex/tendon-top-B-opendihu/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_opendihu . diff --git a/muscle-tendon-complex/tendon-top-B-opendihu/run.sh b/muscle-tendon-complex/tendon-top-B-opendihu/run.sh new file mode 100755 index 000000000..ad204d701 --- /dev/null +++ b/muscle-tendon-complex/tendon-top-B-opendihu/run.sh @@ -0,0 +1,4 @@ +#!/bin/sh +set -e -u + +./tendon-solver settings-tendon-top-B.py diff --git a/muscle-tendon-complex/tendon-top-B-opendihu/settings-tendon-top-B.py b/muscle-tendon-complex/tendon-top-B-opendihu/settings-tendon-top-B.py new file mode 100644 index 000000000..8a6a63b7f --- /dev/null +++ b/muscle-tendon-complex/tendon-top-B-opendihu/settings-tendon-top-B.py @@ -0,0 +1,332 @@ +# Transversely-isotropic Mooney Rivlin on a tendon geometry +# Note, this is not possible to be run in parallel because the fibers +# cannot be initialized without MultipleInstances class. +import sys +import os +import numpy as np +import sys + +# set title of terminal +title = "tendon-top-b" +print('\33]0;{}\a'.format(title), end='', flush=True) + +# add variables subfolder to python path where the variables script is located +script_path = os.path.dirname(os.path.abspath(__file__)) +sys.path.insert(0, script_path) +sys.path.insert(0, os.path.join(script_path, 'variables')) + +import variables +from create_partitioned_meshes_for_settings import * + +# update material parameters +if (variables.tendon_material == "nonLinear"): + c = 9.98 # [N/cm^2=kPa] + ca = 14.92 # [-] + ct = 14.7 # [-] + cat = 9.64 # [-] + ctt = 11.24 # [-] + mu = 3.76 # [N/cm^2=kPa] + k1 = 42.217e3 # [N/cm^2=kPa] + k2 = 411.360e3 # [N/cm^2=kPa] + + variables.material_parameters = [c, ca, ct, cat, ctt, mu, k1, k2] + +if (variables.tendon_material == "SaintVenantKirchoff"): + # material parameters for Saint Venant-Kirchhoff material + # https://www.researchgate.net/publication/230248067_Bulk_Modulus + + youngs_modulus = 7e4 # [N/cm^2 = 10kPa] + shear_modulus = 3e4 + + lambd = shear_modulus * (youngs_modulus - 2 * shear_modulus) / \ + (3 * shear_modulus - youngs_modulus) # Lamé parameter lambda + mu = shear_modulus # Lamé parameter mu or G (shear modulus) + + variables.material_parameters = [lambd, mu] + +# If the fiber geometry data should be loaded completely in the python +# script. If True, this reads the binary file and assigns the node +# positions in the config. If False, the C++ code will read the binary +# file and only extract the local node positions. This is more performant +# for highly parallel runs. +load_fiber_data = False + +# parse arguments +rank_no = (int)(sys.argv[-2]) +n_ranks = (int)(sys.argv[-1]) + +# partitioning +# ------------ + +# compute partitioning +if rank_no == 0: + if n_ranks != variables.n_subdomains_x * variables.n_subdomains_y * variables.n_subdomains_z: + print( + "\n\nError! Number of ranks {} does not match given partitioning {} x {} x {} = {}.\n\n".format( + n_ranks, + variables.n_subdomains_x, + variables.n_subdomains_y, + variables.n_subdomains_z, + variables.n_subdomains_x * + variables.n_subdomains_y * + variables.n_subdomains_z)) + sys.exit(-1) + +# stride for sampling the 3D elements from the fiber data +# here any number is possible +sampling_stride_x = 1 +sampling_stride_y = 1 +sampling_stride_z = 2 + +# create the partitioning using the script in create_partitioned_meshes_for_settings.py +result = create_partitioned_meshes_for_settings( + variables.n_subdomains_x, variables.n_subdomains_y, variables.n_subdomains_z, + variables.fiber_file, load_fiber_data, + sampling_stride_x, sampling_stride_y, sampling_stride_z, True, True) + +[variables.meshes, + variables.own_subdomain_coordinate_x, + variables.own_subdomain_coordinate_y, + variables.own_subdomain_coordinate_z, + variables.n_fibers_x, + variables.n_fibers_y, + variables.n_points_whole_fiber] = result + +n_points_3D_mesh_linear_global_x = sum([n_sampled_points_in_subdomain_x(subdomain_coordinate_x) + for subdomain_coordinate_x in range(variables.n_subdomains_x)]) +n_points_3D_mesh_linear_global_y = sum([n_sampled_points_in_subdomain_y(subdomain_coordinate_y) + for subdomain_coordinate_y in range(variables.n_subdomains_y)]) +n_points_3D_mesh_linear_global_z = sum([n_sampled_points_in_subdomain_z(subdomain_coordinate_z) + for subdomain_coordinate_z in range(variables.n_subdomains_z)]) +n_points_3D_mesh_linear_global = n_points_3D_mesh_linear_global_x * \ + n_points_3D_mesh_linear_global_y * n_points_3D_mesh_linear_global_z +nx = n_points_3D_mesh_linear_global_x - 1 +ny = n_points_3D_mesh_linear_global_y - 1 +nz = n_points_3D_mesh_linear_global_z - 1 + +node_positions = variables.meshes["3Dmesh_quadratic"]["nodePositions"] + +# boundary conditions (for quadratic elements) +# -------------------------------------------- +[mx, my, mz] = variables.meshes["3Dmesh_quadratic"]["nPointsGlobal"] +[nx, ny, nz] = variables.meshes["3Dmesh_quadratic"]["nElements"] + +# set Dirichlet BC, fix top end of tendon that is attached to the bone +variables.elasticity_dirichlet_bc = {} +k = mz - 1 + +# fix the whole x-y plane +for j in range(my): + for i in range(mx): + variables.elasticity_dirichlet_bc[k * mx * my + j * mx + i] = [0.0, 0.0, 0.0, None, None, None] + +# set no Neumann BC +variables.elasticity_neumann_bc = [] + +config_hyperelasticity = { # for both "HyperelasticitySolver" and "DynamicHyperelasticitySolver" + "timeStepWidth": variables.dt_elasticity, # time step width + "endTime": variables.end_time, # end time of the simulation time span + "durationLogKey": "duration_mechanics", # key to find duration of this solver in the log file + "timeStepOutputInterval": 1, # how often the current time step should be printed to console + + "materialParameters": variables.material_parameters, # material parameters of the Mooney-Rivlin material + "density": variables.rho, # density of the material + # scaling factor for displacements, only set to sth. other than 1 only to + # increase visual appearance for very small displacements + "displacementsScalingFactor": 1.0, + # log file where residual norm values of the nonlinear solver will be written + "residualNormLogFilename": "out/tendon_top_b_log_residual_norm.txt", + # whether to use the analytically computed jacobian matrix in the nonlinear solver (fast) + "useAnalyticJacobian": True, + # whether to use the numerically computed jacobian matrix in the nonlinear + # solver (slow), only works with non-nested matrices, if both numeric and + # analytic are enable, it uses the analytic for the preconditioner and the + # numeric as normal jacobian + "useNumericJacobian": False, + + # whether to have extra output of matlab vectors, x,r, jacobian matrix (very slow) + "dumpDenseMatlabVariables": False, + # if useAnalyticJacobian,useNumericJacobian and dumpDenseMatlabVariables + # all all three true, the analytic and numeric jacobian matrices will get + # compared to see if there are programming errors for the analytic + # jacobian + + # mesh + "meshName": "3Dmesh_quadratic", # mesh with quadratic Lagrange ansatz functions + # boundary conditions are specified in global numberings, whereas the mesh is given in local numberings + "inputMeshIsGlobal": True, + + "fiberMeshNames": [], # fiber meshes that will be used to determine the fiber direction + # "fiberDirection": [0,0,1], # if fiberMeshNames is empty, directly set the constant fiber direction, in element coordinate system + # if fiberMeshNames and fiberDirections are empty, directly set the + # constant fiber direction, in element coordinate system + "fiberDirectionInElement": [0, 0, 1], + + # nonlinear solver + "relativeTolerance": 1e-10, # 1e-10 relative tolerance of the linear solver + "absoluteTolerance": 1e-10, # 1e-10 absolute tolerance of the residual of the linear solver + "solverType": "preonly", # type of the linear solver: cg groppcg pipecg pipecgrr cgne nash stcg gltr richardson chebyshev gmres tcqmr fcg pipefcg bcgs ibcgs fbcgs fbcgsr bcgsl cgs tfqmr cr pipecr lsqr preonly qcg bicg fgmres pipefgmres minres symmlq lgmres lcd gcr pipegcr pgmres dgmres tsirm cgls + "preconditionerType": "lu", # type of the preconditioner + "maxIterations": 1e4, # maximum number of iterations in the linear solver + "snesMaxFunctionEvaluations": 1e8, # maximum number of function iterations + "snesMaxIterations": 240, # maximum number of iterations in the nonlinear solver + "snesRelativeTolerance": 1e-2, # relative tolerance of the nonlinear solver + # type of linesearch, possible values: "bt" "nleqerr" "basic" "l2" "cp" "ncglinear" + "snesLineSearchType": "l2", + "snesAbsoluteTolerance": 1e-5, # absolute tolerance of the nonlinear solver + # how often the jacobian should be recomputed, -1 indicates NEVER rebuild, + # 1 means rebuild every time the Jacobian is computed within a single + # nonlinear solve, 2 means every second time the Jacobian is built etc. -2 + # means rebuild at next chance but then never again + "snesRebuildJacobianFrequency": 5, + + # "dumpFilename": "out/r{}/m".format(sys.argv[-1]), # dump system matrix and right hand side after every solve + "dumpFilename": "", # dump disabled + "dumpFormat": "matlab", # default, ascii, matlab + + # "loadFactors": [0.1, 0.2, 0.35, 0.5, 1.0], # load factors for every timestep + # "loadFactors": [0.5, 1.0], # load factors for every timestep + "loadFactors": [], # no load factors, solve problem directly + # a threshold for the load factor, when to abort the solve of the current + # time step. The load factors are adjusted automatically if the nonlinear + # solver diverged. If the load factors get too small, it aborts the solve. + "loadFactorGiveUpThreshold": 1e-3, + "nNonlinearSolveCalls": 1, # how often the nonlinear solve should be called + + # boundary and initial conditions + # the initial Dirichlet boundary conditions that define values for displacements u and velocity v + "dirichletBoundaryConditions": variables.elasticity_dirichlet_bc, + # Neumann boundary conditions that define traction forces on surfaces of elements + "neumannBoundaryConditions": variables.elasticity_neumann_bc, + # if the given Neumann boundary condition values under + # "neumannBoundaryConditions" are total forces instead of surface loads + # and therefore should be scaled by the surface area of all elements where + # Neumann BC are applied + "divideNeumannBoundaryConditionValuesByTotalArea": False, + # function that updates the dirichlet BCs while the simulation is running + "updateDirichletBoundaryConditionsFunction": None, + # every which step the update function should be called, 1 means every time step + "updateDirichletBoundaryConditionsFunctionCallInterval": 1, + + # the initial values for the displacements, vector of values for every node [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesDisplacements": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # the initial values for the velocities, vector of values for every node [[node1-x,y,z], [node2-x,y,z], ...] + "initialValuesVelocities": [[0.0, 0.0, 0.0] for _ in range(mx * my * mz)], + # if the initial values for the dynamic nonlinear problem should be + # computed by extrapolating the previous displacements and velocities + "extrapolateInitialGuess": False, + "constantBodyForce": variables.constant_body_force, # a constant force that acts on the whole body, e.g. for gravity + + # filename for a vtp file that contains the Dirichlet boundary condition + # nodes and their values, set to None to disable + "dirichletOutputFilename": "out/tendon_top_b_dirichlet_boundary_conditions", + # filename of a log file that will contain the total (bearing) forces and + # moments at the top and bottom of the volume + "totalForceLogFilename": "out/tendon_top_b_force.csv", + "totalForceLogOutputInterval": 10, # output interval when to write the totalForceLog file + # global element nos of the bottom elements used to compute the total forces in the log file totalForceLogFilename + "totalForceBottomElementNosGlobal": [j * nx + i for j in range(ny) for i in range(nx)], + # global element nos of the top elements used to compute the total forces + # in the log file totalForceTopElementsGlobal + "totalForceTopElementNosGlobal": [(nz - 1) * ny * nx + j * nx + i for j in range(ny) for i in range(nx)], + + # define which file formats should be written + # 1. main output writer that writes output files using the quadratic + # elements function space. Writes displacements, velocities and PK2 + # stresses. + "OutputWriter": [ + + # Paraview files + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_elasticity * variables.output_timestep_3D), + "filename": "out/tendon_top_b", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + + # Python callback function "postprocess" + # {"format": "PythonCallback", "outputInterval": 1, "callback": postprocess, "onlyNodalValues":True, "filename": ""}, + ], + # 2. additional output writer that writes also the hydrostatic pressure + "pressure": { # output files for pressure function space (linear elements), contains pressure values, as well as displacements and velocities + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/p", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, + # 3. additional output writer that writes virtual work terms + "dynamic": { # output of the dynamic solver, has additional virtual work values + "OutputWriter": [ # output files for displacements function space (quadratic elements) + {"format": "Paraview", + "outputInterval": int(1. / variables.dt_elasticity * variables.output_timestep_3D), + "filename": "out/tendon_top_b_virtual_work", + "binary": True, + "fixedFormat": False, + "onlyNodalValues": True, + "combineFiles": True, + "fileNumbering": "incremental"}, + # {"format": "Paraview", "outputInterval": 1, "filename": "out/"+variables.scenario_name+"/virtual_work", "binary": True, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ], + }, + # 4. output writer for debugging, outputs files after each load increment, + # the geometry is not changed but u and v are written + "LoadIncrements": { + "OutputWriter": [ + # {"format": "Paraview", "outputInterval": 1, "filename": "out/load_increments", "binary": False, "fixedFormat": False, "onlyNodalValues":True, "combineFiles":True, "fileNumbering": "incremental"}, + ] + }, +} + +config = { + "scenarioName": variables.scenario_name, # scenario name to identify the simulation runs in the log file + "logFormat": "csv", # "csv" or "json", format of the lines in the log file, csv gives smaller files + # output file of a diagram that shows data connection between solvers + "solverStructureDiagramFile": "out/tendon_top_b_solver_structure.txt", + "mappingsBetweenMeshesLogFile": "out/tendon_top_b_mappings_between_meshes_log.txt", # log file for mappings + "Meshes": variables.meshes, + + "PreciceAdapter": { # precice adapter for bottom tendon + "timeStepOutputInterval": 100, # interval in which to display current timestep and time in console + "timestepWidth": 1, # coupling time step width, must match the value in the precice config + # if the precice coupling is enabled, if not, it simply calls the nested solver, for debugging + "couplingEnabled": True, + "preciceConfigFilename": variables.precice_config_file, # the preCICE configuration file + # name of the own precice participant, has to match the name given in the precice xml config file + "preciceParticipantName": "Tendon-Top-B", + "scalingFactor": 1, # a factor to scale the exchanged data, prior to communication + # if the output writers should be called only after a time window of + # precice is complete, this means the timestep has converged + "outputOnlyConvergedTimeSteps": True, + "preciceMeshes": [ # the precice meshes get created as the top or bottom surface of the main geometry mesh of the nested solver + { + "preciceMeshName": "Tendon-Top-B-Mesh", # precice name of the 2D coupling mesh + "face": "2-", # face of the 3D mesh where the 2D mesh is located, "2-" = bottom, "2+" = top + } + ], + "preciceData": [ + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "write-displacements-velocities", + # name of the precice coupling surface mesh, as given in the precice xml settings file + "preciceMeshName": "Tendon-Top-B-Mesh", + # name of the displacements "data", i.e. field variable, as given in the precice xml settings file + "displacementsName": "Displacement", + # name of the velocity "data", i.e. field variable, as given in the precice xml settings file + "velocitiesName": "Velocity", + }, + { + # mode is one of "read-displacements-velocities", "read-traction", + # "write-displacements-velocities", "write-traction" + "mode": "read-traction", + # name of the precice coupling surface mesh, as given in the precice xml settings + "preciceMeshName": "Tendon-Top-B-Mesh", + # name of the traction "data", i.e. field variable, as given in the precice xml settings file + "tractionName": "Traction", + } + ], + "HyperelasticitySolver": config_hyperelasticity, + "DynamicHyperelasticitySolver": config_hyperelasticity, + } +} diff --git a/muscle-tendon-complex/tendon-top-B-opendihu/variables/variables.py b/muscle-tendon-complex/tendon-top-B-opendihu/variables/variables.py new file mode 100644 index 000000000..eca352d9d --- /dev/null +++ b/muscle-tendon-complex/tendon-top-B-opendihu/variables/variables.py @@ -0,0 +1,122 @@ +scenario_name = "tendon-top-B" + +# time parameters +# --------------- +dt_elasticity = 1 # [ms] time step width for elasticity +end_time = 20000 # [ms] simulation time +output_timestep_3D = 50 # [ms] output timestep + +# setup +# ----- +constant_body_force = (0, 0, -9.81e-4) # [cm/ms^2], gravity constant for the body force +force = 100.0 # [N] pulling force to the bottom + + +# input files +# ----------- + +import os +input_dir = os.environ.get('OPENDIHU_INPUT_DIR') +fiber_file = input_dir + "/left_biceps_brachii_tendon2b.bin" +cellml_file = input_dir + "/2020_06_03_hodgkin-huxley_shorten_ocallaghan_davidson_soboleva_2007.cellml" +precice_config_file = "../precice-config.xml" +# If the fiber geometry data should be loaded completely in the python +# script. If True, this reads the binary file and assigns the node +# positions in the config. If False, the C++ code will read the binary +# file and only extract the local node positions. This is more performant +# for highly parallel runs. +load_fiber_data = False +debug_output = False # verbose output in this python script, for debugging the domain decomposition +disable_firing_output = True # Disables the initial list of fiber firings on the console to save some console space +paraview_output = False # If the paraview output writer should be enabled +adios_output = False # If the MegaMol/ADIOS output writer should be enabled +python_output = False # If the Python output writer should be enabled +exfile_output = False # If the Exfile output writer should be enabled# material parameters + +# material parameters +# -------------------- +tendon_material = "nonLinear" +rho = 10 + +# solvers +# ------- +diffusion_solver_type = "cg" # solver and preconditioner for the diffusion part of the Monodomain equation +diffusion_preconditioner_type = "none" # preconditioner +# solver and preconditioner for an initial Laplace flow on the domain, from which fiber directions are determined +potential_flow_solver_type = "gmres" +potential_flow_preconditioner_type = "none" # preconditioner +# solver and preconditioner for the 3D static Bidomain equation that solves the intra-muscular EMG signal +emg_solver_type = "cg" +emg_preconditioner_type = "none" # preconditioner +emg_initial_guess_nonzero = False # < If the initial guess for the emg linear system should be set to the previous solution + +# partitioning +# ------------ +# this has to match the total number of processes +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 + +# stride for sampling the 3D elements from the fiber data +# here any number is possible +sampling_stride_x = 2 +sampling_stride_y = 2 +sampling_stride_z = 50 + +mapping_tolerance = 0.1 + + +# further internal variables that will be set by the helper.py script and used in the config in settings_fibers_emg.py +n_fibers_total = None +n_subdomains_xy = None +own_subdomain_coordinate_x = None +own_subdomain_coordinate_y = None +own_subdomain_coordinate_z = None +n_fibers_x = None +n_fibers_y = None +n_points_whole_fiber = None +n_points_3D_mesh_global_x = None +n_points_3D_mesh_global_y = None +n_points_3D_mesh_global_z = None +output_writer_fibers = None +output_writer_emg = None +output_writer_0D_states = None +states_output = False +parameters_used_as_algebraic = None +parameters_used_as_constant = None +parameters_initial_values = None +output_algebraic_index = None +output_state_index = None +nodal_stimulation_current = None +fiber_file_handle = None +fibers = None +fiber_distribution = None +firing_times = None +n_fibers_per_subdomain_x = None +n_fibers_per_subdomain_y = None +n_points_per_subdomain_z = None +z_point_index_start = None +z_point_index_end = None +meshes = None +potential_flow_dirichlet_bc = None +elasticity_dirichlet_bc = None +elasticity_neumann_bc = None +fibers_on_own_rank = None +n_fiber_nodes_on_subdomain = None +fiber_start_node_no = None +generate_linear_3d_mesh = False +generate_quadratic_3d_mesh = True +nx = None +ny = None +nz = None +constant_body_force = None +bottom_traction = None +n_subdomains_x = 1 +n_subdomains_y = 1 +n_subdomains_z = 1 +states_initial_values = [] +enable_coupling = True +enable_force_length_relation = True +lambda_dot_scaling_factor = 1 +mappings = None +vm_value_stimulated = None diff --git a/tools/cleaning-tools.sh b/tools/cleaning-tools.sh index 6edaf3804..80db82242 100755 --- a/tools/cleaning-tools.sh +++ b/tools/cleaning-tools.sh @@ -153,3 +153,13 @@ clean_dumux() { ) } + +clean_opendihu() { + ( + set -e -u + cd "$1" + echo "- Cleaning up OpenDiHu case in $(pwd)" + rm -rfv ./logs/ ./out/ + clean_precice_logs . + ) +}