forked from SimVascular/simvascular.github.io-archive
-
Notifications
You must be signed in to change notification settings - Fork 0
/
archiveFlowSolverSV2.html
2949 lines (2495 loc) · 126 KB
/
archiveFlowSolverSV2.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1">
<meta name="description" content="">
<meta name="author" content="">
<title>SimVascular Docs</title>
<link href="css/bootstrap.css" rel="stylesheet" type="text/css" />
<link href="css/shop-item.css" rel="stylesheet" type="text/css" />
<link href="css/codestyle.css" rel="stylesheet" type="text/css" />
<link rel="stylesheet" href="font-awesome-4.1.0/css/font-awesome.min.css">
<link rel="stylesheet" href="https://code.ionicframework.com/ionicons/1.5.2/css/ionicons.min.css">
<link rel="shortcut icon" href="img/favicon.ico">
<!-- HTML5 Shim and Respond.js IE8 support of HTML5 elements and media queries -->
<!-- WARNING: Respond.js doesn't work if you view the page via file:// -->
<!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/libs/html5shiv/3.7.0/html5shiv.js"></script>
<script src="https://oss.maxcdn.com/libs/respond.js/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body>
<!-- Navigation -->
<nav class="navbar navbar-inverse navbar-fixed-top" role="navigation">
<div class="container">
<!-- Brand and toggle get grouped for better mobile display -->
<div class="navbar-header">
<button type="button" class="navbar-toggle" data-toggle="collapse" data-target=".navbar-main-collapse">
<i class="fa fa-bars" id="barIcon"></i>
</button>
<a class="navbar-brand" id="brandName" href="index.html">
<img src="img/svlogo/svLogoSmallText.png" alt="...">
</a>
<a class="navbar-brand page-scroll">
2.0
</a>
</div>
<!-- Collect the nav links, forms, and other content for toggling -->
<div class="collapse navbar-collapse navbar-right navbar-main-collapse">
<ul class="nav navbar-nav">
<!-- USER GUIDES -->
<li>
<a href="#" id="dropdownMenu1" data-toggle="dropdown">
<b><span class="fa fa-user"></span> User Guides </b>
</a>
<ul class="dropdown-menu" role="menu" aria-labelledby="dropdownMenu1">
<li role="presentation"><a role="menuitem" tabindex="-1" href="archiveInstallationSV2.html"><b><span class="fa fa-sign-in fa-rotate-90"></span> Installation</b></a></li>
<li role="presentation" class="divider"></li>
<li role="presentation"><a role="menuitem" tabindex="-1" href="archiveQuickGuideSV2.html"><b><span class="icon ion-ios7-bolt"></span> Quick Guide</b></a></li>
<li role="presentation" class="divider"></li>
<li role="presentation"><a role="menuitem" tabindex="-1" href="archiveImageGuideSV2.html"><b><span class="icon ion-settings"></span> Image Visualization</b></a></li>
<li role="presentation"><a role="menuitem" tabindex="-1" href="archiveModelGuideSV2.html"><b><span class="icon ion-settings"></span> Modeling</b></a></li>
<li role="presentation"><a role="menuitem" tabindex="-1" href="archiveMeshingSV2.html"><b><span class="icon ion-ios7-keypad"></span> Meshing</b></a></li>
<li role="presentation"><a role="menuitem" tabindex="-1" href="archiveFlowSolverSV2.html"><b><span class="icon ion-play"></span> Simulation</b></a></li>
<li role="presentation" class="divider"></li>
<li role="presentation"><a role="menuitem" tabindex="-1" href="archiveCodeSnippetsSV2.html"><b><span class="icon ion-play"></span> Code Snippets</b></a></li>
</ul>
</li>
<!-- Latest version -->
<li>
<a href="docsQuickGuide.html" id="dropdownMenu1" >
<b>Back To Lastest Version</b>
</a>
</li>
<!-- <li><a href="docsQuickGuide.html" id="btnQuickGuide"><b><span class="icon ion-ios7-bolt"></span> Quick Guide</b></a></li>
<li><a href="docsModelGuide.html" id="btnModelGuide"><b><span class="icon ion-settings"></span> Modeling</b></a></li>
<li><a href="docsMeshing.html" id="btnMeshing"><b><span class="icon ion-ios7-keypad"></span> Meshing</b></a></li>
<li><a href="docsPresolver.html" id="btnPresolver"><b><span class="icon ion-log-in"></span> svPre</b></a></li>
<li><a href="docsFlowSolver.html" id="btnFlowSolver"><b><span class="icon ion-play"></span> svSolver</b></a></li>
<li><a href="docsRefs.html" id="btnRefs"><b><span class="icon ion-document-text"></span> References</b></a></li>
<li><a href="clinicalCase1.html" id="btnRefs"><b>Case Studies</b></a></li> -->
</ul>
</div>
<!-- /.navbar-collapse -->
</div>
<!-- /.container -->
</nav>
<!-- Page Content -->
<!--Nav Bar -->
<div class="row">
<div class="col-xs-1 col-sm-1 hidden-md hidden-lg">
</div>
<!-- ONE COLUMN OF SPACE -->
<nav class="hidden-xs hidden-sm col-md-2 col-lg-2 bs-docs-sidebar">
<ul id="sidebar" class="nav nav-stacked fixed manFlowSolverSV2"> <!--Nav Bar -->
<h3>Simulation Guide</h3>
<li><a href="#solverSec1">Introduction</a></li>
<li><a href="#solverSec2">Overview</a></li>
<li><a href="#solverSec4">Boundary Conditions</a></li>
<li><a href="#solverSec3">Presolver</a>
<ul class="nav nav-stacked">
<li><a href="#solverSec31">Prerequisite Files</a></li>
<li><a href="#solverSec32">Inflow BC</a></li>
<li><a href="#solverSec33">Running Script</a></li>
</ul>
</li>
<li><a href="#solverSec6">Flowsolver</a>
<ul class="nav nav-stacked">
<li><a href="#solverSec61">Prerequisite Files</a></li>
<li><a href="#solverSec62">solver.inp</a></li>
<li><a href="#solverSec63">Running Simulation</a></li>
</ul>
</li>
<li><a href="#solverSec9">Postprocessing</a></li>
<li><a href="#solverSec8">Examples</a>
<ul class="nav nav-stacked">
<li><a href="#solverSec81">Example 2</a></li>
<li><a href="#solverSec82">Example 3</a></li>
<li><a href="#solverSec83">Example 4</a></li>
</ul>
</li>
<li><a href="#solverSec7">Appendix</a>
<ul class="nav nav-stacked">
<li><a href="#solverSec71">svPre Commands</a></li>
<li><a href="#solverSec72">solver.inp Guide</a></li>
<li><a href="#solverSec73">Impedance BC Files</a></li>
<li><a href="#solverSec74">RCR BC File</a></li>
<li><a href="#solverSec75">COR BC File</a></li>
</ul>
</li>
</ul>
</nav>
<!--Main Content -->
<div class="col-xs-10 col-sm-10 col-md-9 col-lg-9" id="manualContent">
<!-- ACTUAL CONTENT -->
<div class="manFlowSolverSV2"> <section id="solverSec1" class="group"><h2>Introduction</h2>
<p><strong>Simulation Module in SimVascular</strong> solves the three-dimensional incompressible Navier-Stokes equations in an arbitrary domain, generally a vascular model reconstructed from image data. This is a complex subject with extensive underlying theory, and therefore this document will focus mainly on the practical aspects of simulation and analysis. This module includes three parts: <strong>Presolver(svPre), Flowsolver(svSolver), Postporcessing(svPost).</strong></p>
<p>The <strong>svSolver</strong> evolved from the academic finite element code PHASTA (Parallel, Hierarchical, Adaptive, Stabilized, Transient Analysis), developed mainly at RPI (Rensselaer Polytechnic Institute, Troy, NY) by Professor Kenneth Jansen. This code was in turned inspired by the Stabilized Finite Element theory developed by Professor Thomas J.R. Hughes during his Stanford years.</p>
<p>Building on the original PHASTA code, there have been a number of important additions and modifications. Professor Charles Taylor’s group at Stanford University developed key additions in the areas of Boundary Conditions and Fluid-Solid Interaction (FSI) coupling. These additions are crucial to represent with a high level of realism the way blood flows in arteries, since this flow is highly dependent on the characteristics of the vascular trees that are downstream of our three-dimensional model, and the compliance of the three-dimensional vascular tree.</p>
<h3>What’s New?</h3>
<p>Building on the above features, the Marsden lab at UCSD has added additional key functionality enabling efficient and stable solutions with models of the circulatory physiology:</p>
<ul>
<li><p><a href="/docsRefs.html#refSec2"><strong>Backflow stabilization.</strong></a> This problem usually arises in large vessels that are exposed to backflow in 3D and 2D flow simulations. This phenomenon may be a cause of divergence of the numerical scheme due to bulk reversal of the flow through an outlet, localized areas of flow reversal or use of a boundary 0D circulation model. </p></li>
<li><p>Custom and efficient <a href="/docsRefs.html#refSec3"><strong>linear solver.</strong></a> Accurate simulation of blood flow in vessels require the repeated solution of linear systems of equations with millions of unknowns. Moreover, use of closed-loop boundary models significantly increases the degree of coupling between boundary degrees of freedoms. The <strong>svLS</strong> linear solver is designed to efficiently handle large cardiovascular simulations with arbitrary boundary conditions and reduce solution times. </p></li>
<li><p>Multiscale Coupling for <a href="/docsRefs.html#refSec2"><strong>closed loop boundary conditions.</strong></a> Coupling a three-dimensional finite element solver with a 0D lumped circulation model drastically improves the possibility of realistically simulate patient-specific hemodynamics and phisiology.</p></li>
</ul>
<h3>About this guide</h3>
<p>This document will teach you the fundamentals of:</p>
<ol>
<li><strong>svPre:</strong> Preparing the necessary svSolver input files</li>
<li><strong>svSolver:</strong> Running a flow analysis </li>
<li><strong>svPost:</strong> Looking and providing interpretation to the results generated by the code </li>
</ol>
<p>In addition, this tutorial will show you a number of good practices that are important to observe during the simulation process. We will do this considering very simple geometry (a straight cylinder) to illustrate different points in a simple way.</p>
<h3>Theory and Implementation</h3>
<p>The theory and implementation details are not covered in this document. For more information about those details, please refer to our <a href="/docsRefs.html">publications page</a>.</p>
</section>
<section id="solverSec2" class="group"><h2>Overview</h2>
<h3>Process Flow of SimVascular Simulation</h3>
<p>The following figure contains a schematic representation of the processes involved in running a simulation using SimVascular.</p>
<figure>
<img class="svImg svImgLg" src="archives/sv2/flowsolver/imgs/Fig_01.png">
<figcaption class="svCaption" >Workflow for generating hemodynamic results of a cylindrical model starting from a stereolithography of its exterior surface</figcaption>
</figure>
<p>We start off with the files coming from the <a href="/docsMeshing.html">meshing</a> of the analysis: these files contain nodal and connectivity information for the finite element mesh, located in the the <em>mesh-complete/mesh-surfaces/</em> folder.</p>
<p>We then run <strong>Presolver(svPre)</strong> using the *.svpre_ file. The *.svpre file contains the set of instructions that define the boundary conditions, initial conditions, and geometrical configuration of our problem. The output of <strong>svPre</strong> is a set of files (<strong>bct.dat, geombc.dat.1, restart.0.1, numstart.dat</strong>) that are ready to be processed by <strong>svSolver</strong> to run a blood flow analysis. Running svSolver also need <strong>solver.inp</strong>, which provide further info for flowsolver.</p>
<p>Once the analysis is finished, the solver outputs files that characterize the finite element solution over a defined time period, typically several cardiac cycles. These files are taken by <strong>svPost</strong> to generate visualization files (typically *.vtu and *.vtp files) that are used to post-process and analyze the desired hemodynamic results. </p>
<p>In the following sections the components of this flow chart will be discussed in detail.</p>
<h3>Units in Simulation</h3>
<p><strong>svSolver</strong>, just like many other Finite Element Programs, does not enforce a consistent set of physical units in the analysis, but it is up to the analyst to make sure that input data are dimensionally consistent.</p>
<p>To have a consistent set of units, users are advised to either work in CGS, MGS, or SI units; see the following table. </p>
<table class="table table-bordered">
<thead>
<tr>
<th>Quantity</th>
<th>CGS Unit</th>
<th>MGS Unit</th>
<th>SI Unit</th>
</tr>
</thead>
<tr>
<td>Length</td>
<td>cm</td>
<td>mm</td>
<td>m</td>
</tr>
<tr>
<td>Mass</td>
<td>gr</td>
<td>gr</td>
<td>Kg</td>
</tr>
<tr>
<td>Time</td>
<td>s</td>
<td>s</td>
<td>s</td>
</tr>
</table>
<h3>Useful constants in Simulation</h3>
<p>The following table gathers several important physical constants of blood given in different unit
systems.</p>
<table class="table table-bordered">
<thead>
<tr>
<th>Property</th>
<th>CGS Unit</th>
<th>MGS Unit</th>
<th>SI Unit</th>
</tr>
</thead>
<tr>
<td>Dynamic viscosity $\mu$ [M· L -1 · T -1 ]</td>
<td>0.04 poise [gr· cm -1 · s -1 ]</td>
<td>0.004 [gr· mm -1 · s -1 ]</td>
<td>0.004 [Pa· s -1 ]</td>
</tr>
<tr>
<td>Blood density $\rho$ [M· L -3 ]</td>
<td>1.06 [gr· cm -3 ]</td>
<td>0.00106 [gr· mm -3 ] </td>
<td>1060 [Kg· m -3 ]</td>
</tr>
</table>
</section>
<section id="solverSec4" class="group"><h2>Boundary Condition Specification: the Physical Side of the Problem</h2>
<p>Boundary conditions are crucial to obtaining high quality cardiovascular simulation results. It is essential that boundary conditions accurately capture the physiology of vascular networks outside of the 3D domain of the model. <strong>SimVascular</strong> provides several different options for boundary condition assignment at inlets and outlets that are described in this section. Typically, we begin with some physiologic information about our problem, for instance:</p>
<ul>
<li>Flow rate info coming from <strong>MRI</strong> or <strong>ultrasound</strong> measurements.</li>
<li>Pressure values in the model obtained with a <strong>catheter</strong> transducer or a pressure cuff.</li>
<li>Vessel wall elastic properties (if we are modeling the vessel walls as deformable).</li>
</ul>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/Fig_15.png">
<figcaption class="svCaption" >Inflow, outflow and wall characterization in a simple cylindrical vessel</figcaption>
</figure>
<p>From a conceptual standpoint, no matter how geometrically complex a vascular model is (we’ll refer to this as $\Omega$), its boundaries can be classified into three groups (see figure above):</p>
<ul>
<li>An <strong>inflow</strong> boundary $\Gamma_g$. This is the set of face(s) of the model where we will prescribe a flow wave as obtained from a clinical measurement.</li>
<li>A vessel <strong>wall</strong> boundary $\Gamma_s$. This boundary represents the interface between the fluid domain and the vessel wall. In the physical world, this boundary is lined by a layer of endothelial cells, the treatment of which can be complex. Many blood flow simulations have traditionally used a <strong>rigid wall assumption</strong>. Under these circumstances, a zero velocity condition is applied on these surfaces. <strong>SimVascular</strong> also offers options for fluid structure interaction as discussed below.</li>
<li>An <strong>outflow</strong> boundary $\Gamma_h$. On this boundary, we will typically prescribe a pressure value that is uniform over the face (i.e. spatially not temporally constant) in a <strong>weak manner</strong>. A <strong>weakly applied</strong> pressure means that we are not prescribing nodal values of pressure at the nodes of the outlet face as Dirichlet boundary conditions. Instead, we apply this pressure by enforcing that the integral of the pressure field on that face must be a certain value.</li>
</ul>
<p>These boundaries have an absolutely critical impact on the numerical simulation results. The SimVascular package contains several options for boundary condition assignment. All of these use a weakly prescribed pressure formulation, with the purpose of accounting for effects of the downstream vasculature on the vascular model (see figure below). These boundary conditions include:</p>
<ul>
<li><p><strong>Resistance Boundary condition</strong>. Here, the condition prescribed on the face is a relationship between flow and pressure of the form
$p = p_0 + R\,Q$, where $R$ is the resistance parameter that characterizes the downstream vasculature, $p$ is the weakly prescribed pressure, $Q$ is the flow rate passing through the face and $p_0$ is a “flag” that sets the boundary as a “weakly-prescribed pressure boundary”. This flag has a “zero” numerical value, so the total value of the pressure on that face is simply given by $R\,Q$.</p></li>
<li><p><strong>Impedance Boundary condition</strong>. Here, the condition prescribed on the face is a relationship of the form:</p></li>
</ul>
<p>$$
p(t)=p_0 + \frac{1}{T}\,\int_{t-T}^{t} Z(t−\tau) \, Q(\tau) \, d\tau
$$</p>
<p>where $Z$ is the Inverse Fourier Transform of an impedance function that characterizes the downstream vasculature, $p$ is the weakly prescribed pressure, $Q$ is the flow rate passing through the face, and $T$ is the cardiac cycle.</p>
<figure>
<img class="svImg svImgLg" src="archives/sv2/flowsolver/imgs/Fig_16.png">
<figcaption class="svCaption" >Cardiovascular model with various boundary conditions</figcaption>
</figure>
<figure>
<img class="svImg svImgLg" src="archives/sv2/flowsolver/imgs/Fig_17.png">
<figcaption class="svCaption" >Frequency content of impedance from the literature for the left external iliac and a canine femoral artery</figcaption>
</figure>
<p>The mathematical definition of an impedance function is:</p>
<p>$$
Z(\omega)=P(\omega)\,Q(\omega),\,\omega=0,1,2,\dots
$$</p>
<p>That is, a relationship between pressure and flow modes for different frequencies. The figure above shows the typical shape of these impedances as a function of the frequency in the human iliac artery under rest and exercise conditions. Getting a good characterization of these impedance functions is once again critical to accurately represent the way blood flows in our computational domain.</p>
<ul>
<li><strong>RCR lumped parameters condition</strong>. In this boundary condition type, we use a <em>reduced-order</em> model of the downstream vasculature, considering an electric circuit analog. In this theory, the behavior of the vasculature is represented by three parameters: a proximal resistance $R_p$, a capacitance $C$, and a distal resistance $R_d$.</li>
</ul>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/Fig_18.png">
<figcaption class="svCaption" >Circuit representation of RCR block</figcaption>
</figure>
<ul>
<li><p><strong>Coronary boundary conditions</strong>. These conditions simulate what happens at the coronary outlets. The implementation in the <strong>svSolver</strong> follows the derivations in <a href="/docsRefs.html#refSec2">this paper</a>.</p></li>
<li><p><strong>Closed-loop boundary circulation model</strong>. The capability of coupling a 3D finite element model with a lumped parameter model is built into the <strong>svSolver</strong>. Documentation on this feature will be available with later releases of the code. </p></li>
</ul>
<h3>Boundary conditions considered in Example 1</h3>
<p>Before we move on, let us recap the type of <em>physical problem</em> (<strong>Example 1</strong>) we are solving: the geometry used in this problem consists of an idealized blood vessel, represented by a cylindrical segment with a radius $r=2$ cm and length $L=30$ cm. We prescribe an idealized constant inlet volumetric flow rate $Q$ of $100$ cc/sec to a parabolic profile at the inlet face of the model ($\Gamma_g$). The dynamic viscosity $\mu$ and density $\rho$ of the blood are 0.04 poise and 1.06 gr/cm$^3$, respectively. The lateral surface of the vessel ($\Gamma_{s}$) is considered to be rigid (therefore, we will apply a zero-velocity condition there). For the outlet boundary ($\Gamma_h$), a spatially-constant pressure boundary condition is weakly enforced via a resistance $R$.
In this problem, we will consider a resistance of $R = 1333.0$ dynes·s/cm$^5$. </p>
<p>This resistance will give a (weakly-applied) pressure at the outlet face of</p>
<p>$$
p=p_0 + R\,Q = 0.0 + 1333.0 \cdot 100.0=133300.00 \approx 100\,\text{mmHg}
$$</p>
<p>(recall that $1.0$ mmHg = $1333.2$ dyn/cm$^2$). For steady flow in a long tube with a circular cross section, the flow will develop a flow profile known as the <em>Poiseuille</em> flow profile assuming the flow remains laminar. The flow will remain laminar in a circular tube assuming that the non-dimensional parameter given by the <em>Reynolds</em> number $Re$ is $Re < 2100$.</p>
<p>The definition of the Reynolds number is given by:</p>
<p>$$
Re = \frac{\rho\,D\,V}{\mu}
$$</p>
<p>where $V$ is a representative velocity of the flow, $D$ is a characteristic dimension of the vessel where the flow is happening (in this case, the diameter), and $\mu$ and $\rho$ are the dynamic viscosity and density, respectively.</p>
<p>For this problem, the Reynolds number is about $884$, so it is well within the laminar flow regime. For a fully developed flow, the axisymmetric profile is parabolic and is given by: </p>
<p>$$
v_z(a) = v_z^{max}\left[1-\left(\frac{a}{r}\right)^2\right]
$$</p>
<p>where $v_z^{max}$ is the maximum velocity at the center of the vessel, a is the radial coordinate from center of the vessel $0\le a \le r$ . The expression for maximum velocity is given by:</p>
<p>$$
v_z^{max} = 2\frac{Q}{\pi\,r^2}
$$</p>
<p>where $Q$ is the volumetric flow rate. The wall shear stress $\tau$, is given by</p>
<p>$$
\tau = \frac{2\,\mu\,v_z^{max}}{r}
$$</p>
<p>For this example, the maximum velocity is $v_z^{max} = 16.68$ cm/s and the wall shear stress is $\tau$ = $0.67$ dynes/cm$^2$.</p>
</section>
<section id="solverSec3" class="group"><h2>Presolver(svPre)</h2>
<p><strong>svPre</strong> is the preprocessor for <strong>svSolver</strong>. The input files to svPre contain a complete description of the discrete model: nodal coordinates, element connectivity, element adjacency information and connectivity of boundary nodes and elements. Running <strong>svPre</strong> with an input *.svpre file with the appropriate commands, generates the input files for <strong>svSolver</strong>.</p>
<p>The <strong>svPre</strong> program is called either from the command line (in terminal) or using the <strong>SimVascular</strong> GUI. The input files for <strong>svPresolver</strong> are those generated by <strong>Simvascular Meshing Module</strong>. We will review this process briefly with a simple example of steady flow through a cylinder (Example 1). Before we start, first set the project folder as the example folder (…/cylinder)</p>
</section>
<section id="solverSec31" class="subgroup"><h3>Prerequisite Files for svPre</h3>
<p>These prerequisite files for svPre are generate by the output from Meshing Module (Click <em>Write Files</em> button in Mesh tab after meshing).</p>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/meshfiles.png">
<figcaption class="svCaption" >Folder structure and file created after clicking on <b>Write Files</b></figcaption>
</figure>
<p>These files are:</p>
<p>in the <strong>mesh-complete/</strong> folder: </p>
<ul>
<li><strong>mesh-complete.mesh.vtu</strong>, vtu file containing the solid mesh generated with TetGen.</li>
<li><strong>mesh-complete.exterior.vtp</strong>, vtp file containig all the exterior elements of the mesh generated with TetGen.</li>
<li><strong>walls_combined.vtp</strong>, vtp file containing all surface elements assigned to the wall, possibily combined from various surfaces. </li>
</ul>
<p>in the <strong>mesh-complete/mesh-surfaces/</strong> folder:</p>
<ul>
<li><strong>inflow.vtp</strong>, vtp file containing the meshed inlet surface.</li>
<li><strong>outlet.vtp</strong>, vtp file containing the meshed outlet surface.</li>
<li><strong>wall.vtp</strong>, vtp file containing the meshed wall surface.</li>
</ul>
<p>The files for Example 1 can be found <a href="documentation/flowsolver/files/examples.zip">here</a>. Create an empty folder on your hard drive to unzip the content of the archive. The following files are contained:</p>
<p><strong>HINT</strong> - It is advisable that you set the project folder as <em>cylinder</em>. <strong>SimVascular</strong> will use this folder name as the default when creating new files. Using a meaningful folder name will make sure that your model files are named consistently. Also store the files containing the inlet flow rates in a folder called <em>flow-files</em>. Your problem may have more that one inflow wave form file. In this case, we only have a single flow file (called inflow.flow).</p>
<p>The format of the <strong>steady.flow</strong> file is as follows:</p>
<pre class="highlight plaintext"><code># Time (sec) Flow (cc/sec)
0 -100.0
1 -100.0
</code></pre>
<p>The first line is a comment line that you don’t need to include, but it helps to remember what units you used in the analysis. Then, two columns of numbers follow. The first column contains time values, and the second column flow values.</p>
<p><strong>WARNING</strong>: please note that flow coming <strong>into</strong> the model (forward flow) will have a negative sign, (like in the example considered here), whereas flow coming <strong>out of</strong> the model (backflow) will be positive. A good way to remember that is that in the case of forward flow, the vector that gives you the direction of the flow and the normal to the face of the model point in opposite directions, and therefore their dot product will be negative.</p>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/Fig_04.png">
<figcaption class="svCaption" >Cylinder with negative inflow</figcaption>
</figure>
<p>On the other hand, in a situation of back flow, the numerical value in the *.flow file with be positive. </p>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/Fig_05.png">
<figcaption class="svCaption" >Cylinder with positive inflow</figcaption>
</figure>
<p>In this problem, since we are running a steady case, our physical time goes from 0.0 to 1.0 seconds, and the flow is constant with a value of 100.0 cc/sec.</p>
<p><strong>HINT</strong>: it is very important that you are absolutely sure about the physical dimensions of your model: every unit (length, time, flow, density, etc.) in your analysis must be dimensionally consistent. You can easily check the size of your model in <a href="http://www.paraview.org/">Paraview</a> before importing it into <strong>SimVascular</strong>.</p>
<p>In this case, our cylinder has a radius $r=2.0$ cm and length $L=30$ cm.</p>
</section>
<section id="solverSec32" class="subgroup"><h3>Inlet Boundary Condition Specification (From GUI)</h3>
<p>First we need to create a file called <strong>bct.dat</strong> (its vtp format file is also created, called <strong>bct.vtp</strong>) that specifies the boundary conditions at the inlet defined by inflow.vtp</p>
<p>In the <strong>SimVascular</strong> GUI window, go to the <em>Inflow BC</em> subtab under <em>Simulations</em>. You will have to enter the following values in the various boxes/buttons of the GUI (see figure below):</p>
<figure>
<img class="svImg scImgLg" src="archives/sv2/flowsolver/imgs/BCT_Creation.png">
<figcaption class="svCaption" >Creating a <b>bct.dat</b> file through the GUI</figcaption>
</figure>
<ul>
<li><p>Under <strong>Analytic Shape of Profile</strong>, select the <strong>parabolic</strong> radio button. This options allows a parabolic velocity profile to be applied at the inlet. Other options are <strong>plug</strong>, which applies a constant velocity profile throughout the inlet face and Womerseley, that uses a closed form solution for the velocity profile of pulsatile flow in arteries. </p></li>
<li><p>In the <strong>Mesh Face File (vtp)</strong> box, click on the <strong>Browse</strong> button and select the <strong>inflow.vtp</strong> file in the <em>mesh-complete/mesh-surfaces/</em> folder.</p></li>
<li><p>In the <strong>Flow Rate File</strong> box, click on <strong>Browse</strong> and select the desired *.flow file. In this case, the file is <strong>steady.flow</strong>.</p></li>
</ul>
<p>Under the <strong>Parameters</strong> menu, enter the following values:</p>
<ul>
<li><p><strong>Period</strong>: $1.0$ sec. For a steady flow problem like ours, the concept of <em>period</em> is somewhat vague. In this case, $1.0$ means the amount of physical time that we are going to run our simulation for.</p></li>
<li><p><strong>viscosity</strong>: $0.04$ Poise (gr/cm/s).</p></li>
</ul>
<p><strong>WARNING</strong>: Be very careful with the units! This is one the points where it is easy to make a mistake with the dimensions. For consistency, besides entering the right numerical value, try to also modify the units listed next to it (see figure below).</p>
<ul>
<li><strong>density</strong>: $1.06$ gr/cm$^3$. Same as before, be very careful with these units!</li>
</ul>
<p>Under the <strong>Output Parameters</strong> menu, enter the following values:</p>
<ul>
<li><p><strong>num of periods</strong>: always enter 1 here. If you want to run your simulation for multiple cardiac cycles, it is possible to ask <strong>svSolver</strong> to loop over the information given by the <strong>bct.dat</strong> file for just one cycle (using the <strong>solver.inp</strong> file). By doing this, you will minimize the size of this file that can potentially be very large. </p></li>
<li><p><strong>num pts in period</strong>: 2. This is the number of <em>temporal</em> data points that you want to have in your bct.dat. This is not necessarily the number of time points in the *.flow file. In this case, they match (2 in both cases), but this is because this is a very simple example using steady flow and two time points is all we need to characterize a constant flow. In general, your *.flow file will have in the order of $20$ data points over the cardiac cycle (that’s about how many points you will be able to reconstruct using <strong>PC-MRI</strong>, for example), and your bct.dat will have on the order of $100$-$200$ points. Whatever is enough to have a smooth representation of the inflow wave mapping to velocity vectors at the inlet face.</p></li>
<li><p><strong>num fourier modes</strong>: 1. Fourier smoothing allows to smooth your inlet flow curve and to make sure that you have a periodic function in the specified interval. </p></li>
</ul>
<p><strong>WARNING</strong>: Be careful with this! <strong>SimVascular</strong> is doing a Fourier Series approximation of the data that you provide in the *.flow file. Since in this case, our data is constant flow, we only need one Fourier mode to capture this appropriately. For pulsatile flow problems, we will need more Fourier Modes to accurately represent the *.flow data (usually, $10$ Fourier modes is enough for a pulsatile problem).</p>
<p>After you are done entering all these parameters, click on the <em>CREATE 3-D FLOW SOLVER BC FILE</em> button to generate the <strong>bct.dat</strong> file. The format of this file is as follows:</p>
<pre class="highlight plaintext"><code>np nl
x1 y1 z1 nl nn
vx1 vy1 vz1 t1
. . . .
. . . .
. . . .
vxnl vynl vznl tnl
. . . .
. . . .
. . . .
xnp ynp znp nl
vx1 vy1 vz1 t1
. . . .
. . . .
. . . .
</code></pre>
<p>The file defines the inflow boundary condition both spatially and in time. The spatial definition is obtained through $n_p$ point-wise velocity input blocks. In this case, $n_p = 102$, this is the total number of nodes on the <strong>inflow.vtp</strong> face. The temporal definition is given by $n_l$ input lines of the values at a certain position for $n_l$ time points, $t_1$ to $t_{n_l}$. In this case, $n_l = 2$ points (this is the value we entered in the <em>num pts in period</em> box. </p>
<p>Each block of data has, for each of the $n_p = 102$ spatial points, the following info:</p>
<ul>
<li><p>The coordinates of the point: $x_1$ $x_2$ $x_3$ and the number of time points $n_l$.</p></li>
<li><p>The list of velocity vectors $v_x$ $v_y$ $v_z$ at time t.</p></li>
</ul>
<p>A vtp file <strong>bct.vtp</strong> can be written using this option <strong>Create Vtp</strong> to graphically visualize the velocity distribution at the inlet surface with Paraview, as shown in the picture below.</p>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/BCT_Cration_VTP.png">
<figcaption class="svCaption" >Visualizing the inlet velocity profile in Paraview</figcaption>
</figure>
</section>
<section id="solverSec33" class="subgroup"><h3>Running Script</h3>
<p>To define the initial and boundary conditions of this problem, svPre need a script file (*.svpre) file. We go to the <em>Simulations->Create 3-D Solver Files</em> tab. In the “Create PreSolver script file” menu (see figure below), make sure that the right *.svpre file is loaded in the box (in this case, it should be cylinder.svpre . Click on the “<strong>Load PreSolver scriptfile</strong>” button. The following screen will appear:</p>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/svpre_gui.png">
<figcaption class="svCaption" >Running <b>svPre</b> through the GUI</figcaption>
</figure>
<p>The contents of the cylinder.svpre script file are:</p>
<pre class="highlight plaintext"><code># Read Mesh info
mesh_and_adjncy_vtu mesh-complete/mesh-complete.mesh.vtu
# Assign IDs to the surfaces
set_surface_id_vtp mesh-complete/mesh-complete.exterior.vtp 1
set_surface_id_vtp mesh-complete/mesh-surfaces/inflow.vtp 2
set_surface_id_vtp mesh-complete/mesh-surfaces/outlet.vtp 3
# Set Inlet BC
prescribed_velocities_vtp mesh-complete/mesh-surfaces/inflow.vtp
# Set BCT for Inlet
fluid_density 1.06
fluid_viscosity 0.04
bct_analytical_shape parabolic
bct_period 1.0
bct_point_number 2
bct_fourier_mode_number 1
bct_create mesh-complete/mesh-surfaces/inflow.vtp flow-files/steady.flow
bct_write_dat
bct_write_vtp
# Set Outlet BC
zero_pressure_vtp mesh-complete/mesh-surfaces/outlet.vtp
# Set Wall BC
noslip_vtp mesh-complete/walls_combined.vtp
# Write geometry and property data to geombc.dat.1
write_geombc
# Write initial values of velocity, pressure, etc to restart.0.1
write_restart
</code></pre>
<p>As we said before, each line of this *.svpre file represents a command that will be executed by <strong>svPre</strong>. This file needs to be edited to incorporate the right parameters/conditions for this problem. A complete list of svpre commands available is <a href="#solverSec71">this section</a>.Here is a description of each line.</p>
<p>The first line is used to define the topology of the finite element mesh. The first command <strong>mesh_and_adjncy_vtu</strong> used a vtk unstructured mesh file to define node coordinates, element connectivities and an adjacency relationship between elements. </p>
<pre class="highlight plaintext"><code>mesh_and_adjncy_vtu mesh-complete/mesh-complete.mesh.vtu
</code></pre>
<p>The following command is used to assign an ID to all model surfaces.</p>
<pre class="highlight plaintext"><code>set_surface_id_vtp cylinder.exterior.vtp 1
</code></pre>
<p><strong>HINT</strong>: This line tags all the exterior faces in the model with an identifier (a Suface ID) , in this case, the number one. The assignment can tell flowsolver later which faces needed for wall stress calculation. We also need to introduce a new command if we want to activate the resistance boundary condition at the outlet face. We had previously determined that a resistance equal to $R = 1333.0$ dynes$\,$s/cm$^5$ needs to be applied at that outlet. </p>
<p>In order to do this, we need to assign a Surface ID that will help us later to identify the face and assign the correct resistance value. This is a trivial case, because we only have one single outflow face, and therefore one single resistance. But imagine one case where many more are needed. In this case, it is very important to meticulously label all the outlet faces with a meaningful name, and to make a sketch that helps you remember the list of Surface IDs that you considered in the model. Each surface ID will have a corresponding Resistance value (or impedance function, or set of RCR parameters etc.).</p>
<p>Going back to our *.svpre file, we need to add a line specifying the surface ID for the outlet face. We also number all other model surfaces in case we need to apply different boundary conditions through the <strong>solver.inp</strong> file.</p>
<pre class="highlight plaintext"><code>set_surface_id_vtp mesh-complete/mesh-surfaces/inflow.vtp 2
set_surface_id_vtp mesh-complete/mesh-surfaces/outlet.vtp 3
</code></pre>
<p><strong>HINT</strong>: Since a face can only have one ID, now ID 1 just represents cylnder wall because the inlet and outlet are assigned with ID 2 and 3.</p>
<p>The following line uses the existing <strong>inflow.vtp</strong> file to define a boundary subregion with applied velocities. </p>
<pre class="highlight plaintext"><code>prescribed_velocities_vtp mesh-complete/mesh-surfaces/inflow.vtp
</code></pre>
<p>Note that we are just pointing to the right file (inflow.vtp) in the mesh-surfaces folder where we want some velocity vectors to be prescribed. These velocity vectors are given by the <strong>bct.dat</strong> file, already created from GUI as shown above. We must use the command prescribed<em>velocities</em>vtp in each surface where we prescribe some flow wave via a Dirichlet condition on the velocity vectors of that face. Instead, we can also create bct.dat (and bct.vtp) here using script commands as below. Similar to the GUI example above, these commands need to provide info about fluid density, fluid viscosiyt, velocity profile shape, period length, number of points in one period, number of fourier modes, inlet face, and flow file.</p>
<pre class="highlight plaintext"><code>fluid_density 1.06
fluid_viscosity 0.04
bct_analytical_shape parabolic
bct_period 1.0
bct_point_number 2
bct_fourier_mode_number 1
bct_create mesh-complete/mesh-surfaces/inflow.vtp flow-files/steady.flow
bct_write_dat
bct_write_vtp
</code></pre>
<p>Like before, we are only pointing to the right path of the surface file where we want to prescribe the nonslip (rigid wall) condition. This is also a Dirichlet condition that makes all the velocity vectors of the nodes of the surface <strong>wall</strong> to be zero.</p>
<pre class="highlight plaintext"><code>zero_pressure_vtp mesh-complete/mesh-surfaces/outlet.vtp
</code></pre>
<p>By using this condition, we are making the face <strong>outlet</strong> into a <strong>weakly-pressure</strong> face. This is mathematically expressed by the expressions we saw before:</p>
<p>$$
p = p_0 + R\,Q
$$</p>
<p>$$
p(t)=p_0 + \frac{1}{T}\,\int_{t-T}^{t} Z(t−\tau) \, Q(\tau) \, d\tau
$$</p>
<p>This expression sets $p_0 = 0$ for the face <strong>outlet</strong>. The total pressure set on that face will be the result of the flow-pressure operator considered (i.e., resistance, impedance, RCR, Coronary etc.).</p>
<p>The following line set no-slip boundary conditions for all walls. Since for this simple cylinder, there is only one wall and We can also use mesh-complete/mesh-surfaces/wall.vtp instead of mesh-complete/walls_combined.vtp.</p>
<pre class="highlight plaintext"><code>noslip_vtp mesh-complete/walls_combined.vtp
</code></pre>
<p>The last two lines are the culmination of all the work you have been doing in <strong>SimVascular</strong> to this point!</p>
<pre class="highlight plaintext"><code>write_geombc
write_restart
</code></pre>
<p>They generate two binary files (geombc.dat.1) and (restart.0.1) that are used as inputs for <strong>svSolver</strong> and are used to run the flow analysis. </p>
<ul>
<li><p><strong>geombc.dat.1</strong> contains the combination of geometry, material properties and boundary conditions specified in the problem.</p></li>
<li><p><strong>restart.0.1</strong> contains the set of initial conditions for our problem. We haven’t said anything about initial conditions till now. If you do not do something different, <strong>SimVascular</strong> will specify an almost-zero velocity initial condition for all the nodes of the mesh and a zero pressure <strong>initial condition</strong>. Here, the number <strong>0</strong> refers to <strong>time step zero</strong>, as it corresponds to the first file of the simulation.</p></li>
</ul>
<p>Now, click on <strong>Run PreSolver</strong>. This command will actually launch the <strong>svPre</strong> application. A window will pop up, and you will see the list of commands of your <em>.svpre file being executed. After a few seconds (or minutes, depending on the size of the problem), the files *</em>geombc.dat.1** and <strong>restart.0.1</strong> will be generated.</p>
<p>You can do the same if, instead of using the <strong>SimVascular</strong> GUI, you edit the *.svpre file like shown above, and then, from the command line, type:</p>
<pre class="highlight plaintext"><code>%svpre cylinder.svpre
</code></pre>
<p><strong>HINT</strong>: In both files, the number “.1” refers to the <strong>partition number</strong> of the file. Remember <strong>svSolver</strong> has the ability of running a problem <em>in parallel</em> (i.e., using multiple processors or computing cores), using MPI (message-passing interface). When we run a job using multiple processors, the first thing that happens is the “splitting” of these two files into as many processors we are going to use in our analysis, so the calculations can be performed faster. For example, if we use $4$ processors later in svSolver, these files will be split as follows:</p>
<pre class="highlight plaintext"><code>geombc.dat.1 => geombc.dat.1 , geombc.dat.2 , geombc.dat.3 , geombc.dat.4
restart.0.1 => restart.0.1 , restart.0.2 , restart.0.3 , restart.0.4
</code></pre>
<p>Roughly speaking, each of the four files is $1⁄4$ of the size of the original un-split file. For a generic time step <strong>n</strong>, the solution will be given by the following files:</p>
<pre class="highlight plaintext"><code>restart.n.1 , restart.n.2 , restart.n.3 , restart.n.4 , ...
</code></pre>
<p>We are almost done. There is only one thing left in the <strong>svPresolver</strong> part: to generate the numstart.dat dat. To do this, click on the <strong>Create File</strong> button under <strong>Analysis Files</strong>. </p>
<p>This file is really simple: it contains the scalar <strong>0</strong>. This number is used by the solver to identify the restart file that should be used as initial condition. In this case, since this file is <strong>restart.0.1</strong>, the file <strong>numstart.dat</strong> should contain a <strong>0</strong>. If for whatever reason, the initial file of our simulation was <strong>restart.300.1</strong>, the numstart.dat file should have a <strong>300</strong> entry. The value stored in this file gets updated as the simulation advances in time (we will see this later one).</p>
<h4>Final recap of the files generated by <strong>svPre</strong></h4>
<p>At this point, we are almost ready to run the flow solver. Using <strong>svPre</strong>, we have generated the following files that are direct inputs to the solver:</p>
<ul>
<li><strong>geombc.dat.1</strong> : this file contains the combination of geometry and boundary conditions specified in the problem.</li>
<li><strong>restart.0.1</strong> : this file contains the set of initial conditions for our problem. </li>
<li><strong>numstart.dat</strong>: this file contains the scalar <strong>0</strong>. This number is used by the solver to identify the restart file that should be used as initial condition.</li>
<li><strong>bct.dat</strong> : this file contains the history of velocity vectors at the inflow face of the model according to a prescribed flow wave coming from a *.flow file. </li>
</ul>
</section>
<section id="solverSec6" class="group"><h2>Flowsolver (svSolver)</h2>
<p><strong>svSolver</strong> is the flowsolver for SimVascular simulation. The <strong>svSolver</strong> program is called either using the <strong>SimVascular</strong> GUI or from the command line (in terminal). The input files to svSolver contain a complete description of model geometry,material properties, initial condition, boundary condition, and various parameter to control simulation. We will review this process briefly by keeping using the above cylinder example. First make sure the project folder is the example folder (…/cylinder)</p>
</section>
<section id="solverSec61" class="subgroup"><h3>Prerequisite Files for svSolver</h3>
<p>Besides bct.dat, geombc.dat.1,restart.0.1 and numstart.dat, we are only missing one file in order to be able to run our analysis. This file is another input file for the solver that controls the actual flow of the numerical simulation, specifying parameters such as time step size, number of time steps, number of nonlinear iterations, boundary condition control, etc. This file needs to have the name <strong>solver.inp</strong> (input file for the solver), and we will characterize it in detail in the following section. A detailed description is also presented in <a href="#solverSec72">this section</a>.</p>
<p>These five files are the bare minimum we need to run an analysis. However, if we want to perform more complicated simulations, involving more complex boundary conditions, we will need more input files.</p>
<h4>Impedance Boundary Condition simulations :</h4>
<p>In addition to the five standard files (geombc.dat.1, restart.0.1, numstart.dat, bct.dat, solver.inp), we will need to provide impedance functions in the time domain for each impedance outlet, as well as a history of flow rates for each outlet. We will have therefore two additional ascii files: <strong>impt.dat</strong> (containing the impedance functions for each of the outlets), and <strong>Qhistor.dat</strong> (containing the flow rate history). A detailed description is <a href="#solverSec73">here</a>.</p>
<h4>RCR Boundary Condition simulations :</h4>
<p>In addition to the five standard files (geombc.dat.1, restart.0.1, numstart.dat, bct.dat, solver.inp), we will need to provide the RCR parameters in a ascii file that will set the relationship between flow and pressure on each outflow face. This is done by defining a file named <strong>rcrt.dat</strong> containing such parameters. A detailed description is <a href="#solverSec74">here</a>.</p>
<h4>Coronary Boundary Condition simulations :</h4>
<p>In addition to the five standard files (geombc.dat.1, restart.0.1, numstart.dat, bct.dat, solver.inp), we will need to provide the coronary model parameters in a ascii file that will set the relationship between flow and pressure on each outflow face. This is done by defining a file named <strong>cort.dat</strong> containing such parameters. A detailed description is <a href="#solverSec75">here</a>.</p>
<h4>Closed-loop boundary conditions:</h4>
<p>This will required an executable that implements a lumped parameter network model for the patient circulation. This will be covered in a later version of this tutorial. Stay tuned!</p>
<p>We have completed the section on preprocessing the model. Let’s move on to <strong>svSolver</strong>, define the solver.inp file and run the analysis.</p>
</section>
<section id="solverSec62" class="subgroup"><h3>solver.inp</h3>
<p>The main goal of this section is to define the file we are missing to run the analysis. This is the <strong>solver.inp</strong> file (i.e., input parameters for the solver). Most parameters are already assigned default values for cardiovascular simulation. Only a very small number of parameters must be set up in solver.inp. For this problem, the file we need will look like this:</p>
<pre class="highlight plaintext"><code># ================
# SOLUTION CONTROL
# ================
Number of Timesteps: 200
Time Step Size: 0.03
# ==============
# OUTPUT CONTROL
# ==============
Number of Timesteps between Restarts: 10
Number of Force Surfaces: 1
Surface ID's for Force Calculation: 1
# ===================
# MATERIAL PROPERTIES
# ===================
Viscosity: 0.04
Density: 1.06
# ==================================
# CARDIOVASCULAR MODELING PARAMETERS
# ==================================
Number of Coupled Surfaces: 1
Number of Resistance Surfaces: 1
List of Resistance Surfaces: 3
Resistance Values : 1333
# =============
# STEP SEQUENCE
# =============
Step Construction: 0 1 0 1
</code></pre>
<p>The file consists of a number of blocks, each block containing a number of lines that are instructions for the solver.</p>
<p><strong>WARNING</strong>: it is very important that the wording of each line is exactly as presented here. Even a slight change (for instance, a change from uppercase to lowercase) will make the solver not understand the command!</p>
<p>The lines preceded by a <strong>#</strong> sign are comments and are ignored by the solver. Likewise, anything placed after a <strong>#</strong> on a given line is also ignored.</p>
<p><br></p>
<h4>Solution Control Block</h4>
<p>In this block, the different commands are:</p>
<p><strong>Number of Timesteps: 200</strong> and <strong>Time Step Size: 0.03</strong> - These two lines control the amount of physical time that you run your problem for. In this case,</p>
<p>$$
\text{Total physical time} = \text{N. time steps} \times \text{Time Step Size} = T = N \times \Delta t = 200 \times 0.03 = 6.0\,\text{sec}
$$</p>
<p>Note that this doesn’t match the <strong>period</strong> options we specified to generate the <strong>bct.dat</strong>. In this case, like we mentioned before, it does not really make sense to talk about a <em>cardiac cycle</em> (this is a steady flow), but if we wanted to run this analysis for <em>six</em> cardiac cycles, we would have to run the problem for $6.0$ seconds of physical time. If we kept our choice of time step size the same ( $\Delta t = 0.03$ sec), we will need a total number of time steps of $N = 200$.</p>
<p><strong>WARNING</strong>: Note that this $N$ is the total number of time steps you need in your numerical simulation to model a certain physical time, given a prescribed $\Delta t$. This is not to be confused with the previous number of time steps you used to generate the bct.dat!</p>
<p><strong>WARNING</strong>: Now the question is: how do you come up with a reasonable value for $\Delta t$? There is not a straightforward answer for this. $\Delta t$ is the parameter that controls your <strong>temporal discretization</strong>, which is something that works in a similar fashion to the <strong>spatial discretization</strong> given by your mesh: the finer, the more accurate the results, but also the bigger the size of the problem and the time to solve it! We don’t want to get into a lot of theoretical details now, so we will just provide you with a reasonable <strong>recipe</strong> to evaluate this parameter. The <strong>recipe</strong> to estimate a reasonable $\Delta t$ is based on a dimensionless parameter called the <strong>CFL</strong> number. The CFL number relates the velocities happening in the fluid domain ($v$), a temporal discretization parameter ($\Delta t$), and a mesh discretization parameter (i.e. mesh size) ($h$) as follows:</p>
<p>$$
\text{CFL} = \frac{v\,\Delta t}{h}
$$</p>
<p>We want this <strong>CFL</strong> number to be around $1.0$. This will mean that, for the velocities present in our fluid domain, the temporal and spatial discretizations are <em>balanced</em>. In our problem, it can be shown that the average expected velocity is about $v = 16.7$ cm/s; the spatial discretization parameter or finite element size is $h = 0.5$. Therefore, if we shoot for a CFL number close to one, we have:</p>
<p>$$
\Delta t = \frac{h}{v} = \frac{0.5}{16.7} = 0.03
$$</p>
<p>Of course, you can imagine that in a real-world problem things are way more complicated to evaluate: it will be much harder to estimate where your model will have the largest velocities, what the mesh element size will be there, etc. The time step size $\Delta t$ is a parameter that will have a very important impact on the performance of the linear solver of equations. The smaller you make it, the <em>easier</em> you will be for the solver to find a solution, but the longer it will take you to reach a certain point in time.</p>
<p><br></p>
<h4>Material Properties Block</h4>
<p>This block contains the values for density and dynamic viscosity of blood: nothing really new here. Be careful though and make sure that you use the same units you have been using through the simulation process!</p>
<p><br></p>
<h4>Output Control Block</h4>
<p>In this block, the meaning of the command is:</p>
<p><strong>Number of Timesteps between Restarts: 10</strong> - This line tells the solver how often it should save solution files. In this problem, you are really calculating $200$ solutions to the problem at $200$ different time points, but in general you do not want to save a solution file for every single time step. Keep in mind that two consecutive solutions are only $\Delta t = 0.03$ seconds apart! In this line, we are asking the solver to save every other $20$ files. Therefore, the output files of the solver will look like this: restart.0.*, restart.10.*, restart.20.*, restart.30.*, …., restart.190.*, restart.200.*</p>
<p><strong>Number of Force Surfaces: 1</strong> - This is the number of surfaces of the model where we are calculating the wall stress.</p>
<p><strong>Surface ID’s for Force Calculation: 1</strong> - This line is list of surface ID’s considered for walls stress calculation. In our case, we only defined one surface ID (the number 1, assigned to the cylinder in svPre).</p>
<p><br></p>
<h4>Cardiovascular Modeling Parameters Block</h4>
<p>This is the block that controls the Boundary Conditions and the other features such as deformable wall parameters. The meaning of each command is:</p>
<ol>
<li><p><strong>Number of Coupled Surfaces: 1</strong> - This is the number of surfaces of the model where we are specifying a relation that couples Flow and Pressure. In order words, this number is the total number of <strong>Resistance, Impedance, RCR and coronary surfaces</strong> we have in our problem. In this case, since we only have one outlet with a resistance boundary condition, we enter a 1 in this line.</p>
<p><strong>WARNING</strong>: this line refers to coupled surfaces. A constant pressure outlet with no coupling between flow and pressure <strong>does not</strong> qualify as a coupled surface! </p></li>
<li><p><strong>Number of Resistance Surfaces: 1</strong> - This line sets the number of resistance surfaces in the model. In our case, we have one resistance surface.</p></li>
<li><p><strong>List of Resistance Surfaces: 3</strong> - This line the list of surface ID’s considered in the model for Boundary Condition specification. In our case, we only defined one surface ID (the number 3), at the outlet face of the model. It is very important that this number matches what you used in your *.svpre file. Otherwise, things will not work!</p></li>
<li><p><strong>Resistance Values : 1333.0</strong> - This line the list of resistancese considered in the outlets of the model. In our case, this resistance is 1333.0. </p></li>
</ol>
<p><strong>WARNING</strong>: Be very careful with the units! It is also very important that ordering of the resistance values in this line and the surface ID’s you provided in the previous line is consistent. This is a very common place to make a mistake. It is also very important that whatever you enter in these last two lines is consistent with want you entered in the *.svpre file. </p>
<p>Let us illustrate this with a more complex problem with 4 outlets (see figure below)</p>
<figure>
<img class="svImg svImgMd" src="archives/sv2/flowsolver/imgs/FourOutlets.png">
<figcaption class="svCaption" >Schematic representation of a model with four outlets</figcaption>
</figure>
<p>The *.svpre file should read something like this:</p>
<pre class="highlight plaintext"><code>.
.
.
zero_pressure_vtp mesh-surfaces/outlet1.vtp
zero_pressure_vtp mesh-surfaces/outlet2.vtp
zero_pressure_vtp mesh-surfaces/outlet3.vtp
zero_pressure_vtp mesh-surfaces/outlet4.vtp
#
set_surface_id_vtp exterior_faces.vtp 1
set_surface_id_vtp mesh-surfaces/outlet1.vtp 2
set_surface_id_vtp mesh-surfaces/outlet2.vtp 3
set_surface_id_vtp mesh-surfaces/outlet3.vtp 4
set_surface_id_vtp mesh-surfaces/outlet4.vtp 5
.
.
.
</code></pre>
<p>And the solver.inp file:</p>
<pre class="highlight plaintext"><code>.
.
.
Number of Resistance Surfaces: 4
List of Resistance Surfaces: 2 3 4 5
Resistance Values : 20000 10000 15000 21000
.
.
.
</code></pre>
<p><br></p>
<h4>Step sequence Block</h4>
<p>This line controls the non-linear iteration loop within the time step. The syntax of the line represents a two nonlinear iteration process for each time step. The <strong>0</strong> tells the solver to make a solve, the <strong>1</strong> to make an update of the solution. Since this sequence 0 1 is repeated, the two iterations are defined. </p>
<p><strong>WARNING</strong>: Deciding on the adequate number of non-linear iterations for a problem is also a non-trivial problem. In principle, we need to iterate until the residual (i.e., the <em>error</em>) of our numerical solution is small enough. But doing many non-linear iterations on each time step is very costly. In general, for steady flow problems, 1 or 2 non-linear iterations are enough. For pulsatile problems, at least three non-linear iterations are needed. For deformable wall problems, 4 or more non-linear iterations are required. This parameter, together with the time step size $\Delta t$ and the quality of the spatial discretization given by the finite element mesh, will completely determine the performance of the linear solver of equations. The better chosen these parameters are, the faster and more accurately our simulation will run. We will talk more about this later.</p>
<p>The set of instructions explained here constitute a very small sample of all the possible instructions the <strong>svSolver</strong> can take via a solver.inp file. A more detailed discussion can be found in <a href="#solverSec72">this section</a>.</p>
</section>
<section id="solverSec63" class="subgroup"><h3>Running Simulation</h3>
<p>At this point we have generated all the files we need for this problem. Start simulation from the GUI.</p>
<pre class="highlight plaintext"><code>Tab: Simulation -> Run Solver -> localhost
Select the project dir for Run Dir
Select Log Dir
Starting Step Number: 0
Click the button "whoami" to set username
Choose the number of processors - localhost num procs:4
Click "Run Simulation"
Wait a few seconds, Click the last "Start Trail" button to track the simulation progress
</code></pre>
<figure>
<img class="svImg scImgLg" src="archives/sv2/flowsolver/imgs/svsolver_gui.png">
<figcaption class="svCaption" >Running simulation through the GUI</figcaption>
</figure>
<p>You cana also run simulation by a command lines.</p>
<pre class="highlight plaintext"><code>%mpiexec -np 4 svsolver
</code></pre>
<p>This will launch a four-processor job in your computer. Therefore, the input file are split as follows: </p>
<pre class="highlight plaintext"><code>geombc.dat.1 => geombc.dat.1, geombc.dat.2, geombc.dat.3, geombc.dat.4
restart.0.1 => restart.0.1, restart.0.2, restart.0.3, restart.0.4
</code></pre>
<p>At the same time, the solver copies all these files to a newly created folder called <strong>4-procs_case/</strong> under the project folder, and this is where all the output files of the analysis will be written to. In general, if you launch a <strong>n</strong> processor job, all the files will be copied to a <strong>n-procs_case/</strong> folder.</p>
<p>You can check the simulation progress in tab Console. It contains containing details that allows to evaluate how well your numerical simulation is doing. Here’s a brief description of what each of those columns means.</p>
<table class='table borderless' id="solverTable">
<thead>
<tr>
<th>a</th>
<th>b</th>
<th>c</th>
<th>d</th>
<th>e</th>
<th>f</th>
<th>g</th>
<th>h</th>
</tr>
</thead>
<tr>
<td> 1 </td>
<td> 1.30E+001 </td>
<td> 1.16E-002 </td>
<td> (0) </td>
<td> 2.10E+002 </td>
<td> 2.62E+028 </td>
<td> < -10474 1 | 15 > </td>
<td> [199-190] </td>
</tr>
<tr>
<td> 1 </td>
<td> 2.50E+001 </td>
<td> 7.35E-003 </td>
<td> (-1) </td>
<td> 2.93E-001 </td>
<td> 5.15E+000 </td>
<td> < -3237 1 | 13> </td>
<td> [117-200] </td>
</tr>
<tr>
<td> 2 </td>
<td> 2.80E+001 </td>
<td> 5.13E-001 </td>
<td> (-16) </td>
<td> 1.75E-001 </td>
<td> 1.69E-001 </td>
<td> < -1357 1 | 5> </td>
<td> [63-1] </td>
</tr>
<tr>
<td> 2 </td>
<td> 3.00E+001 </td>
<td> 2.05E-002 </td>
<td> (-2) </td>
<td> 8.07E-002 </td>
<td> 2.67E-002 </td>
<td> < -3286 1 | 11> </td>
<td> [21-13] </td>
</tr>
<tr>
<td> 3 </td>
<td> 3.20E+001 </td>
<td> 1.20E-001 </td>
<td> (-10) </td>
<td> 8.75E-002 </td>
<td> 2.44E-002 </td>
<td> < -2342 1 | 7> </td>
<td> [36-1] </td>
</tr>
<tr>
<td> 3 </td>
<td> 3.40E+001 </td>
<td> 5.18E-003 </td>
<td> (-3) </td>
<td> 2.13E-002 </td>
<td> 3.59E-003 </td>
<td> < -3277 1 | 10> </td>
<td> [6-6] </td>
</tr>
<tr>
<td> 4 </td>
<td> 3.60E+001 </td>
<td> 2.14E-002 </td>
<td> (-2) </td>
<td> 5.57E-002 </td>
<td> 6.13E-003 </td>
<td> < -3146 1 | 9> </td>
<td> [24-2] </td>
</tr>
<tr>
<td> 4 </td>
<td> 3.80E+001 </td>
<td> 2.18E-003 </td>
<td> (-7) </td>
<td> 7.33E-003 </td>
<td> 3.15E-004 </td>
<td> < -3233 1 | 11> </td>
<td> [9-5] </td>
</tr>
<tr>
<td> 5 </td>
<td> 4.00E+001 </td>
<td> 1.52E-002 </td>
<td> (-1) </td>
<td> 4.22E-002 </td>
<td> 3.45E-004 </td>
<td> < -3141 1 | 10> </td>
<td> [27-3] </td>
</tr>