Merge branch 'ticket/111-cell-model-restart' of localhost:nektar

library/LibUtilities/BasicUtils/SessionReader.cpp

@@ -809,10 +809,13 @@ namespace Nektar
             FunctionVariableMap::const_iterator it2;
             std::string vName = boost::to_upper_copy(pName);
 
-            if ((it1 = m_functions.find(vName))     != m_functions.end() &&
-                (it2 = it1->second.find(pVariable)) != it1->second.end())
+            // Check function exists
+            if ((it1 = m_functions.find(vName))     != m_functions.end())
             {
-                return true;
+                bool varExists =
+                    (it2 = it1->second.find(pVariable)) != it1->second.end() ||
+                    (it2 = it1->second.find("*")) != it1->second.end();
+                return varExists;
             }
             return false;
         }
@@ -826,18 +829,33 @@ namespace Nektar
             const std::string &pVariable) const
         {
             FunctionMap::const_iterator it1;
-            FunctionVariableMap::const_iterator it2;
+            FunctionVariableMap::const_iterator it2, it3;
             std::string vName = boost::to_upper_copy(pName);
 
             ASSERTL0((it1 = m_functions.find(vName)) != m_functions.end(),
                      std::string("No such function '") + pName
                      + std::string("' has been defined in the session file."));
-            ASSERTL0((it2 = it1->second.find(pVariable)) != it1->second.end(),
+
+            // Check for specific and wildcard definitions
+            bool specific = (it2 = it1->second.find(pVariable)) !=
+                            it1->second.end();
+            bool wildcard = (it3 = it1->second.find("*")) !=
+                            it1->second.end();
+
+            // Check function is defined somewhere
+            ASSERTL0(specific || wildcard,
                      std::string("No such variable '") + pVariable
                      + std::string("' defined for function '") + pName
                      + std::string("' in session file."));
+
+            // If not specific, must be wildcard
+            if (!specific)
+            {
+                it2 = it3;
+            }
+
             ASSERTL0((it2->second.m_type == eFunctionTypeExpression),
-                     std::string("Function is defined by a file."));
+                    std::string("Function is defined by a file."));
             return it2->second.m_expression;
         }
 
@@ -862,17 +880,32 @@ namespace Nektar
             const std::string &pVariable) const
         {
             FunctionMap::const_iterator it1;
-            FunctionVariableMap::const_iterator it2;
+            FunctionVariableMap::const_iterator it2, it3;
             std::string vName = boost::to_upper_copy(pName);
 
             it1 = m_functions.find(vName);
             ASSERTL0 (it1 != m_functions.end(),
                       std::string("Function '") + pName
                       + std::string("' not found."));
-            ASSERTL0 ((it2 = it1->second.find(pVariable)) != it1->second.end(),
-                    std::string("No such variable '") + pVariable
-                    + std::string("' defined for function '") + pName
-                    + std::string("' in session file."));
+
+            // Check for specific and wildcard definitions
+            bool specific = (it2 = it1->second.find(pVariable)) !=
+                            it1->second.end();
+            bool wildcard = (it3 = it1->second.find("*")) !=
+                            it1->second.end();
+
+            // Check function is defined somewhere
+            ASSERTL0(specific || wildcard,
+                     std::string("No such variable '") + pVariable
+                     + std::string("' defined for function '") + pName
+                     + std::string("' in session file."));
+
+            // If not specific, must be wildcard
+            if (!specific)
+            {
+                it2 = it3;
+            }
+
             return it2->second.m_type;
         }
 
@@ -897,17 +930,32 @@ namespace Nektar
             const std::string &pVariable) const
         {
             FunctionMap::const_iterator it1;
-            FunctionVariableMap::const_iterator it2;
+            FunctionVariableMap::const_iterator it2, it3;
             std::string vName = boost::to_upper_copy(pName);
 
             it1 = m_functions.find(vName);
             ASSERTL0 (it1 != m_functions.end(),
                       std::string("Function '") + pName
                       + std::string("' not found."));
-            ASSERTL0 ((it2 = it1->second.find(pVariable)) != it1->second.end(),
-                    std::string("No such variable '") + pVariable
-                    + std::string("' defined for function '") + pName
-                    + std::string("' in session file."));
+
+            // Check for specific and wildcard definitions
+            bool specific = (it2 = it1->second.find(pVariable)) !=
+                            it1->second.end();
+            bool wildcard = (it3 = it1->second.find("*")) !=
+                            it1->second.end();
+
+            // Check function is defined somewhere
+            ASSERTL0(specific || wildcard,
+                     std::string("No such variable '") + pVariable
+                     + std::string("' defined for function '") + pName
+                     + std::string("' in session file."));
+
+            // If not specific, must be wildcard
+            if (!specific)
+            {
+                it2 = it3;
+            }
+
             return it2->second.m_filename;
         }
 
@@ -1679,11 +1727,16 @@ namespace Nektar
                     FunctionVariableDefinition funcDef;
                     std::string conditionType = variable->Value();
 
-                    // All function variables must specify VAR
-                    ASSERTL0(variable->Attribute("VAR"),
-                             "Attribute VAR expected for function '"
-                             + functionStr + "'.");
-                    std::string variableStr = variable->Attribute("VAR");
+                    // If no var is specified, assume wildcard
+                    std::string variableStr;
+                    if (!variable->Attribute("VAR"))
+                    {
+                        variableStr = "*";
+                    }
+                    else
+                    {
+                        variableStr = variable->Attribute("VAR");
+                    }
 
                     // Parse list of variables
                     std::vector<std::string> variableList;

library/MultiRegions/ExpList.cpp

@@ -2131,13 +2131,8 @@ namespace Nektar
                 offset += datalen;
             }
             
-            if(i == fielddef->m_fields.size())
-            {
-                cerr << "Field (" << field << ") not found in data file; "
-                     << "Setting it to zero. " << endl;
-                Vmath::Zero(coeffs.num_elements(),coeffs,1);
-                return;
-            }
+            ASSERTL0(i != fielddef->m_fields.size(),
+                     "Field (" + field + ") not found in file.");
 
             // Determine mapping from element ids to location in expansion list
             map<int, int> elmtToExpId;

solvers/CardiacEPSolver/CellModels/CellModel.cpp

@@ -139,7 +139,14 @@ namespace Nektar
             m_gates_tau[i] = Array<OneD, NekDouble>(m_nq);
         }
 
-        v_SetInitialConditions();
+        if (m_session->DefinesFunction("CellModelInitialConditions"))
+        {
+            LoadCellModel();
+        }
+        else
+        {
+            v_SetInitialConditions();
+        }
     }
 
     /**
@@ -304,4 +311,144 @@ namespace Nektar
 
         return outarray;
     }
+
+    void CellModel::LoadCellModel()
+    {
+        const bool root = (m_session->GetComm()->GetRank() == 0);
+        const std::string fncName = "CellModelInitialConditions";
+        std::string varName;
+        Array<OneD, NekDouble> coeffs(m_field->GetNcoeffs());
+        Array<OneD, NekDouble> tmp;
+
+        SpatialDomains::MeshGraphSharedPtr vGraph = m_field->GetGraph();
+
+        if (root)
+        {
+            cout << "Cell model initial conditions: " << endl;
+        }
+
+        // Load each cell model variable
+        // j=0 and j=1 are for transmembrane or intra/extra-cellular volt.
+        Vmath::Zero(m_nq, m_cellSol[0], 1);
+        for(int j = 1; j < m_cellSol.num_elements(); ++j)
+        {
+            // Get the name of the jth variable
+            varName = GetCellVarName(j);
+
+            // Check if this variable is defined in a file or analytically
+            if (m_session->GetFunctionType(fncName, varName) ==
+                    LibUtilities::eFunctionTypeFile)
+            {
+                const std::string file =
+                        m_session->GetFunctionFilename(fncName, varName);
+
+                if (root)
+                {
+                    cout << "  - Field " << varName << ": from file "
+                         << file << endl;
+                }
+
+                std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
+                std::vector<std::vector<NekDouble> > FieldData;
+
+                // Read the restart file containing this variable
+                LibUtilities::Import(file, FieldDef, FieldData);
+
+                LibUtilities::FieldMetaDataMap fieldMetaDataMap;
+                LibUtilities::FieldMetaDataMap::iterator iter;
+                LibUtilities::ImportFieldMetaData(file,fieldMetaDataMap);
+                iter = fieldMetaDataMap.find("Time");
+                if(iter != fieldMetaDataMap.end())
+                {
+                    m_lastTime = iter->second;
+                }
+
+                // Extract the data into the modal coefficients
+                for(int i = 0; i < FieldDef.size(); ++i)
+                {
+                    m_field->ExtractDataToCoeffs(FieldDef[i],
+                                                 FieldData[i],
+                                                 varName,
+                                                 coeffs);
+                }
+
+                // If using nodal cell model then we do a modal->nodal transform
+                // otherwise we do a backward transform onto physical points.
+                if (m_useNodal)
+                {
+                    for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
+                    {
+                        int coef_offset = m_field->GetCoeff_Offset(i);
+                        if (m_field->GetExp(0)->DetShapeType() ==
+                                LibUtilities::eTriangle)
+                        {
+                            m_nodalTri->ModalToNodal(coeffs+coef_offset,
+                                                tmp=m_cellSol[j]+coef_offset);
+                        }
+                        else
+                        {
+                            m_nodalTet->ModalToNodal(coeffs+coef_offset,
+                                                tmp=m_cellSol[j]+coef_offset);
+                        }
+                    }
+                }
+                else
+                {
+                    m_field->BwdTrans(coeffs, m_cellSol[j]);
+                }
+            }
+            else if (m_session->GetFunctionType(fncName, varName) ==
+                    LibUtilities::eFunctionTypeExpression)
+            {
+                LibUtilities::EquationSharedPtr equ =
+                        m_session->GetFunction(fncName, varName);
+
+                if (root)
+                {
+                    cout << "  - Field " << varName << ": "
+                         << equ->GetExpression() << endl;
+                }
+
+                const unsigned int nphys = m_field->GetNpoints();
+                Array<OneD, NekDouble> x0(nphys);
+                Array<OneD, NekDouble> x1(nphys);
+                Array<OneD, NekDouble> x2(nphys);
+                m_field->GetCoords(x0,x1,x2);
+
+                if (m_useNodal)
+                {
+                    Array<OneD, NekDouble> phys(nphys);
+                    equ->Evaluate(x0, x1, x2, phys);
+                    for (unsigned int i = 0; i < m_field->GetNumElmts(); ++i)
+                    {
+                        int phys_offset = m_field->GetPhys_Offset(i);
+                        int coef_offset = m_field->GetCoeff_Offset(i);
+                        if (m_field->GetExp(0)->DetShapeType() ==
+                                LibUtilities::eTriangle)
+                        {
+                            m_field->GetExp(0)->FwdTrans(
+                                            phys + phys_offset,
+                                            m_nodalTri->UpdateCoeffs());
+                            m_nodalTri->ModalToNodal(
+                                            m_nodalTri->GetCoeffs(),
+                                            tmp = m_cellSol[j] + coef_offset);
+                        }
+                        else
+                        {
+                            m_field->GetExp(0)->FwdTrans(
+                                            phys + phys_offset,
+                                            m_nodalTet->UpdateCoeffs());
+                            m_nodalTet->ModalToNodal(
+                                            m_nodalTet->GetCoeffs(),
+                                            tmp = m_cellSol[j] + coef_offset);
+                        }
+                    }
+                }
+                else
+                {
+                    equ->Evaluate(x0, x1, x2, m_cellSol[j]);
+                }
+            }
+        }
+    }
 }

solvers/CardiacEPSolver/CellModels/CellModel.h

@@ -96,6 +96,11 @@ namespace Nektar
             return m_nvar;
         }
 
+        std::string GetCellVarName(unsigned int idx)
+        {
+            return v_GetCellVarName(idx);
+        }
+
         Array<OneD, NekDouble> GetCellSolutionCoeffs(unsigned int idx);
 
     protected:
@@ -139,7 +144,14 @@ namespace Nektar
 
         virtual void v_PrintSummary(std::ostream &out) = 0;
 
+        virtual std::string v_GetCellVarName(unsigned int idx)
+        {
+            return "Var" + boost::lexical_cast<std::string>(idx);
+        }
+
         virtual void v_SetInitialConditions() = 0;
+
+        void LoadCellModel();
     };
 
 }

solvers/CardiacEPSolver/CellModels/CourtemancheRamirezNattel98.cpp

@@ -147,7 +147,6 @@ namespace Nektar
 
         // Variables
         //  0   V    membrane potential
-        //  1   -    unused
         //  2   m    fast sodium current m gate
         //  3   h    fast sodium current h gate
         //  4   j    fast sodium current j gate
@@ -575,4 +574,33 @@ namespace Nektar
         Vmath::Fill(m_nq, 1.488,      m_cellSol[19], 1);
         Vmath::Fill(m_nq, 1.488,      m_cellSol[20], 1);
     }
+
+    std::string CourtemancheRamirezNattel98::v_GetCellVarName(unsigned int idx)
+    {
+        switch (idx)
+        {
+            case 0:  return "u";
+            case 1:  return "m";
+            case 2:  return "h";
+            case 3:  return "j";
+            case 4:  return "o_a";
+            case 5:  return "o_i";
+            case 6:  return "u_a";
+            case 7:  return "u_i";
+            case 8:  return "x_r";
+            case 9: return "x_s";
+            case 10: return "d";
+            case 11: return "f";
+            case 12: return "f_Ca";
+            case 13: return "U";
+            case 14: return "V";
+            case 15: return "W";
+            case 16: return "Na_i";
+            case 17: return "Ca_i";
+            case 18: return "K_i";
+            case 19: return "Ca_rel";
+            case 20: return "Ca_up";
+        }
+    }
+
 }

solvers/CardiacEPSolver/CellModels/CourtemancheRamirezNattel98.h

@@ -76,6 +76,8 @@ namespace Nektar
 
         virtual void v_SetInitialConditions();
 
+        virtual std::string v_GetCellVarName(unsigned int idx);
+
     private:
         NekDouble C_m;
         NekDouble g_Na;

solvers/CardiacEPSolver/Filters/FilterCheckpointCellModel.cpp

@@ -73,12 +73,13 @@ namespace Nektar
 
         m_index = 0;
         m_outputIndex = 0;
+
+        v_Update(pFields, 0.0);
     }
 
     void FilterCheckpointCellModel::v_Update(const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields, const NekDouble &time)
     {
-        m_index++;
-        if (m_index % m_outputFrequency > 0)
+        if (m_index++ % m_outputFrequency > 0)
         {
             return;
         }
@@ -99,19 +100,27 @@ namespace Nektar
         std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
 
         // copy Data into FieldData and set variable
-        for(int j = 0; j < m_cell->GetNumCellVariables(); ++j)
+        std::string varName;
+        for(int j = 1; j < m_cell->GetNumCellVariables(); ++j)
         {
+            varName = m_cell->GetCellVarName(j);
+
             for(int i = 0; i < FieldDef.size(); ++i)
             {
                 // Retrieve data from cell model
                 Array<OneD, NekDouble> data = m_cell->GetCellSolutionCoeffs(j);
 
                 // Could do a search here to find correct variable
-                FieldDef[i]->m_fields.push_back(boost::lexical_cast<std::string>(j));
+                FieldDef[i]->m_fields.push_back(varName);
                 pFields[0]->AppendFieldData(FieldDef[i], FieldData[i], data);
             }
         }
-        LibUtilities::Write(vOutputFilename.str(),FieldDef,FieldData);
+
+        // Update time in field info if required
+        LibUtilities::FieldMetaDataMap fieldMetaDataMap;
+        fieldMetaDataMap["Time"] =  time;
+
+        LibUtilities::Write(vOutputFilename.str(),FieldDef,FieldData,fieldMetaDataMap);
         m_outputIndex++;
     }
 

utilities/PostProcessing/CombinePartitions.cpp

@@ -37,7 +37,7 @@ int main(int argc, char *argv[])
         TiXmlElement *nektar = docInput.FirstChildElement("NEKTAR");
         
         // load up root processor's meta data
-        if(n == 0)
+        if(n == 0 && nektar->FirstChildElement("FIELDMETADATA"))
         {
             TiXmlElement *metadata = nektar->FirstChildElement("FIELDMETADATA");
             master->LinkEndChild(new TiXmlElement(*metadata));