Commit 0fb7dee7 authored by Dmitry Bagaev's avatar Dmitry Bagaev
Browse files

Merge remote-tracking branch 'INMOST-DEV/master'

# Conflicts:
#	Examples/ADFVDiscr/main.cpp
#	Examples/ADMFD/main.cpp
#	Examples/DrawGrid/main.cpp
#	Examples/GridGen/main.cpp
#	Examples/OldDrawGrid/main.cpp
#	Source/Misc/xml.cpp
#	Source/Solver/solver_petsc.cpp
#	Source/Solver/sparse.cpp
parents 96f0d91c eef2ef04
......@@ -64,57 +64,6 @@ The example of the layout of the file is presented below.
<ParallelMesh Number="2" [Layers="2" Element="Face"]>
<!-- Definition of the first mesh-->
<Mesh [Name="Box"] [RepairOrientation="False"]>
<!-- Describe data that is defined on the mesh -->
<!-- Number can be used to match
the number of entries and
optimize memory allocation -->
<Tags [Number="8"]>
<!-- See further text for description of each attribute -->
<Tag Name="GLOBAL_ID"
Size="1"
Type="Integer"
Sparse="Sets"
Definition="Cells,Faces,Nodes,Sets" />
<!-- Size is optional, default: Variable,
Type is optional, default: Real,
Sparsity is optional, default: None,
Name and Definition are mandatory -->
<Tag Name="PERMEABILITY_TENSOR"
[Size="6"]
[Type="Real" ]
[Sparse="None"]
Definition="Cells" />
<Tag Name="BOUNDARY_PRESSURE"
Size="1"
Type="Real"
Sparse="Faces"
Definition="Faces" />
<Tag Name="WELL_INDEX_WELL0"
Size="1"
Type="Real"
Sparse="Cells"
Definition="Cells" />
<Tag Name="BOUNDARY_DISPLACEMENT"
Size="3"
Type="Real"
Sparse="True"
Definition="Nodes" />
<Tag Name="CONNECTIONS"
Size="Variable"
Type="Reference"
Sparse="None"
Definition="Faces" />
<Tag Name="FLUX"
Size="1"
Type="Variable"
Sparse="None"
Definition="Faces" />
<Tag Name="WHATEVER"
Size="Variable"
Type="Bulk"
Sparse="Sets"
Definition="Sets" />
</Tags>
<!-- Define all the nodes of the mesh -->
<Nodes Number="100" [Dimension="3"]>
<![CDATA[
......@@ -202,6 +151,57 @@ The example of the layout of the file is presented below.
</Set>
<Set> ... </Set>
</Sets>
<!-- Describe data that is defined on the mesh -->
<!-- Number can be used to match
the number of entries and
optimize memory allocation -->
<Tags [Number="8"]>
<!-- See further text for description of each attribute -->
<Tag Name="GLOBAL_ID"
Size="1"
Type="Integer"
Sparse="Sets"
Definition="Cells,Faces,Nodes,Sets" />
<!-- Size is optional, default: Variable,
Type is optional, default: Real,
Sparsity is optional, default: None,
Name and Definition are mandatory -->
<Tag Name="PERMEABILITY_TENSOR"
[Size="6"]
[Type="Real" ]
[Sparse="None"]
Definition="Cells" />
<Tag Name="BOUNDARY_PRESSURE"
Size="1"
Type="Real"
Sparse="Faces"
Definition="Faces" />
<Tag Name="WELL_INDEX_WELL0"
Size="1"
Type="Real"
Sparse="Cells"
Definition="Cells" />
<Tag Name="BOUNDARY_DISPLACEMENT"
Size="3"
Type="Real"
Sparse="True"
Definition="Nodes" />
<Tag Name="CONNECTIONS"
Size="Variable"
Type="Reference"
Sparse="None"
Definition="Faces" />
<Tag Name="FLUX"
Size="1"
Type="Variable"
Sparse="None"
Definition="Faces" />
<Tag Name="WHATEVER"
Size="Variable"
Type="Bulk"
Sparse="Sets"
Definition="Sets" />
</Tags>
<!-- Define the data, the number of data sets should be equal to the number of tags -->
<Data Number="11">
<!-- TagName maps data set to the tag -->
......@@ -310,25 +310,20 @@ The example of the layout of the file is presented below.
The purpose of this XML tag is to optionally describe the number of meshes contained in the file in attribute "Number". As an additional information it can provide number of ghost layers between meshes in attribute "Layers" and the type of elements used to compute adjacency for ghost layers in attribute "Element". This information is reserved for the future use.
\section{Mesh}
This xml tag wraps all the data of the mesh. The optional attribute "Name" can be used to address elements of this mesh from another mesh. The optional attribute "RepairOrientation" if set to "True" will correct the orientation of the faces.
\section{Tags}
This xml tag declares the tags of the data that are present on the mesh. Optional attribute "Number" corresponds to the total number of tags defined.
\subsection{Tag}
Each xml tag contains attributes that describe the way the data is stored on the mesh. The following attributes can be defined:
\begin{itemize}
\item[Name] The name assigned to the data. Cannot be the same for two different tags. Mandatory to provide.
\item[Size] Number of records of the data on each element. Can be any positive number or "Variable". "Variable" means that the data may have different number of entries on each element. Optional, default: "Variable".
\item[Type] Type of the data that tag represents, can be either "Real" or "Integer" or "Bulk", or "Variable" or "Reference" or "RemoteReference". "Real" corresponds to floating numbers; "Integer" - to integral numbers; "Bulk" can represent any binary data; "Variable" can store a floating point value and corresponding variations represented by pairs of integral number and floating point coefficient; "Reference" and "RemoteReference" represent links to mesh elements or sets or mesh itself. Optional, default: "Real".
\item[Sparse] Defines that data is present only on some elements of the indicated element type. Can be either Sets or Cells or Faces or Edges or Nodes or None. Optional, default: None.
\item[Definition] Defines types of elements that posses the data. Can be Mesh, Sets, Cells, Faces, Edges or Nodes. Mandatory to provide.
\end{itemize}
The "GLOBAL\_ID" tag in the example defines a single integer entry on all the nodes, faces and cells, and on some sets. The "WHATEVER" tag defines a binary data of arbitrary size only to some sets. When one would like to have just a handful of entries unrelated to the mesh elements one can prescribe them to the mesh, i.e. write Definition="Mesh".
\section{Nodes}
This is a mandatory XML tag for the file that describes all the nodes of the mesh. The "Number" attribute describes the total number of nodes. Optional attribute "Dimension" tells the number of space dimensions, that is the number of coordinates in each entry, defaults to 3. The contents of the XML tag inside of "$<![CDATA[]]>$" can be entered in any format suitable to represent vectors as described in \S~\ref{data_format}. Nodes are mesh elements of dimension 0.
\section{Faces, Edges, Cells}
XML Tags Faces and Edges are optional and could be used to define only some of the elements of the mesh. For example one may introduce only boundary faces in XML tag Faces. The XML Tag Cells is mandatory and it represents cells of the mesh. One can optionally provide an optional attribute "Number" that will be used to check that the number of elements red is correct. Edges, faces and cells are mesh elements of dimension 1, 2 and 3 respectively.
\subsubsection{Connections}
XML tag "Connections" is used to provide the connection of each listed element to elements of lower dimension. Attributes are "Type" that describe types of listed elements; "Number" - total number of listed elements; "Offset" - optional attribute that describes first position of listed elements, default 0. It's required that the type of elements in "Type" has lower dimension than constructed element. It is possible to describe some cells constructed of nodes and some cells constructed of faces by using consecutive XML tags "Connections", see example for details.
XML tag "Connections" is used to provide the connection of each listed element to elements of lower dimension. Attributes are:
\begin{itemize}
\item[Type] Describe types of listed elements. It's required that the type of elements in "Type" has lower dimension than constructed element.
\item[Number] Total number of listed elements, mandatory.
\item[Offset] Optional attribute that describes first position of listed elements, default 0.
\item[Dimensions] This attribute is used to distinguish definition of 3D cells from 2D cells defined by nodes. Set to 2 to define 2D polygons and to 3 to define 3D volumetric cells, default 3.
\end{itemize}
It is possible to describe some cells constructed of nodes and some cells constructed of faces by using consecutive XML tags "Connections", see example for details.
The content of the XML tag inside of "$<![CDATA[]]>$" should start from the number of elements connected followed by a list of positions of elements. For example if we have a record $3 1 2 3$ for a face with connection to nodes then this record specifies that the face consists of 3 nodes namely node 1, node 2 and node 3. Enumeration of nodes here corresponds to the order of nodes in which their coordinates are listed in XML tag "Nodes".
\section{Sets}
......@@ -345,6 +340,19 @@ Describes each set. Sets could be arranged into an arbitrary tree, such as an oc
\item [Offset] Sets a shift in the enumeration of elements. Optional, default: 0.
\end{itemize}
In the contents of the XML tag "Set" inside of "$<![CDATA[]]>$" all the mesh elements are listed in whatever order is needed. They will be reordered internally if ordering was prescribed. The ordering happens after the mesh data is attached to the elements. Each element is listed as "type:position" where type can be either Mesh or Set or Cell or Face or Edge or Node. The position corresponds to the position in the order of elements in which their records are encountered in corresponding XML tags.
\section{Tags}
This xml tag declares the tags of the data that are present on the mesh. Optional attribute "Number" corresponds to the total number of tags defined.
\subsection{Tag}
Each xml tag contains attributes that describe the way the data is stored on the mesh. The following attributes can be defined:
\begin{itemize}
\item[Name] The name assigned to the data. Cannot be the same for two different tags. Mandatory to provide.
\item[Size] Number of records of the data on each element. Can be any positive number or "Variable". "Variable" means that the data may have different number of entries on each element. Optional, default: "Variable".
\item[Type] Type of the data that tag represents, can be either "Real" or "Integer" or "Bulk", or "Variable" or "Reference" or "RemoteReference". "Real" corresponds to floating numbers; "Integer" - to integral numbers; "Bulk" can represent any binary data; "Variable" can store a floating point value and corresponding variations represented by pairs of integral number and floating point coefficient; "Reference" and "RemoteReference" represent links to mesh elements or sets or mesh itself. Optional, default: "Real".
\item[Sparse] Defines that data is present only on some elements of the indicated element type. Can be either Sets or Cells or Faces or Edges or Nodes or None. Optional, default: None.
\item[Definition] Defines types of elements that posses the data. Can be Mesh, Sets, Cells, Faces, Edges or Nodes. Mandatory to provide.
\end{itemize}
The "GLOBAL\_ID" tag in the example defines a single integer entry on all the nodes, faces and cells, and on some sets. The "WHATEVER" tag defines a binary data of arbitrary size only to some sets. When one would like to have just a handful of entries unrelated to the mesh elements one can prescribe them to the mesh, i.e. write Definition="Mesh".
\section{Data}
Encloses all the data of the mesh. Optional attribute "Number" can be added to control the number of "DataSet" XML tags encountered.
\subsection{DataSet} \label{tag_dataset}
......
......@@ -31,17 +31,24 @@ Storage::real dot_prod(Storage::real v1[3], Storage::real v2[3])
Storage::real func(Storage::real x[3], Storage::real tmp)
{
<<<<<<< HEAD
// return x[0] + 2 * x[1] + 3 * x[2];
double s0 = sin (M_PI * x[0]);
double s1 = sin (M_PI * x[1]);
double s2 = sin (M_PI * x[2]);
=======
// return x[0] + 2 * x[1] + 3 * x[2];
double s0 = sin (M_PI * x[0]);
double s1 = sin (M_PI * x[1]);
double s2 = sin (M_PI * x[2]);
>>>>>>> INMOST-DEV/master
return s0 * s1 * s2;
(void) tmp;
}
Storage::real func_rhs(Storage::real x[3], Storage::real tmp)
{
// return 0;
// return 0;
return -3 * tmp * M_PI * M_PI * sin (M_PI * x[0]) * sin (M_PI * x[1]) * sin (M_PI * x[2]);
}
......@@ -73,7 +80,12 @@ int main(int argc,char ** argv)
if( m->GetProcessorRank() == 0 )
m->Load(argv[1]); // Load mesh from the serial file format
}
<<<<<<< HEAD
BARRIER
=======
BARRIER;
>>>>>>> INMOST-DEV/master
if( m->GetProcessorRank() == 0 ) std::cout << "Processors: " << m->GetProcessorsNumber() << std::endl;
if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_File): " << Timer()-ttt << std::endl;
......@@ -87,21 +99,26 @@ int main(int argc,char ** argv)
//~ if( m->GetProcessorRank() == 0 ) std::cout << "Load(MPI_Scatter): " << Timer()-ttt2 << std::endl;
#if defined(USE_PARTITIONER)
<<<<<<< HEAD
if (!repartition)
{ // currently only non-distributed meshes are supported by Inner_RCM partitioner
=======
if (!repartition)
{ // currently only non-distributed meshes are supported by Inner_RCM partitioner
>>>>>>> INMOST-DEV/master
ttt = Timer();
Partitioner * p = new Partitioner(m);
p->SetMethod(Partitioner::Inner_RCM,Partitioner::Partition); // Specify the partitioner
p->Evaluate(); // Compute the partitioner and store new processor ID in the mesh
delete p;
BARRIER
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Evaluate: " << Timer()-ttt << std::endl;
ttt = Timer();
m->Redistribute(); // Redistribute the mesh data
m->ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
BARRIER
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Redistribute: " << Timer()-ttt << std::endl;
}
......@@ -109,15 +126,24 @@ int main(int argc,char ** argv)
ttt = Timer();
m->AssignGlobalID(CELL | EDGE | FACE | NODE);
<<<<<<< HEAD
BARRIER
if( m->GetProcessorRank() == 0 ) std::cout << "Assign id: " << Timer()-ttt << std::endl;
=======
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Assign id: " << Timer()-ttt << std::endl;
>>>>>>> INMOST-DEV/master
id = m->GlobalIDTag(); // Get the tag of the global ID
//m->Save("solution_check_0.vtk");
phi = m->CreateTag("Solution",DATA_REAL,CELL,NONE,1); // Create a new tag for the solution phi
tensor_K = m->CreateTag("K",DATA_REAL,CELL,NONE,1); // Create a new tag for K tensor
//m->Save("solution_check_1.vtk");
//m->Save("solution_check_1.vtk");
<<<<<<< HEAD
=======
>>>>>>> INMOST-DEV/master
for( Mesh::iteratorCell cell = m->BeginCell(); cell != m->EndCell(); ++cell ) // Loop over mesh cells
if( cell->GetStatus() != Element::Ghost ) // If the cell is an own one
cell->Real(tensor_K) = 1.0; // Store the tensor K value into the tag
......@@ -125,16 +151,24 @@ int main(int argc,char ** argv)
ttt = Timer();
m->ExchangeGhost(1,FACE);
m->ExchangeData(tensor_K,CELL,0); // Exchange the tensor_K data over processors
BARRIER
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Exchange ghost: " << Timer()-ttt << std::endl;
ttt = Timer();
<<<<<<< HEAD
Solver S("inner_ilu2"); // Specify the linear solver to ASM+ILU2+BiCGStab one
S.SetParameter("absolute_tolerance", "1e-8");
Residual R; // Residual vector
Sparse::LockService Locks;
=======
Solver S(Solver::INNER_ILU2); // Specify the linear solver to ASM+ILU2+BiCGStab one
S.SetParameterReal("absolute_tolerance",1e-8);
S.SetParameterEnum("schwartz_overlap",2);
Residual R; // Residual vector
Sparse::LockService Locks;
>>>>>>> INMOST-DEV/master
Sparse::Vector Update; // Declare the solution and the right-hand side vectors
Mesh::GeomParam table;
......@@ -148,6 +182,7 @@ int main(int argc,char ** argv)
//~ BARRIER
//~ if( m->GetProcessorRank() == 0 ) std::cout << "Prepare geometric data: " << Timer()-ttt << std::endl;
<<<<<<< HEAD
{
Automatizator aut(m);
Automatizator::MakeCurrent(&aut);
......@@ -162,71 +197,88 @@ int main(int argc,char ** argv)
dynamic_variable Phi(aut,iphi);
// Solve \nabla \cdot \nabla phi = f equation
//for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
=======
{
Automatizator aut(m);
Automatizator::MakeCurrent(&aut);
INMOST_DATA_ENUM_TYPE iphi = aut.RegisterDynamicTag(phi,CELL);
aut.EnumerateDynamicTags();
// Set the indeces intervals for the matrix and vectors
R.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
R.InitLocks();
Update.SetInterval(aut.GetFirstIndex(),aut.GetLastIndex());
//~ std::cout << m->GetProcessorRank() << " A,x,b interval " << idmin << ":" << idmax << " size " << idmax-idmin << std::endl;
dynamic_variable Phi(aut,iphi);
// Solve \nabla \cdot \nabla phi = f equation
//for( Mesh::iteratorFace face = m->BeginFace(); face != m->EndFace(); ++face )
>>>>>>> INMOST-DEV/master
#if defined(USE_OMP)
#pragma omp parallel
#endif
{
variable flux; //should be more efficient to define here to avoid multiple memory allocations if storage for variations should be expanded
{
variable flux; //should be more efficient to define here to avoid multiple memory allocations if storage for variations should be expanded
#if defined(USE_OMP)
#pragma omp for
#endif
for(Storage::integer iface = 0; iface < m->FaceLastLocalID(); ++iface ) if( m->isValidFace(iface) )
{
Face face = Face(m,ComposeFaceHandle(iface));
Element::Status s1,s2;
Cell r1 = face->BackCell();
Cell r2 = face->FrontCell();
if( ((!r1->isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1) +
((!r2->isValid() || (s2 = r2->GetStatus()) == Element::Ghost)?0:1) == 0) continue;
Storage::integer i1 = aut.GetDynamicIndex(r1,iphi), i2;
Storage::real f_nrm[3], r1_cnt[3], r2_cnt[3], f_cnt[3], d1, d2, D, v[3], T;
Storage::real f_area = face->Area(); // Get the face area
face->UnitNormal(f_nrm); // Get the face normal
r1->Centroid(r1_cnt); // Get the barycenter of the cell
face->Centroid(f_cnt); // Get the barycenter of the face
if( !r2->isValid() ) // boundary condition
{
Storage::real bnd_pnt[3], dist;
make_vec(f_cnt,r1_cnt,v);
dist = dot_prod(f_nrm,v);
// bnd_pnt is a projection of the cell center to the face
bnd_pnt[0] = r1_cnt[0] + dist * f_nrm[0];
bnd_pnt[1] = r1_cnt[1] + dist * f_nrm[1];
bnd_pnt[2] = r1_cnt[2] + dist * f_nrm[2];
T = r1->Real(tensor_K) * f_area / dist;
//flux = T * (func(bnd_pnt,0) - variable(aut,r1,iphi));
Locks.Lock(i1);
R[i1] -= T * (func(bnd_pnt,0) - Phi(r1));
Locks.UnLock(i1);
}
else
{
i2 = aut.GetDynamicIndex(r2,iphi);
r2->Centroid(r2_cnt);
D = dot_prod(f_nrm,f_cnt);
d1 = fabs(dot_prod(r1_cnt,f_nrm) - D);
d2 = fabs(dot_prod(r2_cnt,f_nrm) - D);
T = 1.0 / (d1/r1->Real(tensor_K) + d2/r2->Real(tensor_K)) * f_area;
//flux = T * (variable(aut,r2,iphi) - variable(aut,r1,iphi));//(unknown(aut,r2,iphi) - unknown(aut,r1,iphi));
flux = T * (Phi(r2) - Phi(r1));
if( s1 != Element::Ghost )
{
Locks.Lock(i1);
R[i1] -= flux;
Locks.UnLock(i1);
}
if( s2 != Element::Ghost )
{
Locks.Lock(i2);
R[i2] += flux;
Locks.UnLock(i2);
}
}
}
}
for(Storage::integer iface = 0; iface < m->FaceLastLocalID(); ++iface ) if( m->isValidFace(iface) )
{
Face face = Face(m,ComposeFaceHandle(iface));
Element::Status s1,s2;
Cell r1 = face->BackCell();
Cell r2 = face->FrontCell();
if( ((!r1->isValid() || (s1 = r1->GetStatus()) == Element::Ghost)?0:1) +
((!r2->isValid() || (s2 = r2->GetStatus()) == Element::Ghost)?0:1) == 0) continue;
Storage::integer i1 = aut.GetDynamicIndex(r1,iphi), i2;
Storage::real f_nrm[3], r1_cnt[3], r2_cnt[3], f_cnt[3], d1, d2, D, v[3], T;
Storage::real f_area = face->Area(); // Get the face area
face->UnitNormal(f_nrm); // Get the face normal
r1->Centroid(r1_cnt); // Get the barycenter of the cell
face->Centroid(f_cnt); // Get the barycenter of the face
if( !r2->isValid() ) // boundary condition
{
Storage::real bnd_pnt[3], dist;
make_vec(f_cnt,r1_cnt,v);
dist = dot_prod(f_nrm,v);
// bnd_pnt is a projection of the cell center to the face
bnd_pnt[0] = r1_cnt[0] + dist * f_nrm[0];
bnd_pnt[1] = r1_cnt[1] + dist * f_nrm[1];
bnd_pnt[2] = r1_cnt[2] + dist * f_nrm[2];
T = r1->Real(tensor_K) * f_area / dist;
//flux = T * (func(bnd_pnt,0) - variable(aut,r1,iphi));
R.Lock(i1);
R[i1] -= T * (func(bnd_pnt,0) - Phi(r1));
R.UnLock(i1);
}
else
{
i2 = aut.GetDynamicIndex(r2,iphi);
r2->Centroid(r2_cnt);
D = dot_prod(f_nrm,f_cnt);
d1 = fabs(dot_prod(r1_cnt,f_nrm) - D);
d2 = fabs(dot_prod(r2_cnt,f_nrm) - D);
T = 1.0 / (d1/r1->Real(tensor_K) + d2/r2->Real(tensor_K)) * f_area;
//flux = T * (variable(aut,r2,iphi) - variable(aut,r1,iphi));//(unknown(aut,r2,iphi) - unknown(aut,r1,iphi));
flux = T * (Phi(r2) - Phi(r1));
if( s1 != Element::Ghost )
{
R.Lock(i1);
R[i1] -= flux;
R.UnLock(i1);
}
if( s2 != Element::Ghost )
{
R.Lock(i2);
R[i2] += flux;
R.UnLock(i2);
}
}
}
}
#if defined(USE_OMP)
#pragma omp parallel for
#endif
<<<<<<< HEAD
for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell ) if( m->isValidCell(icell) )
{
Cell cell = Cell(m,ComposeCellHandle(icell));
......@@ -264,48 +316,87 @@ int main(int argc,char ** argv)
Tag error = m->CreateTag("error",DATA_REAL,CELL,NONE,1);
Storage::real err_C = 0.0, err_L2 = 0.0;
=======
for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell ) if( m->isValidCell(icell) )
{
Cell cell = Cell(m,ComposeCellHandle(icell));
if( cell->GetStatus() != Element::Ghost )
R[cell->Integer(id)] += cell->Mean(func_rhs, cell->Real(tensor_K)) * cell->Volume();
}
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Matrix assemble: " << Timer()-ttt << std::endl;
m->RemoveGeometricData(table); // Clean the computed geometric data
if( argc > 3 ) // Save the matrix and RHS if required
{
ttt = Timer();
R.GetJacobian().Save(std::string(argv[2])); // "A.mtx"
R.GetResidual().Save(std::string(argv[3])); // "b.rhs"
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Save matrix \"" << argv[2] << "\" and RHS \"" << argv[3] << "\": " << Timer()-ttt << std::endl;
}
ttt = Timer();
S.SetMatrix(R.GetJacobian()); // Compute the preconditioner for the original matrix
S.Solve(R.GetResidual(),Update); // Solve the linear system with the previously computted preconditioner
BARRIER;
if( m->GetProcessorRank() == 0 )
{
std::cout << S.Residual() << " " << S.Iterations() << " " << S.GetReason() << std::endl;
std::cout << "Solve system: " << Timer()-ttt << std::endl;
}
ttt = Timer();
Tag error = m->CreateTag("error",DATA_REAL,CELL,NONE,1);
Storage::real err_C = 0.0, err_L2 = 0.0;
>>>>>>> INMOST-DEV/master
#if defined(USE_OMP)
#pragma omp parallel
#endif
{
Storage::real local_err_C = 0;
{
Storage::real local_err_C = 0;
#if defined(USE_OMP)
#pragma omp for reduction(+:err_L2)
#endif
for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell )
{
Cell cell = Cell(m,ComposeCellHandle(icell));
if( cell->GetStatus() != Element::Ghost )
{
Storage::real old = cell->Real(phi);
Storage::real exact = cell->Mean(func, 0); // Compute the mean value of the function over the cell
Storage::real res = Update[aut.GetDynamicIndex(cell->self(),iphi)];
Storage::real sol = old-res;
Storage::real err = fabs (sol - exact);
if (err > local_err_C) local_err_C = err;
err_L2 += err * err * cell->Volume();
cell->Real(error) = err;
cell->Real(phi) = sol;
}
}
for( Storage::integer icell = 0; icell < m->CellLastLocalID(); ++icell )
{
Cell cell = Cell(m,ComposeCellHandle(icell));
if( cell->GetStatus() != Element::Ghost )
{
Storage::real old = cell->Real(phi);
Storage::real exact = cell->Mean(func, 0); // Compute the mean value of the function over the cell
Storage::real res = Update[aut.GetDynamicIndex(cell->self(),iphi)];
Storage::real sol = old-res;
Storage::real err = fabs (sol - exact);
if (err > local_err_C) local_err_C = err;
err_L2 += err * err * cell->Volume();
cell->Real(error) = err;
cell->Real(phi) = sol;
}
}
#if defined(USE_OMP)
#pragma omp critical
#endif
{
if( local_err_C > err_C ) err_C = local_err_C;
}
}
err_C = m->AggregateMax(err_C); // Compute the maximal C norm for the error
err_L2 = sqrt(m->Integrate(err_L2)); // Compute the global L2 norm for the error
if( m->GetProcessorRank() == 0 ) std::cout << "err_C = " << err_C << std::endl;
if( m->GetProcessorRank() == 0 ) std::cout << "err_L2 = " << err_L2 << std::endl;
}
BARRIER
{
if( local_err_C > err_C ) err_C = local_err_C;
}
}
err_C = m->AggregateMax(err_C); // Compute the maximal C norm for the error
err_L2 = sqrt(m->Integrate(err_L2)); // Compute the global L2 norm for the error
if( m->GetProcessorRank() == 0 ) std::cout << "err_C = " << err_C << std::endl;
if( m->GetProcessorRank() == 0 ) std::cout << "err_L2 = " << err_L2 << std::endl;
}
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Compute true residual: " << Timer()-ttt << std::endl;
ttt = Timer();
m->ExchangeData(phi,CELL,0); // Data exchange over processors
BARRIER
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Exchange phi: " << Timer()-ttt << std::endl;
std::string filename = "result";
......@@ -315,8 +406,8 @@ int main(int argc,char ** argv)
filename += ".pvtk";
ttt = Timer();
m->Save(filename);
m->Save("result.pmf");
BARRIER
m->Save("result.pmf");
BARRIER;
if( m->GetProcessorRank() == 0 ) std::cout << "Save \"" << filename << "\": " << Timer()-ttt << std::endl;
......
......@@ -54,7 +54,11 @@ int main(int argc,char ** argv)
std::cout << "J: "; J.Print();
return 0;
*/
<<<<<<< HEAD
if( argc > 2 )
=======
if( argc > 2 )
>>>>>>> INMOST-DEV/master
{
if( std::string(argv[2]) == "MFD" )
rt0 = false;
......@@ -83,8 +87,13 @@ int main(int argc,char ** argv)
std::cout << "Running MHFE RT0" << std::endl;
else
std::cout << "Running MFD" << std::endl;
<<<<<<< HEAD
=======
>>>>>>> INMOST-DEV/master
#if defined(USE_PARTITIONER)
if (m->GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner
{
......@@ -198,10 +207,17 @@ int main(int argc,char ** argv)
Mesh * mesh = m;