From e2e0785e61b040906d3d75cd845a3d4d6e2e04ee Mon Sep 17 00:00:00 2001
From: David Moxey <d.moxey@imperial.ac.uk>
Date: Tue, 12 Jun 2012 13:01:35 +0000
Subject: [PATCH] Fixing reading of initial conditions and L2 error for
 homogeneous fields. More verbose output for initial conditions in
 EquationSystem.

git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3918 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
---
 library/MultiRegions/ExpList.cpp              | 45 +++++++-------
 .../MultiRegions/ExpList3DHomogeneous1D.cpp   | 20 +++---
 library/MultiRegions/ExpListHomogeneous1D.cpp | 62 ++++++++++++-------
 library/MultiRegions/ExpListHomogeneous1D.h   |  2 +-
 library/SolverUtils/EquationSystem.cpp        | 13 +++-
 5 files changed, 85 insertions(+), 57 deletions(-)

diff --git a/library/MultiRegions/ExpList.cpp b/library/MultiRegions/ExpList.cpp
index 53d30240b1..2c96a48889 100644
--- a/library/MultiRegions/ExpList.cpp
+++ b/library/MultiRegions/ExpList.cpp
@@ -1758,7 +1758,7 @@ namespace Nektar
                 (*m_exp)[i]->SetPhys(m_phys+m_phys_offset[i]);
                 err  = std::max(err,(*m_exp)[i]->Linf());
             }
-            m_comm->AllReduce(err, LibUtilities::ReduceMax);
+            m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
             return err;
         }
 
@@ -1826,35 +1826,36 @@ namespace Nektar
                 errl2 = (*m_exp)[i]->L2();
                 err += errl2*errl2;
             }
-            m_comm->AllReduce(err, LibUtilities::ReduceSum);
+            m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
+            
             return sqrt(err);
         }
 		
-		Array<OneD, NekDouble> ExpList::v_HomogeneousEnergy (void)
+        Array<OneD, NekDouble> ExpList::v_HomogeneousEnergy (void)
         {
-			ASSERTL0(false,
+            ASSERTL0(false,
                      "This method is not defined or valid for this class type");
-			Array<OneD, NekDouble> NoEnergy(1,0.0);
+            Array<OneD, NekDouble> NoEnergy(1,0.0);
             return NoEnergy;
         }
 		
-		Array<OneD, unsigned int> ExpList::v_GetZIDs(void)
-		{
-			ASSERTL0(false,
+        Array<OneD, unsigned int> ExpList::v_GetZIDs(void)
+        {
+            ASSERTL0(false,
                      "This method is not defined or valid for this class type");
-			Array<OneD, unsigned int> NoModes(1);
+            Array<OneD, unsigned int> NoModes(1);
 			
-			return NoModes;
-		}
+            return NoModes;
+        }
 		
-		Array<OneD, unsigned int> ExpList::v_GetYIDs(void)
-		{
-			ASSERTL0(false,
+        Array<OneD, unsigned int> ExpList::v_GetYIDs(void)
+        {
+            ASSERTL0(false,
                      "This method is not defined or valid for this class type");
-			Array<OneD, unsigned int> NoModes(1);
+            Array<OneD, unsigned int> NoModes(1);
 			
-			return NoModes;
-		}
+            return NoModes;
+        }
 		
 
         /**
@@ -2050,11 +2051,11 @@ namespace Nektar
             v_ExtractDataToCoeffs(fielddef,fielddata,field,m_coeffs);
         }
 		
-		//3D-Base Flow (implementation in the homogeneous classes)
-		void ExpList::v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field, bool BaseFlow3D)
-		{
-			ASSERTL0(false, "This method is not defined or valid for this class type");
-		}
+        //3D-Base Flow (implementation in the homogeneous classes)
+        void ExpList::v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field, bool BaseFlow3D)
+        {
+            ASSERTL0(false, "This method is not defined or valid for this class type");
+        }
         
         void ExpList::v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field, Array<OneD, NekDouble> &coeffs)
         {     	
diff --git a/library/MultiRegions/ExpList3DHomogeneous1D.cpp b/library/MultiRegions/ExpList3DHomogeneous1D.cpp
index 5b68c750ef..3c6bd7e3bf 100644
--- a/library/MultiRegions/ExpList3DHomogeneous1D.cpp
+++ b/library/MultiRegions/ExpList3DHomogeneous1D.cpp
@@ -395,25 +395,25 @@ namespace Nektar
             for(int n = 0; n < m_planes.num_elements(); ++n)
             {
                 errL2 = m_planes[n]->L2(soln + cnt);
-                cnt += m_planes[n]->GetTotPoints();
-                err += errL2*errL2*local_w[n]*m_lhom*0.5;
+                cnt  += m_planes[n]->GetTotPoints();
+                err  += errL2*errL2*local_w[n]*m_lhom*0.5;
             }
             
             m_comm->GetColumnComm()->AllReduce(err, LibUtilities::ReduceSum);
-
+            
             return sqrt(err);
         }
-		
+	
         NekDouble ExpList3DHomogeneous1D::v_L2(void)
         {
             NekDouble errL2,err = 0;
             Array<OneD, const NekDouble> w = m_homogeneousBasis->GetW();
-			Array<OneD, NekDouble> local_w(m_planes.num_elements());
-			
-			for(int n = 0; n < m_planes.num_elements(); n++)
-			{
-				local_w[n] = w[m_transposition->GetPlaneID(n)];
-			}
+            Array<OneD, NekDouble> local_w(m_planes.num_elements());
+            
+            for(int n = 0; n < m_planes.num_elements(); n++)
+            {
+                local_w[n] = w[m_transposition->GetPlaneID(n)];
+            }
 
             for(int n = 0; n < m_planes.num_elements(); ++n)
             {
diff --git a/library/MultiRegions/ExpListHomogeneous1D.cpp b/library/MultiRegions/ExpListHomogeneous1D.cpp
index 628c568f97..ef43e94bc0 100644
--- a/library/MultiRegions/ExpListHomogeneous1D.cpp
+++ b/library/MultiRegions/ExpListHomogeneous1D.cpp
@@ -634,7 +634,7 @@ namespace Nektar
         {
             int i,n;
             int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
-
+            
             // Determine mapping from element ids to location in
             // expansion list
             map<int, int> ElmtID_to_ExpID;
@@ -655,21 +655,31 @@ namespace Nektar
             }
         }
 		
-		void ExpListHomogeneous1D::v_AppendFieldData(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
+        void ExpListHomogeneous1D::v_AppendFieldData(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
         {
            v_AppendFieldData(fielddef,fielddata,m_coeffs);
         }
 
         //Extract the data in fielddata into the m_coeff list
-        void ExpListHomogeneous1D::v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field)
+        void ExpListHomogeneous1D::v_ExtractDataToCoeffs(
+            SpatialDomains::FieldDefinitionsSharedPtr &fielddef,
+            std::vector<NekDouble>                    &fielddata,
+            std::string                               &field,
+            Array<OneD, NekDouble>                    &coeffs)
         {
             int i,n;
             int offset = 0;
-            int nzmodes; 
+            int nzmodes;
             int datalen = fielddata.size()/fielddef->m_fields.size();
             int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
             
-
+            // Build map of plane IDs lying on this process.
+            std::map<int,int> homoZids;
+            for (i = 0; i < m_planes.num_elements(); ++i)
+            {
+                homoZids[m_transposition->GetPlaneID(i)] = i;
+            }
+            
             for(i = 0; i < fielddef->m_basis.size(); ++i)
             {
                 if(fielddef->m_basis[i] == m_homogeneousBasis->GetBasisType())
@@ -678,7 +688,7 @@ namespace Nektar
                     break;
                 }
             }
-            ASSERTL1(i != fielddef->m_basis.size()," Failed to determine number of modes");
+            ASSERTL1(i != fielddef->m_basis.size(),"Failed to determine number of modes");
             
             // Find data location according to field definition
             for(i = 0; i < fielddef->m_fields.size(); ++i)
@@ -689,11 +699,10 @@ namespace Nektar
                 }
                 offset += datalen;
             }
+            ASSERTL0(i != fielddef->m_fields.size(),
+                     "Field " + field + " not found in data file");
             
-            ASSERTL0(i!= fielddef->m_fields.size(),"Field not found in data file");
-            
-            // Determine mapping from element ids to location in
-            // expansion list
+            // Determine mapping from element ids to location in expansion list.
             map<int, int> ElmtID_to_ExpID;
             for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
             {
@@ -701,36 +710,44 @@ namespace Nektar
             }
 
             int modes_offset = 0;
-			int planes_offset = 0;
-            Array<OneD, NekDouble> coeff_tmp;             
+            int planes_offset = 0;
+            Array<OneD, NekDouble> coeff_tmp;
+            
             for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
             {
                 int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
-                int datalen = (*m_exp)[eid]->CalcNumberOfCoefficients(fielddef->m_numModes,modes_offset);
+                int datalen = (*m_exp)[eid]->CalcNumberOfCoefficients(
+                    fielddef->m_numModes,modes_offset);
+                
                 if(fielddef->m_uniOrder == true) // reset modes_offset to zero
                 {
                     modes_offset = 0;
                 }
-
-                for(n = 0; n < nzmodes; ++n)
+                
+                for(n = 0; n < nzmodes; ++n, offset += datalen)
                 {
-					planes_offset = fielddef->m_homogeneousZIDs[n];
+                    std::map<int,int>::iterator it = homoZids.find(
+                        fielddef->m_homogeneousZIDs[n]);
+                    
+                    // Check to make sure this mode number lies in this field.
+                    if (it == homoZids.end())
+                    {
+                        continue;
+                    }
+                    
+                    planes_offset = it->second;
                     if(datalen == (*m_exp)[eid]->GetNcoeffs())
                     {
-                        Vmath::Vcopy(datalen,&fielddata[offset],1,&m_coeffs[m_coeff_offset[eid]+ planes_offset*ncoeffs_per_plane],1);
+                        Vmath::Vcopy(datalen,&fielddata[offset],1,&coeffs[m_coeff_offset[eid]+planes_offset*ncoeffs_per_plane],1);
                     }
                     else // unpack data to new order
                     {
-                        (*m_exp)[eid]->ExtractDataToCoeffs(fielddata, offset, fielddef->m_numModes,modes_offset,coeff_tmp = m_coeffs + m_coeff_offset[eid] + planes_offset*ncoeffs_per_plane);
+                        (*m_exp)[eid]->ExtractDataToCoeffs(fielddata, offset, fielddef->m_numModes,modes_offset,coeff_tmp = coeffs + m_coeff_offset[eid] + planes_offset*ncoeffs_per_plane);
                     }
-                    
-                    offset += datalen;
                 }
             }
         }
 		
-		
-		
         //Extract the data in fielddata into the m_coeff list (for 2D files into 3D cases)
         void ExpListHomogeneous1D::v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field, bool BaseFlow3D)
         {
@@ -739,7 +756,6 @@ namespace Nektar
             int nzmodes = m_homogeneousBasis->GetNumModes();
             int datalen = fielddata.size()/fielddef->m_fields.size();
             int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
-		
 			
             // Find data location according to field definition
             for(i = 0; i < fielddef->m_fields.size(); ++i)
diff --git a/library/MultiRegions/ExpListHomogeneous1D.h b/library/MultiRegions/ExpListHomogeneous1D.h
index 2e290012d6..6447b44195 100644
--- a/library/MultiRegions/ExpListHomogeneous1D.h
+++ b/library/MultiRegions/ExpListHomogeneous1D.h
@@ -184,7 +184,7 @@ namespace Nektar
 			
 			virtual void v_AppendFieldData(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs);
 
-            virtual void v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field);
+            virtual void v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field, Array<OneD, NekDouble> &coeffs);
 			
 			virtual void v_ExtractDataToCoeffs(SpatialDomains::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field,bool BaseFlow3D);
 
diff --git a/library/SolverUtils/EquationSystem.cpp b/library/SolverUtils/EquationSystem.cpp
index a6f6fa4a47..bdf6a4e7e9 100644
--- a/library/SolverUtils/EquationSystem.cpp
+++ b/library/SolverUtils/EquationSystem.cpp
@@ -778,13 +778,24 @@ namespace Nektar
                 LibUtilities::EquationSharedPtr ffunc
                     = m_session->GetFunction(pFunctionName, pFieldName);
 
+                if (m_comm->GetRank() == 0)
+                {
+                    cout << "- Field " << pFieldName << ": " 
+                         << ffunc->GetExpression() << endl;
+                }
+
                 ffunc->Evaluate(x0,x1,x2,pTime,pArray);
             }
             else if (vType == LibUtilities::eFunctionTypeFile)
             {
                 std::string filename
                     = m_session->GetFunctionFilename(pFunctionName, pFieldName);
-                cout << pFunctionName << " from file: " << filename << endl;
+                
+                if (m_comm->GetRank() == 0)
+                {
+                    cout << "- Field " << pFieldName << ": " 
+                         << "from file " << filename << endl;
+                }
 
                 std::vector<SpatialDomains::FieldDefinitionsSharedPtr> FieldDef;
                 std::vector<std::vector<NekDouble> > FieldData;
-- 
GitLab