51 Loki::NoDestroy > Type;
52 return Type::Instance();
75 m_substeps = pSession->GetParameter(
"Substeps");
81 std::set<enum LibUtilities::ShapeType> s;
82 for (
unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
84 s.insert(m_field->GetExp(i)->DetShapeType());
99 m_nq = pField->GetNcoeffs();
100 int order = m_field->GetExp(0)->GetBasis(0)->GetNumModes();
120 m_nq = pField->GetTotPoints();
128 void CellModel::Initialise()
130 ASSERTL1(m_nvar > 0,
"Cell model must have at least 1 variable.");
134 for (
unsigned int i = 0; i < m_nvar; ++i)
140 for (
unsigned int i = 0; i < m_gates.size(); ++i)
145 if (m_session->DefinesFunction(
"CellModelInitialConditions"))
151 v_SetInitialConditions();
163 void CellModel::TimeIntegrate(
170 int nvar = inarray.num_elements();
177 if (!m_nodalTmp.num_elements())
180 for (
unsigned int k = 0; k < nvar; ++k)
189 for (
unsigned int k = 0; k < nvar; ++k)
191 for (
unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
193 phys_offset = m_field->GetPhys_Offset(i);
194 coef_offset = m_field->GetCoeff_Offset(i);
197 m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, tmpCoeffs);
198 m_nodalTri->ModalToNodal(tmpCoeffs, tmp=m_nodalTmp[k]+coef_offset);
202 m_field->GetExp(0)->FwdTrans(inarray[k] + phys_offset, tmpCoeffs);
203 m_nodalTet->ModalToNodal(tmpCoeffs, tmp=m_nodalTmp[k]+coef_offset);
217 NekDouble delta_t = (time - m_lastTime)/m_substeps;
221 for (
unsigned int i = 0; i < m_substeps - 1; ++i)
223 Update(m_cellSol, m_wsp, time);
225 Vmath::Svtvp(m_nq, delta_t, m_wsp[0], 1, m_cellSol[0], 1, m_cellSol[0], 1);
227 for (
unsigned int j = 0; j < m_concentrations.size(); ++j)
229 Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
232 for (
unsigned int j = 0; j < m_gates.size(); ++j)
234 Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
235 Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
236 Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
237 Vmath::Vvtvp(m_nq, m_cellSol[m_gates[j]], 1, m_gates_tau[j], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
242 Update(m_cellSol, m_wsp, time);
250 for (
unsigned int k = 0; k < nvar; ++k)
252 for (
unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
254 int phys_offset = m_field->GetPhys_Offset(i);
255 int coef_offset = m_field->GetCoeff_Offset(i);
258 m_nodalTri->NodalToModal(m_wsp[k]+coef_offset, tmpCoeffs);
259 m_field->GetExp(0)->BwdTrans(tmpCoeffs, tmp=outarray[k] + phys_offset);
263 m_nodalTet->NodalToModal(m_wsp[k]+coef_offset, tmpCoeffs);
264 m_field->GetExp(0)->BwdTrans(tmpCoeffs, tmp=outarray[k] + phys_offset);
275 for (
unsigned int j = 0; j < m_concentrations.size(); ++j)
277 Vmath::Svtvp(m_nq, delta_t, m_wsp[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1, m_cellSol[m_concentrations[j]], 1);
281 for (
unsigned int j = 0; j < m_gates.size(); ++j)
283 Vmath::Sdiv(m_nq, -delta_t, m_gates_tau[j], 1, m_gates_tau[j], 1);
284 Vmath::Vexp(m_nq, m_gates_tau[j], 1, m_gates_tau[j], 1);
285 Vmath::Vsub(m_nq, m_cellSol[m_gates[j]], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
286 Vmath::Vvtvp(m_nq, m_cellSol[m_gates[j]], 1, m_gates_tau[j], 1, m_wsp[m_gates[j]], 1, m_cellSol[m_gates[j]], 1);
294 ASSERTL0(idx < m_nvar,
"Index out of range for cell model.");
301 for (
unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
303 int coef_offset = m_field->GetCoeff_Offset(i);
306 m_nodalTri->NodalToModal(m_cellSol[idx]+coef_offset, tmp=outarray+coef_offset);
310 m_nodalTet->NodalToModal(m_cellSol[idx]+coef_offset, tmp=outarray+coef_offset);
316 m_field->FwdTrans_IterPerExp(m_cellSol[idx], outarray);
324 return m_cellSol[idx];
327 void CellModel::LoadCellModel()
329 const bool root = (m_session->GetComm()->GetRank() == 0);
330 const std::string fncName =
"CellModelInitialConditions";
331 const int nvar = m_cellSol[0].num_elements();
341 cout <<
"Cell model initial conditions: " << endl;
345 std::set<std::string> filelist;
346 for (j = 1; j < nvar; ++j)
349 varName = GetCellVarName(j);
351 if (m_session->GetFunctionType(fncName, varName) ==
354 filelist.insert(m_session->GetFunctionFilename(fncName,
360 typedef std::vector<LibUtilities::FieldDefinitionsSharedPtr> FDef;
361 typedef std::vector<std::vector<NekDouble> > FData;
362 std::map<std::string, FDef> FieldDef;
363 std::map<std::string, FData> FieldData;
366 std::set<std::string>::const_iterator setIt;
370 for (setIt = filelist.begin(); setIt != filelist.end(); ++setIt)
374 cout <<
" - Reading file: " << *setIt << endl;
376 FieldDef[*setIt] = FDef(0);
377 FieldData[*setIt] = FData(0);
378 fld->Import(*setIt, FieldDef[*setIt], FieldData[*setIt],
383 iter = fieldMetaDataMap.find(
"Time");
384 if(iter != fieldMetaDataMap.end())
386 m_lastTime = boost::lexical_cast<
NekDouble>(iter->second);
392 for(j = 1; j < m_cellSol.num_elements(); ++j)
395 varName = GetCellVarName(j);
398 if (m_session->GetFunctionType(fncName, varName) ==
401 const std::string file =
402 m_session->GetFunctionFilename(fncName, varName);
406 cout <<
" - Field " << varName <<
": from file "
411 for(
int i = 0; i < FieldDef[file].size(); ++i)
413 m_field->ExtractDataToCoeffs(FieldDef[file][i],
423 for (
unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
425 int coef_offset = m_field->GetCoeff_Offset(i);
426 if (m_field->GetExp(0)->DetShapeType() ==
429 m_nodalTri->ModalToNodal(coeffs+coef_offset,
430 tmp=m_cellSol[j]+coef_offset);
434 m_nodalTet->ModalToNodal(coeffs+coef_offset,
435 tmp=m_cellSol[j]+coef_offset);
441 m_field->BwdTrans(coeffs, m_cellSol[j]);
444 else if (m_session->GetFunctionType(fncName, varName) ==
448 m_session->GetFunction(fncName, varName);
452 cout <<
" - Field " << varName <<
": "
453 << equ->GetExpression() << endl;
456 const unsigned int nphys = m_field->GetNpoints();
460 m_field->GetCoords(x0,x1,x2);
467 equ->Evaluate(x0, x1, x2, phys);
468 for (
unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
470 int phys_offset = m_field->GetPhys_Offset(i);
471 int coef_offset = m_field->GetCoeff_Offset(i);
472 if (m_field->GetExp(0)->DetShapeType() ==
475 m_field->GetExp(0)->FwdTrans(
476 phys + phys_offset, tmpCoeffs);
477 m_nodalTri->ModalToNodal(
479 tmp = m_cellSol[j] + coef_offset);
483 m_field->GetExp(0)->FwdTrans(
484 phys + phys_offset, tmpCoeffs);
485 m_nodalTet->ModalToNodal(
487 tmp = m_cellSol[j] + coef_offset);
493 equ->Evaluate(x0, x1, x2, m_cellSol[j]);
LibUtilities::NekFactory< std::string, CellModel, const LibUtilities::SessionReaderSharedPtr &, const MultiRegions::ExpListSharedPtr & > CellModelFactory
Datatype of the NekFactory used to instantiate classes derived from the EquationSystem class...
#define ASSERTL0(condition, msg)
Principle Modified Functions .
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Principle Modified Functions .
void Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
std::map< std::string, std::string > FieldMetaDataMap
Gauss Radau pinned at x=-1, .
Principle Modified Functions .
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
void Vexp(int n, const T *x, const int incx, T *y, const int incy)
Defines a specification for a set of points.
boost::shared_ptr< FieldIO > FieldIOSharedPtr
3D Evenly-spaced points on a Tetrahedron
boost::shared_ptr< Equation > EquationSharedPtr
CellModelFactory & GetCellModelFactory()
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.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Gauss Radau pinned at x=-1, .
2D Evenly-spaced points on a Triangle
void Zero(int n, T *x, const int incx)
Zero vector.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Describes the specification for a Basis.
1D Gauss-Lobatto-Legendre quadrature points
Provides a generic Factory class.