Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PulseWaveSystem.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File PulseWaveSystem.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Generic timestepping for Pulse Wave Solver
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
38 #include <MultiRegions/ContField1D.h>
39 #include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>
40 #include <PulseWaveSolver/EquationSystems/PulseWaveSystem.h>
41 #include <LibUtilities/BasicUtils/Timer.h>
42 
43 
44 namespace Nektar
45 {
46 
47  /**
48  * @class PulseWaveSystem
49  *
50  * Initialises the arterial subdomains in m_vessels and sets up
51  * all domain-linking conditions (bifurcations, junctions,
52  * merging flows). Detects the network structure and assigns
53  * boundary conditons. Also provides the underlying timestepping
54  * framework for pulse wave solvers including the general
55  * timestepping routines.
56  */
57 
58  /**
59  * Processes SolverInfo parameters from the session file and sets
60  * up timestepping-specific code.
61  *
62  * @param m_Session Session object to read parameters from.
63  */
64  PulseWaveSystem::PulseWaveSystem(const LibUtilities::SessionReaderSharedPtr& m_session)
65  : UnsteadySystem(m_session)
66  {
67  }
68 
69  /**
70  * Destructor
71  */
72  PulseWaveSystem::~PulseWaveSystem()
73  {
74  }
75 
76  /**
77  * Initialisation routine for multidomain solver. Sets up the
78  * expansions for every arterial segment (m_vessels) and for one
79  * complete field m_outfield which is needed to write the
80  * postprocessing output. Also determines which upwind strategy
81  * is used (currently only upwinding scheme available) and reads
82  * blodd flow specific parameters from the inputfile
83  *
84  */
85  void PulseWaveSystem::v_InitObject()
86  {
87  // Initialise base class
88  UnsteadySystem::v_InitObject();
89 
90  // Read the geometry and the expansion information
91  m_nDomains = m_graph->GetDomain().size();
92 
93  // Determine projectiontype
94  ASSERTL0(m_session->MatchSolverInfo("Projection","DisContinuous"),"Pulse solver only set up for Discontinuous projections");
95  m_projectionType = MultiRegions::eDiscontinuous;
96  ASSERTL0(m_graph->GetMeshDimension() == 1,"Pulse wave solver only set up for expansion dimension equal to 1");
97 
98  int i;
99  m_nVariables = m_session->GetVariables().size();
100 
101  m_fields = Array<OneD, MultiRegions::ExpListSharedPtr> (m_nVariables);
102  m_vessels = Array<OneD, MultiRegions::ExpListSharedPtr> (m_nVariables*m_nDomains);
103 
104  const std::vector<SpatialDomains::CompositeMap> domain = m_graph->GetDomain();
105 
106  SpatialDomains::BoundaryConditions Allbcs(m_session, m_graph);
107 
108  // Set up domains and put geometry to be only one space dimension.
109  int cnt = 0;
110  bool SetToOneSpaceDimension = true;
111 
112  if(m_session->DefinesCmdLineArgument("SetToOneSpaceDimension"))
113  {
114  std::string cmdline = m_session->GetCmdLineArgument<std::string>("SetToOneSpaceDimension");
115  if(boost::to_upper_copy(cmdline) == "FALSE")
116  {
117  SetToOneSpaceDimension = false;
118  }
119  }
120 
121  for(i = 0 ; i < m_nDomains; ++i)
122  {
123  for(int j = 0; j < m_nVariables; ++j)
124  {
125  m_vessels[cnt++] = MemoryManager<MultiRegions::DisContField1D>
126  ::AllocateSharedPtr(m_session, m_graph, domain[i],
127  Allbcs,
128  m_session->GetVariable(j),
129  SetToOneSpaceDimension);
130  }
131  }
132 
133  // Reset coeff and phys space to be continuous over all domains
134  int totcoeffs = 0;
135  int totphys = 0;
136  m_fieldPhysOffset = Array<OneD, int>(m_nDomains+1,0);
137  for(i = 0; i < m_nDomains; ++i)
138  {
139  totcoeffs += m_vessels[i*m_nVariables]->GetNcoeffs();
140 
141  m_fieldPhysOffset[i] = totphys;
142  totphys += m_vessels[i*m_nVariables]->GetTotPoints();
143  }
144  m_fieldPhysOffset[m_nDomains] = totphys;
145 
146  for(int n = 0; n < m_nVariables; ++n)
147  {
148  Array<OneD, NekDouble> coeffs(totcoeffs,0.0);
149  Array<OneD, NekDouble> phys(totphys,0.0);
150  Array<OneD, NekDouble> tmpcoeffs,tmpphys;
151 
152  m_vessels[n]->SetCoeffsArray(coeffs);
153  m_vessels[n]->SetPhysArray(phys);
154 
155  int cnt = m_vessels[n]->GetNcoeffs();
156  int cnt1 = m_vessels[n]->GetTotPoints();
157 
158  for(i = 1; i < m_nDomains; ++i)
159  {
160  m_vessels[i*m_nVariables+n]->SetCoeffsArray(tmpcoeffs = coeffs + cnt);
161  m_vessels[i*m_nVariables+n]->SetPhysArray (tmpphys = phys + cnt1);
162  cnt += m_vessels[i*m_nVariables+n]->GetNcoeffs();
163  cnt1 += m_vessels[i*m_nVariables+n]->GetTotPoints();
164  }
165  }
166 
167  // Set Default Parameter
168  m_session->LoadParameter("Time", m_time, 0.0);
169  m_session->LoadParameter("TimeStep", m_timestep, 0.01);
170  m_session->LoadParameter("NumSteps", m_steps, 0);
171  m_session->LoadParameter("IO_CheckSteps", m_checksteps, m_steps);
172  m_session->LoadParameter("FinTime", m_fintime, 0);
173  m_session->LoadParameter("NumQuadPointsError", m_NumQuadPointsError, 0);
174 
175  m_fields[0] = m_vessels[0];
176  m_fields[1] = m_vessels[1];
177 
178  // Zero all physical fields initially.
179  ZeroPhysFields();
180 
181  // If Discontinuous Galerkin determine upwinding method to use
182  for (int i = 0; i < (int)SIZE_UpwindTypePulse; ++i)
183  {
184  bool match;
185  m_session->MatchSolverInfo("UPWINDTYPEPULSE", UpwindTypeMapPulse[i], match, false);
186  if (match)
187  {
188  m_upwindTypePulse = (UpwindTypePulse) i;
189  break;
190  }
191  }
192  // Load solver- specific parameters
193  m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
194  // Load blood density
195  m_session->LoadParameter("rho", m_rho, 0.5);
196  // Load external pressure
197  m_session->LoadParameter("pext", m_pext, 0.0);
198 
199  int nq = 0;
200  /**
201  * Gets the Material Properties of each arterial segment
202  * specified in the inputfile from section MaterialProperties
203  * Also gets the Area at static equilibrium A_0 specified in the
204  * inputfile.
205  *
206  * Having found these points also extract the values at the
207  * trace points and the normal direction consistent with the
208  * left adjacent definition of Fwd and Bwd
209  */
210  m_beta = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
211  m_beta_trace = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
212  m_A_0 = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
213  m_A_0_trace = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
214  m_trace_fwd_normal = Array<OneD, Array<OneD, NekDouble> >(m_nDomains);
215 
216  for (int omega = 0; omega < m_nDomains; omega++)
217  {
218  nq = m_vessels[2*omega]->GetNpoints();
219  m_fields[0] = m_vessels[2*omega];
220 
221  m_beta[omega] = Array<OneD, NekDouble>(nq);
222  EvaluateFunction("beta", m_beta[omega],"MaterialProperties",0.0,omega);
223 
224  m_A_0[omega] = Array<OneD, NekDouble>(nq);
225  EvaluateFunction("A_0", m_A_0[omega],"A_0",0.0,omega);
226 
227  int nqTrace = GetTraceTotPoints();
228 
229  m_beta_trace[omega] = Array<OneD, NekDouble>(nqTrace);
230  m_fields[0]->ExtractTracePhys(m_beta[omega],m_beta_trace[omega]);
231 
232  m_A_0_trace[omega] = Array<OneD, NekDouble>(nqTrace);
233  m_fields[0]->ExtractTracePhys(m_A_0[omega],m_A_0_trace[omega]);
234 
235 
236  if(SetToOneSpaceDimension)
237  {
238  m_trace_fwd_normal[omega] = Array<OneD, NekDouble>(nqTrace,0.0);
239 
240  MultiRegions::ExpListSharedPtr trace = m_fields[0]->GetTrace();
241  int nelmt_trace = trace->GetExpSize();
242 
243  Array<OneD, Array<OneD, NekDouble> > normals(nelmt_trace);
244 
245  for(int i = 0 ; i < nelmt_trace; ++i)
246  {
247  normals[i] = m_trace_fwd_normal[omega]+i;
248  }
249 
250  // need to set to 1 for consistency since boundary
251  // conditions may not have coordim=1
252  trace->GetExp(0)->GetGeom()->SetCoordim(1);
253 
254  trace->GetNormals(normals);
255  }
256  }
257 
258  SetUpDomainInterfaces();
259 
260  }
261 
262  void PulseWaveSystem::SetUpDomainInterfaces(void)
263  {
264  map<int,std::vector<InterfacePointShPtr> > VidToDomain;
265  map<int,std::vector<InterfacePointShPtr> >::iterator iter;
266 
267  // loop over domain and find out if we have any undefined
268  // boundary conditions representing interfaces. If so make a
269  // map based around vid and storing the domains that are
270  // part of interfaces.
271  for(int omega = 0; omega < m_nDomains; ++omega)
272  {
273  int vesselID = omega*m_nVariables;
274 
275  for(int i = 0; i < 2; ++i)
276  {
277  if(m_vessels[vesselID]->GetBndConditions()[i]->GetBoundaryConditionType() == SpatialDomains::eNotDefined)
278  {
279  // Get Vid of interface
280  int vid = m_vessels[vesselID]->UpdateBndCondExpansion(i)->GetExp(0)->GetGeom()->GetVid(0);
281  cout<<"Shared vertex id: "<<vid<<endl;
282  MultiRegions::ExpListSharedPtr trace = m_vessels[vesselID]->GetTrace();
283  InterfacePointShPtr Ipt;
284 
285  bool finish = false;
286  // find which elmt, the lcoal vertex and the data offset of point
287  for(int n = 0; n < m_vessels[vesselID]->GetExpSize(); ++n)
288  {
289  for(int p = 0; p < 2; ++p)
290  {
291  if(m_vessels[vesselID]->GetTraceMap()->GetElmtToTrace()[n][p]->as<LocalRegions::Expansion>()->GetGeom()->GetVid(0) == vid)
292  {
293  int eid = m_vessels[vesselID]->GetTraceMap()->GetElmtToTrace()[n][p]->GetElmtId();
294 
295  int tid = m_vessels[vesselID]->GetTrace()->GetCoeff_Offset(eid);
296  Ipt = MemoryManager<InterfacePoint>::AllocateSharedPtr(vid,omega,n,p,tid,i);
297 
298  cout<<"Global Vid of interface point: "<<vid<<endl;
299  cout<<"Domain interface point belongs to: "<<omega<<endl;
300  cout<<"Element id of vertex: "<<n<<endl;
301  cout<<"Vertex id within local element: "<<p<<endl;
302  cout<<"Element id within the trace: "<<tid<<endl;
303  cout<<"Position of boundary condition in region: "<<i<<endl;
304 
305  finish = true;
306  break;
307  }
308  }
309  if(finish == true)
310  {
311  break;
312  }
313  }
314 
315  VidToDomain[vid].push_back(Ipt);
316 
317  // finally reset boundary condition to Dirichlet
318  m_vessels[vesselID]->GetBndConditions()[i]
319  ->SetBoundaryConditionType(SpatialDomains::eDirichlet);
320 
321  m_vessels[vesselID+1]->GetBndConditions()[i]
322  ->SetBoundaryConditionType(SpatialDomains::eDirichlet);
323  }
324  }
325  }
326 
327  // loop over map and set up Interface information;
328  for(iter = VidToDomain.begin(); iter != VidToDomain.end(); ++iter)
329  {
330  if(iter->second.size() == 2) // Vessel jump interface
331  {
332  m_vesselJcts.push_back(iter->second);
333  }
334  else if(iter->second.size() == 3) // Bifurcation or Merging junction.
335  {
336  int nbeg = 0;
337  int nend = 0;
338 
339  // Determine if bifurcation or merging junction
340  // through number of elemnt vertices that meet at
341  // junction. Only one vertex using a m_elmtVert=1
342  // indicates a bifurcation
343  for(int i = 0; i < 3; ++i)
344  {
345  if(iter->second[i]->m_elmtVert == 0)
346  {
347  nbeg += 1;
348  }
349  else
350  {
351  nend += 1;
352  }
353  }
354 
355  // Set up Bifurcation information
356  if(nbeg == 2)
357  {
358  // ensure first InterfacePoint is parent
359  if(iter->second[0]->m_elmtVert == 1) //m_elmtVert: Vertex id in local element
360  {
361  m_bifurcations.push_back(iter->second);
362  }
363  else
364  {
365  //order points according to Riemann solver convention
366  InterfacePointShPtr I;
367  //find merging vessel
368  if(iter->second[1]->m_elmtVert == 1)
369  {
370  I = iter->second[0];
371  iter->second[0] = iter->second[1];
372  iter->second[1] = I;
373  }
374  else if (iter->second[2]->m_elmtVert == 1)
375  {
376  I = iter->second[0];
377  iter->second[0] = iter->second[2];
378  iter->second[2] = I;
379  }
380  NEKERROR(ErrorUtil::ewarning,"This routine has not been checked");
381  }
382  }
383  else
384  {
385  // ensure last InterfacePoint is merged vessel
386  if(iter->second[0]->m_elmtVert == 0)
387  {
388  m_mergingJcts.push_back(iter->second);
389  }
390  else
391  {
392  //order points according to Riemann solver convention
393  InterfacePointShPtr I;
394  //find merging vessel
395  if(iter->second[1]->m_elmtVert == 0)
396  {
397  I = iter->second[0];
398  iter->second[0] = iter->second[1];
399  iter->second[1] = I;
400  }
401  else if (iter->second[2]->m_elmtVert == 0)
402  {
403  I = iter->second[0];
404  iter->second[0] = iter->second[2];
405  iter->second[2] = I;
406  }
407  NEKERROR(ErrorUtil::ewarning,"This routine has not been checked");
408  }
409  }
410 
411  }
412  else
413  {
414  ASSERTL0(false,"Unknown junction type");
415  }
416  }
417  }
418 
419 
420  /**
421  * Initialisation routine for multiple subdomain case. Sets the
422  * initial conditions for all arterial subdomains read from the
423  * inputfile. Sets the material properties and the A_0 area for
424  * all subdomains and fills the domain-linking boundary
425  * conditions with the initial values of their domain.
426  */
427  void PulseWaveSystem::v_DoInitialise()
428  {
429 
430  if (m_session->GetComm()->GetRank() == 0)
431  {
432  cout << "Initial Conditions:" << endl;
433  }
434 
435  /* Loop over all subdomains to initialize all with the Initial
436  * Conditions read from the inputfile*/
437  for (int omega = 0; omega < m_nDomains; omega++)
438  {
439  m_fields[0] = m_vessels[m_nVariables*omega];
440  m_fields[1] = m_vessels[m_nVariables*omega+1];
441 
442  if (m_session->GetComm()->GetRank() == 0)
443  {
444  cout << "Subdomain = " <<omega<<endl;
445  }
446 
447  SetInitialConditions(0.0,0,omega);
448  }
449  // Reset to first definition
450  m_fields[0] = m_vessels[0];
451  m_fields[1] = m_vessels[1];
452 
453  }
454 
455  /**
456  * NEEDS Updating:
457  *
458  * DoSolve routine for PulseWavePropagation with multiple
459  * subdomains taken from UnsteadySystem and modified for
460  * multidomain case. Initialises the time integration scheme (as
461  * specified in the session file), and perform the time
462  * integration. Within the timestepping loop the following is
463  * done: 1. Link all arterial segments according to the network
464  * structure, solve the Riemann problem between different
465  * arterial segments and assign the values to the boundary
466  * conditions (LinkSubdomains) 2. Every arterial segment is
467  * solved independentl for this timestep. This is done by handing
468  * the solution vector \f$ \mathbf{u} \f$ and the right hand side
469  * m_ode, which is the PulseWavePropagation class in this example
470  * over to the time integration scheme
471  */
472  void PulseWaveSystem::v_DoSolve()
473  {
474  NekDouble IntegrationTime = 0.0;
475  int i,n,nchk = 1;
476 
477  Array<OneD, Array<OneD,NekDouble> > fields(m_nVariables);
478 
479  for(int i = 0; i < m_nVariables; ++i)
480  {
481  fields[i] = m_vessels[i]->UpdatePhys();
482  m_fields[i]->SetPhysState(false);
483  }
484 
485  m_intSoln = m_intScheme->InitializeScheme(
486  m_timestep,fields,m_time,m_ode);
487 
488  // Time loop
489  for(n = 0; n < m_steps; ++n)
490  {
491  Timer timer;
492  timer.Start();
493  fields = m_intScheme->TimeIntegrate(n,m_timestep,m_intSoln,m_ode);
494  //cout<<"integration: "<<fields[0][fields[0].num_elements()-1]<<endl;
495  m_time += m_timestep;
496  timer.Stop();
497  IntegrationTime += timer.TimePerTest(1);
498 
499  // Write out status information.
500  if(m_session->GetComm()->GetRank() == 0 && !((n+1)%m_infosteps))
501  {
502  cout << "Steps: " << n+1
503  << "\t Time: " << m_time
504  << "\t Time-step: " << m_timestep << "\t" << endl;
505  }
506 
507  // Transform data if needed
508  if(!((n+1)%m_checksteps))
509  {
510  for (i = 0; i < m_nVariables; ++i)
511  {
512  int cnt = 0;
513  for (int omega = 0; omega < m_nDomains; omega++)
514  {
515  m_vessels[omega*m_nVariables+i]->FwdTrans(fields[i]+cnt,
516  m_vessels[omega*m_nVariables+i]->UpdateCoeffs());
517  cnt += m_vessels[omega*m_nVariables+i]->GetTotPoints();
518  }
519  }
520  CheckPoint_Output(nchk++);
521  }
522 
523  }//end of timeintegration
524 
525  //Copy Array To Vessel Phys Fields
526  for(int i = 0; i < m_nVariables; ++i)
527  {
528  Vmath::Vcopy(fields[i].num_elements(), fields[i],1,m_vessels[i]->UpdatePhys(),1);
529  }
530 
531  cout <<"Time-integration timing : "
532  << IntegrationTime << " s" << endl << endl;
533  }
534 
535 
536 
537  void PulseWaveSystem::FillDataFromInterfacePoint(InterfacePointShPtr &I,
538  const Array<OneD, const Array<OneD, NekDouble> >&fields,
539  NekDouble &A, NekDouble &u,
540  NekDouble &beta, NekDouble &A_0)
541  {
542  int omega, traceId, eid, vert, phys_offset, vesselID;
543 
544 
545  // Parent vessel
546  omega = I->m_domain;
547  traceId = I->m_traceId;
548  eid = I->m_elmt;
549  vert = I->m_elmtVert;
550  vesselID = omega*m_nVariables;
551 
552  phys_offset = m_vessels[vesselID]->GetPhys_Offset(eid);
553 
554  m_vessels[vesselID]->GetExp(eid)->
555  GetVertexPhysVals(vert, fields[0]+m_fieldPhysOffset[omega]+phys_offset, A);
556  m_vessels[vesselID]->GetExp(eid)->
557  GetVertexPhysVals(vert, fields[1]+m_fieldPhysOffset[omega]+phys_offset, u);
558 
559  beta = m_beta_trace[omega][traceId];
560  A_0 = m_A_0_trace [omega][traceId];
561  }
562 
563  void PulseWaveSystem::EnforceInterfaceConditions(const Array<OneD, const Array<OneD, NekDouble> > &fields)
564  {
565  int dom, bcpos;
566  Array<OneD, NekDouble> Au(3),uu(3),beta(3),A_0(3);
567 
568 
569  // Enfore Bifurcations;
570  for(int n = 0; n < m_bifurcations.size(); ++n)
571  {
572  FillDataFromInterfacePoint(m_bifurcations[n][0],fields,
573  Au[0],uu[0],beta[0],A_0[0]);
574  FillDataFromInterfacePoint(m_bifurcations[n][1],fields,
575  Au[1],uu[1],beta[1],A_0[1]);
576  FillDataFromInterfacePoint(m_bifurcations[n][2],fields,
577  Au[2],uu[2],beta[2],A_0[2]);
578 
579  // Solve the Riemann problem for a bifurcation
580  BifurcationRiemann(Au, uu, beta, A_0);
581 
582  // Store the values into the right positions:
583  for(int i = 0; i < 3; ++i)
584  {
585  dom = m_bifurcations[n][i]->m_domain;
586  bcpos = m_bifurcations[n][i]->m_bcPosition;
587  m_vessels[dom*m_nVariables] ->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = Au[i];
588  m_vessels[dom*m_nVariables+1]->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = uu[i];
589  }
590  }
591 
592  // Enfore Bifurcations;
593  for(int n = 0; n < m_mergingJcts.size(); ++n)
594  {
595  // Merged vessel
596  FillDataFromInterfacePoint(m_mergingJcts[n][0],fields,
597  Au[0],uu[0],beta[0],A_0[0]);
598  FillDataFromInterfacePoint(m_mergingJcts[n][1],fields,
599  Au[1],uu[1],beta[1],A_0[1]);
600  FillDataFromInterfacePoint(m_mergingJcts[n][2],fields,
601  Au[2],uu[2],beta[2],A_0[2]);
602 
603  // Solve the Riemann problem for a merging vessel
604  MergingRiemann(Au, uu, beta, A_0);
605 
606  // Store the values into the right positions:
607  for(int i = 0; i < 3; ++i)
608  {
609  int dom = m_mergingJcts[n][i]->m_domain;
610  int bcpos = m_mergingJcts[n][i]->m_bcPosition;
611  m_vessels[dom*m_nVariables] ->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = Au[i];
612  m_vessels[dom*m_nVariables+1]->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = uu[i];
613  }
614  }
615 
616  for(int n = 0; n < m_vesselJcts.size(); ++n)
617  {
618 
619  FillDataFromInterfacePoint(m_vesselJcts[n][0],fields,
620  Au[0],uu[0],beta[0],A_0[0]);
621  FillDataFromInterfacePoint(m_vesselJcts[n][1],fields,
622  Au[1],uu[1],beta[1],A_0[1]);
623 
624  JunctionRiemann(Au, uu, beta, A_0);
625 
626  // Store the values into the right positions:
627  for(int i = 0; i < 2; ++i)
628  {
629  int dom = m_vesselJcts[n][i]->m_domain;
630  int bcpos = m_vesselJcts[n][i]->m_bcPosition;
631  m_vessels[dom*m_nVariables] ->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = Au[i];
632  m_vessels[dom*m_nVariables+1]->UpdateBndCondExpansion(bcpos)->UpdatePhys()[0] = uu[i];
633  //cout<<"Au: "<<Au[i]<<endl;
634  //cout<<"uu: "<<uu[i]<<endl;
635  //cout<<endl;
636  }
637 
638  }
639  }
640 
641 
642  /**
643  * Solves the Riemann problem at a bifurcation by assuming
644  * subsonic flow at both sides of the boundary and by
645  * applying conservation of mass and continuity of the total
646  * pressure \f$ \frac{p}{rho} + \frac{u^{2}}{2}. \f$ The
647  * other 3 missing equations come from the characteristic
648  * variables. For further information see "Pulse
649  * WavePropagation in the human vascular system" Section
650  * 3.4.4
651  */
652  void PulseWaveSystem::BifurcationRiemann(Array<OneD, NekDouble> &Au, Array<OneD, NekDouble> &uu,
653  Array<OneD, NekDouble> &beta, Array<OneD, NekDouble> &A_0)
654  {
655  NekDouble rho = m_rho;
656  Array<OneD, NekDouble> W(3);
657  Array<OneD, NekDouble> W_Au(3);
658  Array<OneD, NekDouble> P_Au(3);
659  Array<OneD, NekDouble> f(6);
660  Array<OneD, NekDouble> g(6);
661  Array<OneD, NekDouble> tmp(6);
662  Array<OneD, Array<OneD, NekDouble> > inv_J(6);
663  for (int i=0; i<6; i++)
664  {
665  inv_J[i] = Array<OneD, NekDouble> (6);
666  }
667  NekDouble k = 0.0;
668  NekDouble k1 = 0.0;
669  NekDouble k2 = 0.0;
670  NekDouble k3 = 0.0;
671 
672  int proceed = 1;
673  int iter = 0;
674  int MAX_ITER = 7;
675 
676  // Calculated from input
677  W[0] = uu[0] + 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
678  W[1] = uu[1] - 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
679  W[2] = uu[2] - 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
680 
681  // Tolerances for the algorithm
682  NekDouble Tol = 1.0e-10;
683 
684  // Newton Iteration
685  while ((proceed) && (iter < MAX_ITER))
686  {
687  iter = iter+1;
688 
689  // Calculate the constraint vector, six equations:
690  // 3 characteristic variables, mass conservation,
691  // total pressure
692  W_Au[0] = 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
693  W_Au[1] = 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
694  W_Au[2] = 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
695 
696  P_Au[0] = beta[0]*(sqrt(Au[0]) - sqrt(A_0[0]));
697  P_Au[1] = beta[1]*(sqrt(Au[1]) - sqrt(A_0[1]));
698  P_Au[2] = beta[2]*(sqrt(Au[2]) - sqrt(A_0[2]));
699 
700  f[0] = uu[0] + W_Au[0] - W[0];
701  f[1] = uu[1] - W_Au[1] - W[1];
702  f[2] = uu[2] - W_Au[2] - W[2];
703  f[3] = Au[0]*uu[0] - Au[1]*uu[1] - Au[2]*uu[2];
704  f[4] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[1]*uu[1] - 2.0/rho*P_Au[1];
705  f[5] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[2]*uu[2] - 2.0/rho*P_Au[2];
706 
707  // Calculate the wave speed at each vessel
708  NekDouble c1 = sqrt(beta[0]/(2*rho))*sqrt(sqrt(Au[0]));
709  NekDouble c2 = sqrt(beta[1]/(2*rho))*sqrt(sqrt(Au[1]));
710  NekDouble c3 = sqrt(beta[2]/(2*rho))*sqrt(sqrt(Au[2]));
711 
712  // Inverse Jacobian matrix J(x[n])^(-1), is
713  // already inverted here analytically
714  k = c1*Au[1]*c3+Au[0]*c3*c2+Au[2]*c1*c2;
715  k1 = (c1-uu[0])*k;
716  inv_J[0][0] = (-c2*uu[0]*c3*Au[0]+Au[2]*c2*c1*c1+Au[1]*c1*c1*c3)/k1;
717  inv_J[0][1] = Au[1]*(c2-uu[1])*c1*c3/k1;
718  inv_J[0][2] = Au[2]*(c3-uu[2])*c1*c2/k1;
719  inv_J[0][3] = c1*c2*c3/k1;
720  inv_J[0][4] = -0.5*c1*Au[1]*c3/k1;
721  inv_J[0][5] = -0.5*Au[2]*c1*c2/k1;
722 
723  k2 = (c2+uu[1])*k;
724  inv_J[1][0] = Au[0]*(c1+uu[0])*c2*c3/k2;
725  inv_J[1][1] = (c1*uu[1]*c3*Au[1]+Au[2]*c1*c2*c2+c3*c2*c2*Au[0])/k2;
726  inv_J[1][2] = -Au[2]*(c3-uu[2])*c1*c2/k2;
727  inv_J[1][3] = -c1*c2*c3/k2;
728  inv_J[1][4] = -0.5*(c1*Au[2]+Au[0]*c3)*c2/k2;
729  inv_J[1][5] = 0.5*Au[2]*c1*c2/k2;
730 
731  k3 = (c3+uu[2])*k;
732  inv_J[2][0] = Au[0]*(c1+uu[0])*c2*c3/k3;
733  inv_J[2][1] = -Au[1]*(c2-uu[1])*c1*c3/k3;
734  inv_J[2][2] = (c1*c2*uu[2]*Au[2]+c1*Au[1]*c3*c3+c2*c3*c3*Au[0])/k3;
735  inv_J[2][3] = -c1*c2*c3/k3;
736  inv_J[2][4] = 0.5*c1*Au[1]*c3/k3;
737  inv_J[2][5] = -0.5*(Au[1]*c1+c2*Au[0])*c3/k3;
738 
739  inv_J[3][0] = Au[0]*(Au[0]*c3*c2-uu[0]*c3*Au[1]-uu[0]*c2*Au[2])/k1;
740  inv_J[3][1] = -Au[0]*Au[1]*(c2-uu[1])*c3/k1;
741  inv_J[3][2] = -Au[0]*Au[2]*(c3-uu[2])*c2/k1;
742  inv_J[3][3] = -Au[0]*c3*c2/k1;
743  inv_J[3][4] = 0.5*Au[0]*Au[1]*c3/k1;
744  inv_J[3][5] = 0.5*Au[0]*c2*Au[2]/k1;
745 
746  inv_J[4][0] = Au[0]*Au[1]*(c1+uu[0])*c3/k2;
747  inv_J[4][1] = -Au[1]*(c1*Au[1]*c3+c1*uu[1]*Au[2]+c3*uu[1]*Au[0])/k2;
748  inv_J[4][2] = -Au[2]*Au[1]*(c3-uu[2])*c1/k2;
749  inv_J[4][3] = -c1*Au[1]*c3/k2;
750  inv_J[4][4] = -0.5*Au[1]*(c1*Au[2]+Au[0]*c3)/k2;
751  inv_J[4][5] = 0.5*Au[2]*Au[1]*c1/k2;
752 
753  inv_J[5][0] = Au[0]*Au[2]*(c1+uu[0])*c2/k3;
754  inv_J[5][1] = -Au[2]*Au[1]*(c2-uu[1])*c1/k3;
755  inv_J[5][2] = -Au[2]*(Au[2]*c1*c2+c1*uu[2]*Au[1]+c2*uu[2]*Au[0])/k3;
756  inv_J[5][3] = -Au[2]*c1*c2/k3;
757  inv_J[5][4] = 0.5*Au[2]*Au[1]*c1/k3;
758  inv_J[5][5] = -0.5*Au[2]*(Au[1]*c1+c2*Au[0])/k3;
759 
760 
761  // Solve the system by multiplying the Jacobian with the vector f:
762  // g = (inv_J)*f
763  for (int j=0; j<6; j++)
764  {
765  tmp[j] =0.0;
766  g[j] = 0.0;
767  }
768 
769  for (int j=0; j<6; j++)
770  {
771  for (int i=0; i<6; i++)
772  {
773  tmp[j] = inv_J[j][i]*f[i];
774  g[j] += tmp[j];
775  }
776  }
777 
778  // Update the solution: x_new = x_old - dx
779  uu[0] = uu[0] - g[0];
780  uu[1] = uu[1] - g[1];
781  uu[2] = uu[2] - g[2];
782  Au[0] = Au[0] - g[3];
783  Au[1] = Au[1] - g[4];
784  Au[2] = Au[2] - g[5];
785 
786  // Check if the error of the solution is smaller than Tol
787  if ((g[0]*g[0] + g[1]*g[1] + g[2]*g[2] + g[3]*g[3]+ g[4]*g[4] + g[5]*g[5]) < Tol)
788  {
789  proceed = 0;
790  }
791 
792  // Check if solver converges
793  if (iter >= MAX_ITER)
794  {
795  ASSERTL0(false,"Riemann solver for Bifurcation did not converge");
796  }
797  }
798  }
799 
800 
801 
802  /**
803  * Solves the Riemann problem at an merging flow condition by
804  * assuming subsonic flow at both sides of the boundary and by
805  * applying conservation of mass and continuity of the total
806  * pressure \f$ \frac{p}{rho} + \frac{u^{2}}{2}. \f$ The other 3
807  * missing equations come from the characteristic variables. For
808  * further information see "Pulse WavePropagation in the human
809  * vascular system" Section 3.4.4
810  */
811  void PulseWaveSystem::MergingRiemann(Array<OneD, NekDouble> &Au, Array<OneD, NekDouble> &uu, Array<OneD, NekDouble> &beta, Array<OneD, NekDouble> &A_0)
812  {
813  NekDouble rho = m_rho;
814  Array<OneD, NekDouble> W(3);
815  Array<OneD, NekDouble> W_Au(3);
816  Array<OneD, NekDouble> P_Au(3);
817  Array<OneD, NekDouble> f(6);
818  Array<OneD, NekDouble> g(6);
819  Array<OneD, NekDouble> tmp(6);
820  Array<OneD, Array<OneD, NekDouble> > inv_J(6);
821 
822  for (int i=0; i<6; i++)
823  {
824  inv_J[i] = Array<OneD, NekDouble> (6);
825  }
826 
827  NekDouble k = 0.0;
828  NekDouble k1 = 0.0;
829  NekDouble k2 = 0.0;
830  NekDouble k3 = 0.0;
831 
832  int proceed = 1;
833  int iter = 0;
834  int MAX_ITER = 7;
835 
836  // Calculated from input
837  W[0] = uu[0] - 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
838  W[1] = uu[1] + 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
839  W[2] = uu[2] + 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
840 
841  // Tolerances for the algorithm
842  NekDouble Tol = 1.0e-10;
843 
844  // Newton Iteration
845  while ((proceed) && (iter < MAX_ITER))
846  {
847  iter = iter+1;
848 
849  // Calculate the constraint vector, six equations:
850  // 3 characteristic variables, mass conservation,
851  // total pressure
852  W_Au[0] = 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
853  W_Au[1] = 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
854  W_Au[2] = 4*sqrt(beta[2]/(2*rho))*(sqrt(sqrt(Au[2])) - sqrt(sqrt(A_0[2])));
855 
856  P_Au[0] = beta[0]*(sqrt(Au[0]) - sqrt(A_0[0]));
857  P_Au[1] = beta[1]*(sqrt(Au[1]) - sqrt(A_0[1]));
858  P_Au[2] = beta[2]*(sqrt(Au[2]) - sqrt(A_0[2]));
859 
860  f[0] = uu[0] - W_Au[0] - W[0];
861  f[1] = uu[1] + W_Au[1] - W[1];
862  f[2] = uu[2] + W_Au[2] - W[2];
863  f[3] = Au[0]*uu[0] - Au[1]*uu[1] - Au[2]*uu[2];
864  f[4] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[1]*uu[1] - 2.0/rho*P_Au[1];
865  f[5] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[2]*uu[2] - 2.0/rho*P_Au[2];
866 
867  // Calculate the wave speed at each vessel
868  NekDouble c1 = sqrt(beta[0]/(2*rho))*sqrt(sqrt(Au[0]));
869  NekDouble c2 = sqrt(beta[1]/(2*rho))*sqrt(sqrt(Au[1]));
870  NekDouble c3 = sqrt(beta[2]/(2*rho))*sqrt(sqrt(Au[2]));
871 
872  // Inverse Jacobian matrix J(x[n])^(-1), is already inverted here analytically
873  k = c1*Au[1]*c3+Au[0]*c3*c2+Au[2]*c1*c2;
874  k1 = (c1+uu[0])*k;
875  inv_J[0][0] = (c2*uu[0]*c3*Au[0]+Au[2]*c2*c1*c1+Au[1]*c1*c1*c3)/k1;
876  inv_J[0][1] = Au[1]*(c2+uu[1])*c1*c3/k1;
877  inv_J[0][2] = Au[2]*(c3+uu[2])*c1*c2/k1;
878  inv_J[0][3] = c1*c2*c3/k1;
879  inv_J[0][4] = 0.5*Au[1]*c1*c3/k1;
880  inv_J[0][5] = 0.5*Au[2]*c1*c2/k1;
881 
882  k2 = (c2-uu[1])*k;
883  inv_J[1][0] = Au[0]*(c1-uu[0])*c2*c3/k2;
884  inv_J[1][1] = (-c1*uu[1]*c3*Au[1]+Au[2]*c1*c2*c2+c3*c2*c2*Au[0])/k2;
885  inv_J[1][2] = -Au[2]*(c3+uu[2])*c1*c2/k2;
886  inv_J[1][3] = -c1*c2*c3/k2;
887  inv_J[1][4] = 0.5*(c1*Au[2]+Au[0]*c3)*c2/k2;
888  inv_J[1][5] = -0.5*Au[2]*c1*c2/k2;
889 
890  k3 = (c3-uu[2])*k;
891  inv_J[2][0] = Au[0]*(c1-uu[0])*c2*c3/k3;
892  inv_J[2][1] = -Au[1]*(c2+uu[1])*c1*c3/k3;
893  inv_J[2][2] = -(c1*uu[2]*c2*Au[2]-Au[1]*c1*c3*c3-c2*c3*c3*Au[0])/k3;
894  inv_J[2][3] = -c1*c2*c3/k3;
895  inv_J[2][4] = -0.5*Au[1]*c1*c3/k3;
896  inv_J[2][5] = 0.5*(Au[1]*c1+Au[0]*c2)*c3/k3;
897 
898  inv_J[3][0] = -Au[0]*(Au[0]*c3*c2+uu[0]*c3*Au[1]+uu[0]*c2*Au[2])/k1;
899  inv_J[3][1] = Au[0]*Au[1]*(c2+uu[1])*c3/k1;
900  inv_J[3][2] = Au[0]*Au[2]*(c3+uu[2])*c2/k1;
901  inv_J[3][3] = Au[0]*c3*c2/k1;
902  inv_J[3][4] = 0.5*Au[0]*Au[1]*c3/k1;
903  inv_J[3][5] = 0.5*Au[0]*c2*Au[2]/k1;
904 
905  inv_J[4][0] = -Au[0]*Au[1]*(c1-uu[0])*c3/k2;
906  inv_J[4][1] = Au[1]*(Au[1]*c1*c3-c1*uu[1]*Au[2]-c3*uu[1]*Au[0])/k2;
907  inv_J[4][2] = Au[2]*Au[1]*(c3+uu[2])*c1/k2;
908  inv_J[4][3] = Au[1]*c1*c3/k2;
909  inv_J[4][4] = -0.5*Au[1]*(c1*Au[2]+Au[0]*c3)/k2;
910  inv_J[4][5] = 0.5*Au[2]*Au[1]*c1/k2;
911 
912  inv_J[5][0] = -Au[0]*Au[2]*(c1-uu[0])*c2/k3;
913  inv_J[5][1] = Au[2]*Au[1]*(c2+uu[1])*c1/k3;
914  inv_J[5][2] = Au[2]*(Au[2]*c1*c2-c1*uu[2]*Au[1]-c2*uu[2]*Au[0])/k3;
915  inv_J[5][3] = Au[2]*c1*c2/k3;
916  inv_J[5][4] = 0.5*Au[2]*Au[1]*c1/k3;
917  inv_J[5][5] = -0.5*Au[2]*(Au[1]*c1+Au[0]*c2)/k3;
918 
919  // Solve the system by multiplying the Jacobian with the vector f:
920  // g = (inv_J)*f
921  for (int j=0; j<6; j++)
922  {
923  tmp[j] =0.0;
924  g[j] = 0.0;
925  }
926 
927  for (int j=0; j<6; j++)
928  {
929  for (int i=0; i<6; i++)
930  {
931  tmp[j] = inv_J[j][i]*f[i];
932  g[j] += tmp[j];
933  }
934  }
935 
936  // Update the solution: x_new = x_old - dx
937  uu[0] = uu[0] - g[0];
938  uu[1] = uu[1] - g[1];
939  uu[2] = uu[2] - g[2];
940  Au[0] = Au[0] - g[3];
941  Au[1] = Au[1] - g[4];
942  Au[2] = Au[2] - g[5];
943 
944  // Check if the error of the solution is smaller than Tol
945  if ((g[0]*g[0] + g[1]*g[1] + g[2]*g[2] + g[3]*g[3]+ g[4]*g[4] + g[5]*g[5]) < Tol)
946  {
947  proceed = 0;
948  }
949 
950  // Check if solver converges
951  if (iter >= MAX_ITER)
952  {
953  ASSERTL0(false,"Riemann solver for Merging Flow did not converge");
954  }
955 
956  }
957  }
958 
959  /**
960  * Solves the Riemann problem at an interdomain junction by
961  * assuming subsonic flow at both sides of the boundary and
962  * by applying conservation of mass and continuity of the
963  * total pressure \f$ \frac{p}{rho} + \frac{u^{2}}{2}. \f$
964  * The other 2 missing equations come from the characteristic
965  * variables. For further information see "Pulse
966  * WavePropagation in the human vascular system" Section 3.4.
967  */
968  void PulseWaveSystem::JunctionRiemann(Array<OneD, NekDouble> &Au, Array<OneD, NekDouble> &uu,
969  Array<OneD, NekDouble> &beta, Array<OneD, NekDouble> &A_0)
970  {
971  NekDouble rho = m_rho;
972  Array<OneD, NekDouble> W(2);
973  Array<OneD, NekDouble> W_Au(2);
974  Array<OneD, NekDouble> P_Au(2);
975  Array<OneD, NekDouble> f(4);
976  Array<OneD, NekDouble> g(4);
977  Array<OneD, NekDouble> tmp(4);
978  Array<OneD, Array<OneD, NekDouble> > inv_J(4);
979  for (int i=0; i<4; i++)
980  {
981  inv_J[i] = Array<OneD, NekDouble> (4);
982  }
983  NekDouble k = 0.0;
984  NekDouble k1 = 0.0;
985  NekDouble k2 = 0.0;
986 
987  int proceed = 1;
988  int iter = 0;
989  int MAX_ITER = 7;
990  NekDouble Tol = 1.0e-10;
991 
992  // Calculated from input
993  W[0] = uu[0] + 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
994  W[1] = uu[1] - 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
995 
996  while((proceed) && (iter < MAX_ITER))
997  {
998  iter = iter+1;
999 
1000  // Calculate the constraint vector, 4 equations:
1001  // 2 characteristic variables, mass conservation,
1002  // total pressure
1003  W_Au[0] = 4*sqrt(beta[0]/(2*rho))*(sqrt(sqrt(Au[0])) - sqrt(sqrt(A_0[0])));
1004  W_Au[1] = 4*sqrt(beta[1]/(2*rho))*(sqrt(sqrt(Au[1])) - sqrt(sqrt(A_0[1])));
1005 
1006  P_Au[0] = beta[0]*(sqrt(Au[0]) - sqrt(A_0[0]));
1007  P_Au[1] = beta[1]*(sqrt(Au[1]) - sqrt(A_0[1]));
1008 
1009  f[0] = uu[0] + W_Au[0] - W[0];
1010  f[1] = uu[1] - W_Au[1] - W[1];
1011  f[2] = Au[0]*uu[0] - Au[1]*uu[1];
1012  f[3] = uu[0]*uu[0] + 2.0/rho*P_Au[0] - uu[1]*uu[1] - 2.0/rho*P_Au[1];
1013 
1014  // Calculate the wave speed at each vessel
1015  NekDouble cl = sqrt(beta[0]/(2*rho))*sqrt(sqrt(Au[0]));
1016  NekDouble cr = sqrt(beta[1]/(2*rho))*sqrt(sqrt(Au[1]));
1017 
1018  // Inverse Jacobian matrix J(x[n])^(-1), is already inverted here analytically
1019  k = (cl*Au[1]+Au[0]*cr);
1020  k1 = (cl-uu[0])*k;
1021  inv_J[0][0] = (Au[1]*cl*cl-cr*uu[0]*Au[0])/k1;
1022  inv_J[0][1] = Au[1]*(cr-uu[1])*cl/k1;
1023  inv_J[0][2] = cl*cr/k1;
1024  inv_J[0][3] = -0.5*cl*Au[1]/k1;
1025 
1026  k2 = (cr+uu[1])*k;
1027  inv_J[1][0] = Au[0]*(cl+uu[0])*cr/k2;
1028  inv_J[1][1] = (cl*uu[1]*Au[1]+cr*cr*Au[0])/k2;
1029  inv_J[1][2] = -cl*cr/k2;
1030  inv_J[1][3] = -0.5*Au[0]*cr/k2;
1031 
1032  inv_J[2][0] = Au[0]*(Au[0]*cr-uu[0]*Au[1])/k1;
1033  inv_J[2][1] = -Au[0]*Au[1]*(cr-uu[1])/k1;
1034  inv_J[2][2] = -Au[0]*cr/k1;
1035  inv_J[2][3] = 0.5*Au[1]*Au[0]/k1;
1036 
1037  inv_J[3][0] = Au[0]*Au[1]*(cl+uu[0])/k2;
1038  inv_J[3][1] = -Au[1]*(cl*Au[1]+uu[1]*Au[0])/k2;
1039  inv_J[3][2] = -cl*Au[1]/k2;
1040  inv_J[3][3] = -0.5*Au[1]*Au[0]/k2;
1041 
1042  // Solve the system by multiplying the Jacobian with the vector f:
1043  // g = (inv_J)*f
1044  for (int j=0; j<4; j++)
1045  {
1046  tmp[j] =0.0;
1047  g[j] = 0.0;
1048  }
1049 
1050  for (int j=0; j<4; j++)
1051  {
1052  for (int i=0; i<4; i++)
1053  {
1054  tmp[j] = inv_J[j][i]*f[i];
1055  g[j] += tmp[j];
1056  }
1057  }
1058 
1059  // Update solution: x_new = x_old - dx
1060  uu[0] -= g[0];
1061  uu[1] -= g[1];
1062  Au[0] -= g[2];
1063  Au[1] -= g[3];
1064 
1065  // Check if the error of the solution is smaller than Tol.
1066  if((g[0]*g[0] + g[1]*g[1] + g[2]*g[2] + g[3]*g[3]) < Tol)
1067  proceed = 0;
1068  }
1069 
1070  if(iter >= MAX_ITER)
1071  {
1072  ASSERTL0(false,"Riemann solver for Junction did not converge");
1073  }
1074  }
1075 
1076 
1077 
1078  /**
1079  * Writes the .fld file at the end of the simulation. Similar to the normal
1080  * v_Output however the Multidomain output has to be prepared.
1081  */
1082  void PulseWaveSystem::v_Output(void)
1083  {
1084  /**
1085  * Write the field data to file. The file is named according to the session
1086  * name with the extension .fld appended.
1087  */
1088  std::string outname = m_sessionName + ".fld";
1089 
1090  WriteVessels(outname);
1091  }
1092 
1093 
1094  /**
1095  * Writes the .fld file at the end of the simulation. Similar to the normal
1096  * v_Output however the Multidomain output has to be prepared.
1097  */
1098  void PulseWaveSystem::CheckPoint_Output(const int n)
1099  {
1100  std::stringstream outname;
1101  outname << m_sessionName << "_" << n << ".chk";
1102 
1103  WriteVessels(outname.str());
1104  }
1105 
1106 
1107  /**
1108  * Writes the field data to a file with the given filename.
1109  * @param outname Filename to write to.
1110  */
1111  void PulseWaveSystem::WriteVessels(const std::string &outname)
1112  {
1113 
1114  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1115  std::vector<std::string> variables = m_session->GetVariables();
1116 
1117  for(int n = 0; n < m_nDomains; ++n)
1118  {
1119  m_vessels[n*m_nVariables]->GetFieldDefinitions(FieldDef);
1120  }
1121 
1122  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1123 
1124  int nFieldDefPerDomain = FieldDef.size()/m_nDomains;
1125  int cnt;
1126  // Copy Data into FieldData and set variable
1127  for(int n = 0; n < m_nDomains; ++n)
1128  {
1129  for(int j = 0; j < m_nVariables; ++j)
1130  {
1131  for(int i = 0; i < nFieldDefPerDomain; ++i)
1132  {
1133  cnt = n*nFieldDefPerDomain+i;
1134  FieldDef[cnt]->m_fields.push_back(variables[j]);
1135  m_vessels[n*m_nVariables]->AppendFieldData(FieldDef[cnt], FieldData[cnt], m_vessels[n*m_nVariables+j]->UpdateCoeffs());
1136  }
1137  }
1138  }
1139 
1140  // Update time in field info if required
1141  if(m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
1142  {
1143  m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
1144  }
1145 
1146  //LibUtilities::CombineFields(FieldDef, FieldData);
1147 
1148  LibUtilities::Write(outname, FieldDef, FieldData, m_fieldMetaDataMap);
1149  }
1150 
1151  /* Compute the error in the L2-norm
1152  * @param field The field to compare.
1153  * @param exactsoln The exact solution to compare with.
1154  * @param Normalised Normalise L2-error.
1155  * @returns Error in the L2-norm.
1156  */
1157  NekDouble PulseWaveSystem::v_L2Error(unsigned int field,
1158  const Array<OneD, NekDouble> &exactsoln,
1159  bool Normalised)
1160  {
1161  NekDouble L2error = 0.0;
1162  NekDouble L2error_dom;
1163  NekDouble Vol = 0.0;
1164 
1165  if(m_NumQuadPointsError == 0)
1166  {
1167  for (int omega = 0; omega < m_nDomains; omega++)
1168  {
1169  int vesselid = field + omega*m_nVariables;
1170 
1171  if(m_vessels[vesselid]->GetPhysState() == false)
1172  {
1173  m_vessels[vesselid]->BwdTrans(m_vessels[vesselid]->GetCoeffs(),
1174  m_vessels[vesselid]->UpdatePhys());
1175  }
1176 
1177  if(exactsoln.num_elements())
1178  {
1179  L2error_dom = m_vessels[vesselid]->L2(
1180  m_vessels[vesselid]->GetPhys(),
1181  exactsoln);
1182  }
1183  else if (m_session->DefinesFunction("ExactSolution"))
1184  {
1185  Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
1186 
1187  LibUtilities::EquationSharedPtr vEqu
1188  = m_session->GetFunction("ExactSolution",field,omega);
1189  EvaluateFunction(m_session->GetVariable(field),exactsoln,"ExactSolution",
1190  m_time);
1191 
1192  L2error_dom = m_vessels[vesselid]->L2(
1193  m_vessels[vesselid]->GetPhys(),
1194  exactsoln);
1195 
1196  }
1197  else
1198  {
1199  L2error_dom = m_vessels[vesselid]->L2(
1200  m_vessels[vesselid]->GetPhys());
1201  }
1202 
1203  L2error += L2error_dom*L2error_dom;
1204 
1205  if(Normalised == true)
1206  {
1207  Array<OneD, NekDouble> one(m_vessels[vesselid]->GetNpoints(), 1.0);
1208 
1209  Vol += m_vessels[vesselid]->PhysIntegral(one);
1210  }
1211 
1212  }
1213  }
1214  else
1215  {
1216  ASSERTL0(false,"Not set up");
1217  }
1218 
1219 
1220  if(Normalised == true)
1221  {
1222  m_comm->AllReduce(Vol, LibUtilities::ReduceSum);
1223 
1224  L2error = sqrt(L2error/Vol);
1225  }
1226  else
1227  {
1228  L2error = sqrt(L2error);
1229  }
1230 
1231  return L2error;
1232  }
1233 
1234 
1235  /**
1236  * Compute the error in the L_inf-norm
1237  * @param field The field to compare.
1238  * @param exactsoln The exact solution to compare with.
1239  * @returns Error in the L_inft-norm.
1240  */
1241  NekDouble PulseWaveSystem::v_LinfError(unsigned int field,
1242  const Array<OneD, NekDouble> &exactsoln)
1243  {
1244  NekDouble LinferrorDom, Linferror = -1.0;
1245 
1246  for (int omega = 0; omega < m_nDomains; omega++)
1247  {
1248  int vesselid = field + omega*m_nVariables;
1249 
1250  if(m_NumQuadPointsError == 0)
1251  {
1252  if(m_vessels[vesselid]->GetPhysState() == false)
1253  {
1254  m_vessels[vesselid]->BwdTrans(m_vessels[vesselid]->GetCoeffs(),
1255  m_vessels[vesselid]->UpdatePhys());
1256  }
1257 
1258  if(exactsoln.num_elements())
1259  {
1260  LinferrorDom = m_vessels[vesselid]->Linf(
1261  m_vessels[vesselid]->GetPhys(),
1262  exactsoln);
1263  }
1264  else if (m_session->DefinesFunction("ExactSolution"))
1265  {
1266  Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
1267 
1268  EvaluateFunction(m_session->GetVariable(field),exactsoln,"ExactSolution",
1269  m_time);
1270 
1271  LinferrorDom = m_vessels[vesselid]->Linf(
1272  m_vessels[vesselid]->GetPhys(),
1273  exactsoln);
1274  }
1275  else
1276  {
1277  LinferrorDom = 0.0;
1278  }
1279 
1280  Linferror = (Linferror > LinferrorDom)? Linferror:LinferrorDom;
1281 
1282  }
1283  else
1284  {
1285  ASSERTL0(false,"ErrorExtraPoints not allowed for this solver");
1286  }
1287  }
1288  return Linferror;
1289  }
1290 
1291  void PulseWaveSystem::CalcCharacteristicVariables(int omega)
1292  {
1293  int nq = m_vessels[omega]->GetTotPoints();
1294  Array<OneD, NekDouble> A = m_vessels[omega]->UpdatePhys();
1295  Array<OneD, NekDouble> u = m_vessels[omega+1]->UpdatePhys();
1296  Array<OneD, NekDouble> c(nq);
1297 
1298  //Calc 4*c
1299  Vmath::Vsqrt(nq,A,1,c,1);
1300  Vmath::Vmul(nq,m_beta[omega],1,c,1,c,1);
1301  Vmath::Smul(nq,16.0/(2*m_rho),c,1,c,1);
1302  Vmath::Vsqrt(nq,c,1,c,1);
1303 
1304  // Characteristics
1305  Vmath::Vadd(nq,u,1,c,1,A,1);
1306  Vmath::Vsub(nq,u,1,c,1,u,1);
1307  }
1308 
1309 
1310 
1311 }
Generated on Sun Mar 15 2015 00:11:34 for Nektar++ by   doxygen 1.8.1.2