Chaste  Build::
AlarconHaematocritSolver.cpp
1 //
3 //Copyright (c) 2005-2015, University of Oxford.
4 //All rights reserved.
5 //
6 //University of Oxford means the Chancellor, Masters and Scholars of the
7 //University of Oxford, having an administrative office at Wellington
8 //Square, Oxford OX1 2JD, UK.
9 //
10 //This file is part of Chaste.
11 //
12 //Redistribution and use in source and binary forms, with or without
13 //modification, are permitted provided that the following conditions are met:
14 // * Redistributions of source code must retain the above copyright notice,
15 // this list of conditions and the following disclaimer.
16 // * Redistributions in binary form must reproduce the above copyright notice,
17 // this list of conditions and the following disclaimer in the documentation
18 // and/or other materials provided with the distribution.
19 // * Neither the name of the University of Oxford nor the names of its
20 // contributors may be used to endorse or promote products derived from this
21 // software without specific prior written permission.
22 //
23 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 //AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 //IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 //ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 //LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 //CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 //GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS unsignedERRUPTION)
30 //HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 //LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 //OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 #include "AlarconHaematocritSolver.hpp"
35 #include "LinearSystem.hpp"
36 #include "VesselNode.hpp"
37 #include "Vessel.hpp"
38 #include "Exception.hpp"
39 #include "ReplicatableVector.hpp"
40 
41 template<unsigned DIM>
43  mTHR(2.5),
44  mAlpha(0.5),
45  mHaematocrit(0.45)
46 {
47 
48 }
49 
50 template<unsigned DIM>
52 {
53 
54 }
55 
56 template<unsigned DIM>
58 {
59  // Give the vessels unique Ids
60  std::vector<boost::shared_ptr<Vessel<DIM> > > vessels = this->mpNetwork->GetVessels();
61  for(unsigned idx=0; idx<vessels.size(); idx++)
62  {
63  vessels[idx]->SetId(idx);
64  }
65 
66  // Set up the linear system
67  PetscInt lhsVectorSize = vessels.size();
68  unsigned max_vessels_per_branch = 5; // this could be increased if needed
69 
70  // If there are very few vessels, e.g. unit tests, set the number of non zeros to number of vessels
71  if(vessels.size() < max_vessels_per_branch)
72  {
73  max_vessels_per_branch = unsigned(lhsVectorSize);
74  }
75  LinearSystem linearSystem(lhsVectorSize, max_vessels_per_branch);
76 
77  // Currently `LinearSystem` defaults to an iterative solver for <6 dof. This check just highlights that this
78  // is happening.
79  if(lhsVectorSize > 6)
80  {
81  linearSystem.SetPcType("lu");
82 
83  // Use hypre if it is installed
84  #ifndef PETSC_HAVE_HYPRE
85  linearSystem.SetPcType("hypre");
86  #endif //PETSC_HAVE_HYPRE
87  linearSystem.SetKspType("preonly");
88  }
89 
90  for(unsigned idx=0; idx<vessels.size(); idx++)
91  {
92  // Always have a diagonal entry for system, this sets zero haematocrit by default
93  linearSystem.SetMatrixElement(idx, idx, 1);
94 
95  // Set arterial haematocrit for input nodes
96  if(vessels[idx]->GetStartNode()->GetFlowProperties()->IsInputNode() or vessels[idx]->GetEndNode()->GetFlowProperties()->IsInputNode())
97  {
98  linearSystem.SetRhsVectorElement(idx, mHaematocrit);
99  }
100  // Set rhs to zero for no flow vessels. It should already be zero, but this explicitly sets it for clarity
101  else if(vessels[idx]->GetFlowProperties()->GetFlowRate()==0.0*unit::metre_cubed_per_second)
102  {
103  linearSystem.SetRhsVectorElement(idx, 0.0);
104  }
105  else
106  {
107  // Identify the inflow node for this vessel
108  boost::shared_ptr<VesselNode<DIM> > p_inflow_node;
109  units::quantity<unit::flow_rate> flow_rate = vessels[idx]->GetFlowProperties()->GetFlowRate();
110  if(flow_rate >0.0 * unit::metre_cubed_per_second)
111  {
112  p_inflow_node = vessels[idx]->GetStartNode();
113  }
114  else
115  {
116  p_inflow_node = vessels[idx]->GetEndNode();
117  }
118 
119  // Identify the number of 'parent' and 'competitor' vessels. Parent vessels feed into the current vessel. Competitor vessels share
120  // a parent vessel and do not feed in, i.e. they compete for haematocrit
121  if(p_inflow_node->GetNumberOfSegments()>1)
122  {
123  std::vector<boost::shared_ptr<Vessel<DIM> > > parent_vessels;
124  std::vector<boost::shared_ptr<Vessel<DIM> > > competitor_vessels;
125  for(unsigned jdx=0; jdx<p_inflow_node->GetSegments().size(); jdx++)
126  {
127  // if not this vessel
128  if(p_inflow_node->GetSegment(jdx)->GetVessel()!=vessels[idx])
129  {
130  units::quantity<unit::flow_rate> inflow_rate = p_inflow_node->GetSegment(jdx)->GetVessel()->GetFlowProperties()->GetFlowRate();
131  if(p_inflow_node->GetSegment(jdx)->GetVessel()->GetEndNode()==p_inflow_node)
132  {
133  if(inflow_rate>0.0 * unit::metre_cubed_per_second)
134  {
135  parent_vessels.push_back(p_inflow_node->GetSegment(jdx)->GetVessel());
136  }
137  else if(inflow_rate<0.0 * unit::metre_cubed_per_second)
138  {
139  competitor_vessels.push_back(p_inflow_node->GetSegment(jdx)->GetVessel());
140  }
141  }
142  if(p_inflow_node->GetSegment(jdx)->GetVessel()->GetStartNode()==p_inflow_node)
143  {
144  if(inflow_rate>0.0 * unit::metre_cubed_per_second)
145  {
146  competitor_vessels.push_back(p_inflow_node->GetSegment(jdx)->GetVessel());
147  }
148  else if(inflow_rate<0.0 * unit::metre_cubed_per_second)
149  {
150  parent_vessels.push_back(p_inflow_node->GetSegment(jdx)->GetVessel());
151  }
152  }
153  }
154  }
155 
156  // If there are no competitor vessels the haematocrit is just the sum of the parent values
157  if(competitor_vessels.size()==0)
158  {
159  for(unsigned jdx=0; jdx<parent_vessels.size();jdx++)
160  {
161  linearSystem.SetMatrixElement(idx, parent_vessels[jdx]->GetId(), -1);
162  }
163  }
164  else
165  {
166  // This is for compatibility with the paper, we could allow for more if needed
167  if(competitor_vessels.size()>1 or parent_vessels.size()>1)
168  {
169  EXCEPTION("This solver can only work with branches with connectivity 3");
170  }
171 
172  // There is a bifurcation, apply a haematocrit splitting rule
173  units::quantity<unit::length> my_radius = vessels[idx]->GetRadius();
174  units::quantity<unit::length> competitor_radius = competitor_vessels[0]->GetRadius();
175  units::quantity<unit::velocity> my_velocity = units::fabs(flow_rate)/(M_PI * my_radius * my_radius);
176  units::quantity<unit::velocity> competitor_velocity = units::fabs(competitor_vessels[0]->GetFlowProperties()->GetFlowRate())/(M_PI * competitor_radius * competitor_radius);
177 
178  if(my_velocity>mTHR*competitor_velocity)
179  {
180  // I get all the haematocrit
181  linearSystem.SetMatrixElement(idx, parent_vessels[0]->GetId(), -1);
182  }
183  else if(my_velocity*mTHR>competitor_velocity)
184  {
185  // I get some haematocrit
186  if(my_velocity>competitor_velocity)
187  {
188  linearSystem.SetMatrixElement(idx, parent_vessels[0]->GetId(), -1.0/(1.0+(competitor_velocity/(my_velocity*mAlpha))));
189  }
190  else if(my_velocity<competitor_velocity)
191  {
192  linearSystem.SetMatrixElement(idx, parent_vessels[0]->GetId(), -1.0/(1.0+(mAlpha*competitor_velocity/(my_velocity))));
193  }
194  else
195  {
196  // Velocities are the same, but the partioining rule is not symmetric. This is not covered in the paper.
197  // Choose side for partitioning by vessel id so that there is some notion of conservation of H
198  if(vessels[idx]->GetId()>competitor_vessels[0]->GetId())
199  {
200  linearSystem.SetMatrixElement(idx, parent_vessels[0]->GetId(), -1.0/(1.0+(mAlpha*competitor_velocity/(my_velocity))));
201  }
202  else
203  {
204  linearSystem.SetMatrixElement(idx, parent_vessels[0]->GetId(), -1.0/(1.0+(competitor_velocity/(my_velocity*mAlpha))));
205  }
206  }
207  }
208  }
209  }
210  }
211  }
212 
213  // Do the solve
214  Vec solution = PetscTools::CreateVec(vessels.size());
215  linearSystem.AssembleFinalLinearSystem();
216  solution = linearSystem.Solve();
217  ReplicatableVector a(solution);
218 
219  // assign haematocrit levels to vessels
220  for (unsigned idx = 0; idx < vessels.size(); idx++)
221  {
222  for (unsigned jdx = 0; jdx < vessels[idx]->GetNumberOfSegments(); jdx++)
223  {
224  vessels[idx]->GetSegments()[jdx]->GetFlowProperties()->SetHaematocrit(a[idx]);
225  }
226  }
227 
228  PetscTools::Destroy(solution);
229 }
230 
231 template<unsigned DIM>
232 void AlarconHaematocritSolver<DIM>::SetTHR(units::quantity<unit::dimensionless> THR)
233 {
234  mTHR = THR;
235  assert(mTHR > 1);
236 }
237 
238 template<unsigned DIM>
239 void AlarconHaematocritSolver<DIM>::SetAlpha(units::quantity<unit::dimensionless> Alpha)
240 {
241  mAlpha = Alpha;
242  assert(mAlpha < 1);
243  assert(mAlpha > 0);
244 }
245 
246 template<unsigned DIM>
247 void AlarconHaematocritSolver<DIM>::SetHaematocrit(units::quantity<unit::dimensionless> haematocrit)
248 {
249  mHaematocrit = haematocrit;
250 }
251 
252 // Explicit instantiation
253 template class AlarconHaematocritSolver<2>;
254 template class AlarconHaematocritSolver<3>;
units::quantity< unit::dimensionless > mAlpha
boost::shared_ptr< VesselNetwork< DIM > > mpNetwork
void SetHaematocrit(units::quantity< unit::dimensionless > haematocrit)
void SetAlpha(units::quantity< unit::dimensionless > alpha)
units::quantity< unit::dimensionless > mHaematocrit
units::quantity< unit::dimensionless > mTHR
void SetTHR(units::quantity< unit::dimensionless > thr)