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");
105 const std::vector<SpatialDomains::CompositeMap> domain =
m_graph->GetDomain();
111 bool SetToOneSpaceDimension =
true;
113 if(
m_session->DefinesCmdLineArgument(
"SetToOneSpaceDimension"))
115 std::string cmdline =
m_session->GetCmdLineArgument<std::string>(
"SetToOneSpaceDimension");
116 if(boost::to_upper_copy(cmdline) ==
"FALSE")
118 SetToOneSpaceDimension =
false;
130 SetToOneSpaceDimension);
153 m_vessels[n]->SetCoeffsArray(coeffs);
154 m_vessels[n]->SetPhysArray(phys);
156 int cnt = m_vessels[n]->GetNcoeffs();
157 int cnt1 = m_vessels[n]->GetTotPoints();
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();
217 for (
int omega = 0; omega <
m_nDomains; omega++)
219 nq = m_vessels[2*omega]->GetNpoints();
234 m_fields[0]->ExtractTracePhys(m_A_0[omega],m_A_0_trace[omega]);
237 if(SetToOneSpaceDimension)
242 int nelmt_trace = trace->GetExpSize();
246 for(
int i = 0 ; i < nelmt_trace; ++i)
248 normals[i] = m_trace_fwd_normal[omega]+i;
253 trace->GetExp(0)->GetGeom()->SetCoordim(1);
255 trace->GetNormals(normals);
265 map<int,std::vector<InterfacePointShPtr> > VidToDomain;
266 map<int,std::vector<InterfacePointShPtr> >
::iterator iter;
272 for(
int omega = 0; omega <
m_nDomains; ++omega)
276 for(
int i = 0; i < 2; ++i)
281 int vid =
m_vessels[vesselID]->UpdateBndCondExpansion(i)->GetExp(0)->GetGeom()->GetVid(0);
282 cout<<
"Shared vertex id: "<<vid<<endl;
288 for(
int n = 0; n <
m_vessels[vesselID]->GetExpSize(); ++n)
290 for(
int p = 0;
p < 2; ++
p)
294 int eid =
m_vessels[vesselID]->GetTraceMap()->GetElmtToTrace()[n][
p]->GetElmtId();
296 int tid =
m_vessels[vesselID]->GetTrace()->GetCoeff_Offset(eid);
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;
316 VidToDomain[vid].push_back(Ipt);
319 m_vessels[vesselID]->GetBndConditions()[i]
322 m_vessels[vesselID+1]->GetBndConditions()[i]
329 for(iter = VidToDomain.begin(); iter != VidToDomain.end(); ++iter)
331 if(iter->second.size() == 2)
335 else if(iter->second.size() == 3)
344 for(
int i = 0; i < 3; ++i)
346 if(iter->second[i]->m_elmtVert == 0)
360 if(iter->second[0]->m_elmtVert == 1)
369 if(iter->second[1]->m_elmtVert == 1)
372 iter->second[0] = iter->second[1];
375 else if (iter->second[2]->m_elmtVert == 1)
378 iter->second[0] = iter->second[2];
387 if(iter->second[0]->m_elmtVert == 0)
396 if(iter->second[1]->m_elmtVert == 0)
399 iter->second[0] = iter->second[1];
402 else if (iter->second[2]->m_elmtVert == 0)
405 iter->second[0] = iter->second[2];
415 ASSERTL0(
false,
"Unknown junction type");
431 if (
m_session->GetComm()->GetRank() == 0)
433 cout <<
"Initial Conditions:" << endl;
438 for (
int omega = 0; omega <
m_nDomains; omega++)
443 if (
m_session->GetComm()->GetRank() == 0)
445 cout <<
"Subdomain = " <<omega<<endl;
503 cout <<
"Steps: " << n+1
505 <<
"\t Time-step: " <<
m_timestep <<
"\t" << endl;
514 for (
int omega = 0; omega <
m_nDomains; omega++)
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();
532 cout <<
"Time-integration timing : "
533 << IntegrationTime <<
" s" << endl << endl;
543 int omega, traceId, eid, vert, phys_offset, vesselID;
548 traceId = I->m_traceId;
550 vert = I->m_elmtVert;
553 phys_offset =
m_vessels[vesselID]->GetPhys_Offset(eid);
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]);
584 for(
int i = 0; i < 3; ++i)
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]);
608 for(
int i = 0; i < 3; ++i)
621 Au[0],uu[0],beta[0],A_0[0]);
623 Au[1],uu[1],beta[1],A_0[1]);
628 for(
int i = 0; i < 2; ++i)
664 for (
int i=0; i<6; i++)
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])));
686 while ((proceed) && (iter < MAX_ITER))
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])));
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]));
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];
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]));
715 k = c1*Au[1]*c3+Au[0]*c3*c2+Au[2]*c1*c2;
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;
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;
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;
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;
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;
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;
764 for (
int j=0; j<6; j++)
770 for (
int j=0; j<6; j++)
772 for (
int i=0; i<6; i++)
774 tmp[j] = inv_J[j][i]*f[i];
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];
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)
794 if (iter >= MAX_ITER)
796 ASSERTL0(
false,
"Riemann solver for Bifurcation did not converge");
823 for (
int i=0; i<6; i++)
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])));
846 while ((proceed) && (iter < MAX_ITER))
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])));
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]));
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];
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]));
874 k = c1*Au[1]*c3+Au[0]*c3*c2+Au[2]*c1*c2;
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;
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;
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;
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;
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;
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;
922 for (
int j=0; j<6; j++)
928 for (
int j=0; j<6; j++)
930 for (
int i=0; i<6; i++)
932 tmp[j] = inv_J[j][i]*f[i];
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];
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)
952 if (iter >= MAX_ITER)
954 ASSERTL0(
false,
"Riemann solver for Merging Flow did not converge");
980 for (
int i=0; i<4; i++)
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])));
997 while((proceed) && (iter < MAX_ITER))
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])));
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]));
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];
1016 NekDouble cl = sqrt(beta[0]/(2*rho))*sqrt(sqrt(Au[0]));
1017 NekDouble cr = sqrt(beta[1]/(2*rho))*sqrt(sqrt(Au[1]));
1020 k = (cl*Au[1]+Au[0]*cr);
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;
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;
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;
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;
1045 for (
int j=0; j<4; j++)
1051 for (
int j=0; j<4; j++)
1053 for (
int i=0; i<4; i++)
1055 tmp[j] = inv_J[j][i]*f[i];
1067 if((g[0]*g[0] + g[1]*g[1] + g[2]*g[2] + g[3]*g[3]) < Tol)
1071 if(iter >= MAX_ITER)
1073 ASSERTL0(
false,
"Riemann solver for Junction did not converge");
1101 std::stringstream outname;
1115 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1116 std::vector<std::string> variables =
m_session->GetVariables();
1123 std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
1125 int nFieldDefPerDomain = FieldDef.size()/
m_nDomains;
1132 for(
int i = 0; i < nFieldDefPerDomain; ++i)
1134 cnt = n*nFieldDefPerDomain+i;
1135 FieldDef[cnt]->m_fields.push_back(variables[j]);
1168 for (
int omega = 0; omega <
m_nDomains; omega++)
1172 if(
m_vessels[vesselid]->GetPhysState() ==
false)
1178 if(exactsoln.num_elements())
1184 else if (
m_session->DefinesFunction(
"ExactSolution"))
1189 =
m_session->GetFunction(
"ExactSolution",field,omega);
1204 L2error += L2error_dom*L2error_dom;
1206 if(Normalised ==
true)
1210 Vol +=
m_vessels[vesselid]->PhysIntegral(one);
1221 if(Normalised ==
true)
1225 L2error = sqrt(L2error/Vol);
1229 L2error = sqrt(L2error);
1245 NekDouble LinferrorDom, Linferror = -1.0;
1247 for (
int omega = 0; omega <
m_nDomains; omega++)
1253 if(
m_vessels[vesselid]->GetPhysState() ==
false)
1259 if(exactsoln.num_elements())
1261 LinferrorDom =
m_vessels[vesselid]->Linf(
1265 else if (
m_session->DefinesFunction(
"ExactSolution"))
1272 LinferrorDom =
m_vessels[vesselid]->Linf(
1281 Linferror = (Linferror > LinferrorDom)? Linferror:LinferrorDom;
1286 ASSERTL0(
false,
"ErrorExtraPoints not allowed for this solver");
1294 int nq =
m_vessels[omega]->GetTotPoints();
Array< OneD, int > m_fieldPhysOffset
Array< OneD, Array< OneD, NekDouble > > m_beta_trace
#define ASSERTL0(condition, msg)
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...
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)
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
virtual void v_InitObject()
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
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.
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.
void SetUpDomainInterfaces(void)
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.
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.
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 ...
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.
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 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.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
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.
boost::shared_ptr< InterfacePoint > InterfacePointShPtr