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