Commit a71f2eaf authored by Kirill Terekhov's avatar Kirill Terekhov

Ethier-Steinman test for MFD

Added test for MFD scheme as projection step for Navier-Stokes problem
with analytic solution of Ethier-Steinman
parent 4a3d3ca2
......@@ -15,6 +15,7 @@ endif(USE_SOLVER AND USE_MESH)
if(USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
add_subdirectory(ADFVDiscr)
add_subdirectory(ADMFD)
add_subdirectory(MFD-ES)
endif(USE_AUTODIFF AND USE_SOLVER AND USE_MESH)
#add_subdirectory(OctreeCutcell)
if(USE_MESH AND USE_PARTITIONER)
......
project(MFDES)
set(SOURCE main.cpp)
add_executable(MFDES ${SOURCE})
target_link_libraries(MFDES inmost)
if(USE_SOLVER)
if(USE_SOLVER_ANI)
message("linking MFDES with ani3d and BLAS")
target_link_libraries(MFDES ani3d ${BLAS_LIBRARIES})
if(BLAS_LINKER_FLAGS)
set_target_properties(MFDES PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}")
endif()
endif()
if(USE_SOLVER_PETSC)
message("linking MFDES with PETSc")
target_link_libraries(MFDES ${PETSC_LIBRARIES})
endif()
if(USE_SOLVER_TRILINOS)
message("linking MFDES with Trilinos")
target_link_libraries(MFDES ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_SOLVER_METIS)
message("linking MFDES with Metis")
target_link_libraries(MFDES ${METIS_LIBRARIES})
endif()
if(USE_SOLVER_MONDRIAAN)
message("linking MFDES with Mondriaan")
target_link_libraries(MFDES ${MONDRIAAN_LIBRARIES})
endif()
if(USE_SOLVER_SUPERLU)
message("linking MFDES with SuperLU")
target_link_libraries(MFDES ${SUPERLU_LIBRARIES})
endif()
endif()
if(USE_PARTITIONER)
if(USE_PARTITIONER_ZOLTAN)
message("linking MFDES with Zoltan")
target_link_libraries(MFDES ${ZOLTAN_LIBRARIES})
endif()
if(USE_PARTITIONER_PARMETIS)
message("linking MFDES with ParMETIS")
target_link_libraries(MFDES ${PARMETIS_LIBRARIES})
endif()
endif()
if(USE_MPI)
message("linking MFDES with MPI")
target_link_libraries(MFDES ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(MFDES PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS MFDES EXPORT inmost-targets RUNTIME DESTINATION bin)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "inmost.h"
using namespace INMOST;
#ifndef M_PI
#define M_PI 3.141592653589
#endif
#if defined(USE_MPI)
#define BARRIER MPI_Barrier(MPI_COMM_WORLD);
#else
#define BARRIER
#endif
bool report_residual = false;
bool write_steps = false;
bool two_point = false;
//shortcuts
typedef Storage::bulk bulk;
typedef Storage::real real;
typedef Storage::integer integer;
typedef Storage::enumerator enumerator;
typedef Storage::real_array real_array;
typedef Storage::var_array var_array;
variable refU(real xyz[3], real nu, real t)
{
unknown x(xyz[0],0);
unknown y(xyz[1],1);
unknown z(xyz[2],2);
const real a = M_PI/4.0;
const real d = M_PI/2.0;
return -a*(exp(a*x)*sin(a*y+d*z)+exp(a*z)*cos(a*x+d*y))*exp(-nu*d*d*t);
}
variable refV(real xyz[3], real nu, real t)
{
unknown x(xyz[0],0);
unknown y(xyz[1],1);
unknown z(xyz[2],2);
const real a = M_PI/4.0;
const real d = M_PI/2.0;
return -a*(exp(a*y)*sin(a*z+d*x)+exp(a*x)*cos(a*y+d*z))*exp(-nu*d*d*t);
}
variable refW(real xyz[3], real nu, real t)
{
unknown x(xyz[0],0);
unknown y(xyz[1],1);
unknown z(xyz[2],2);
const real a = M_PI/4.0;
const real d = M_PI/2.0;
return -a*(exp(a*z)*sin(a*x+d*y)+exp(a*y)*cos(a*z+d*x))*exp(-nu*d*d*t);
}
real refP(real xyz[3], real nu, real t)
{
real x = xyz[0];
real y = xyz[1];
real z = xyz[2];
const real a = M_PI/4.0;
const real d = M_PI/2.0;
return -a*a/2.0*(exp(2*a*x)+exp(2*a*y)+exp(2*a*z)+
2*sin(a*x+d*y)*cos(a*z+d*x)*exp(a*(y+z))+
2*sin(a*y+d*z)*cos(a*x+d*y)*exp(a*(z+x))+
2*sin(a*z+d*x)*cos(a*y+d*z)*exp(a*(x+y)))*exp(-2*nu*d*d*t);
}
// (done)todo: adapt for parallel use
void CheckResidual(Mesh & m, const TagReal & tag_P, const TagReal & tag_nU, real visc, real T)
{
{
TagReal tag_refP = m.CreateTag("refP" ,DATA_REAL,CELL|FACE,NONE,1);
real p_ref_min = 1.0e20, p_ref_max = -1.0e20;
real p_min = 1.0e20, p_max = -1.0e20;
for(ElementType etype = FACE; etype <= CELL; etype = NextElementType(etype) )
for(integer it = 0; it < m.LastLocalID(etype); ++it) if( m.isValidElement(etype,it) )
{
Element e = m.ElementByLocalID(etype,it);
real x[3],p;
e.Centroid(x);
p = refP(x,visc,T);
if( p < p_ref_min ) p_ref_min = p;
if( p > p_ref_max ) p_ref_max = p;
p = tag_P[e];
if( p < p_min ) p_min = p;
if( p > p_max ) p_max = p;
}
p_max = m.AggregateMax(p_max);
p_min = m.AggregateMin(p_min);
p_ref_max = m.AggregateMax(p_ref_max);
p_ref_min = m.AggregateMin(p_ref_min);
for(ElementType etype = FACE; etype <= CELL; etype = NextElementType(etype) )
for(integer it = 0; it < m.LastLocalID(etype); ++it) if( m.isValidElement(etype,it) )
{
Element e = m.ElementByLocalID(etype,it);
real x[3];
e.Centroid(x);
tag_refP[e] = refP(x,visc,T) - p_ref_min;
}
if( m.GetProcessorRank() == 0 )
{
std::cout << "p_ref " << p_ref_min << ":" << p_ref_max << std::endl;
std::cout << "p " << p_min << ":" << p_max << std::endl;
}
TagReal tag_E = m.CreateTag("ERROR_P" ,DATA_REAL,CELL|FACE,NONE,1);
TagRealArray tag_EU = m.CreateTag("ERROR_U" ,DATA_REAL,CELL,NONE,3);
TagRealArray tag_U = m.CreateTag("U" ,DATA_REAL,CELL,NONE,3);
TagReal tag_EnU = m.CreateTag("ERROR_nU",DATA_REAL,FACE,NONE,1);
real C, L2, volume, Cu[3], L2u[3], Cnu, L2nu;
C = L2 = volume = 0.0;
Cu[0] = Cu[1] = Cu[2] = L2u[0] = L2u[1] = L2u[2] = 0;
for(integer q = 0; q < m.CellLastLocalID(); ++q ) if( m.isValidCell(q) )
{
Cell c = m.CellByLocalID(q);
if( c.GetStatus() == Element::Ghost ) continue;
real x[3], z[3];
c.Centroid(x);
real err = tag_P[c] - p_min - tag_refP[c];
real vol = c.Volume();
if( C < fabs(err) ) C = fabs(err);
L2 += err*err*vol;
volume += vol;
tag_E[c] = err;
real V[3] = {0,0,0}; //restored velocity
real a,s,v;
ElementArray<Face> cfaces = c.getFaces();
for(integer kt = 0; kt < cfaces.size(); ++kt)
{
cfaces[kt].Centroid(z);
a = cfaces[kt].Area();
s = cfaces[kt].FaceOrientedOutside(c) ? 1.0 : -1.0;
v = tag_nU[cfaces[kt]]*s*a;
V[0] += v*(z[0]-x[0]);
V[1] += v*(z[1]-x[1]);
V[2] += v*(z[2]-x[2]);
}
V[0] /= c.Volume();
V[1] /= c.Volume();
V[2] /= c.Volume();
tag_U[c][0] = V[0];
tag_U[c][1] = V[1];
tag_U[c][2] = V[2];
err = V[0] - refU(x,visc,T).GetValue();
if( Cu[0] < fabs(err) ) Cu[0] = fabs(err);
L2u[0] += err*err*vol;
tag_EU[c][0] = err;
err = V[1] - refV(x,visc,T).GetValue();
if( Cu[1] < fabs(err) ) Cu[1] = fabs(err);
L2u[1] += err*err*vol;
tag_EU[c][1] = err;
err = V[2] - refW(x,visc,T).GetValue();
if( Cu[2] < fabs(err) ) Cu[2] = fabs(err);
L2u[2] += err*err*vol;
tag_EU[c][2] = err;
}
m.ExchangeData(tag_U,CELL,0);
m.ExchangeData(tag_EU,CELL,0);
volume = m.Integrate(volume);
L2 = sqrt(m.Integrate(L2)/volume);
m.Integrate(L2u,3);
L2u[0] = sqrt(L2u[0]/volume);
L2u[1] = sqrt(L2u[1]/volume);
L2u[2] = sqrt(L2u[2]/volume);
C = m.AggregateMax(C);
m.AggregateMax(Cu,3);
if( m.GetProcessorRank() == 0 )
{
std::cout << "Pressure error on cells, C-norm " << C << " L2-norm " << L2 << std::endl;
std::cout << "Reconstructed U error on cells, C-norm " << Cu[0] << " L2-norm " << L2u[0] << std::endl;
std::cout << "Reconstructed V error on cells, C-norm " << Cu[1] << " L2-norm " << L2u[1] << std::endl;
std::cout << "Reconstructed W error on cells, C-norm " << Cu[2] << " L2-norm " << L2u[2] << std::endl;
}
C = L2 = volume = 0.0;
Cnu = L2nu = 0.0;
for( int q = 0; q < m.FaceLastLocalID(); ++q ) if( m.isValidFace(q) )
{
Face f = m.FaceByLocalID(q);
if( f.GetStatus() == Element::Ghost ) continue;
real x[3],y[3],n[3];
f.Centroid(x);
f.UnitNormal(n);
f.BackCell().Centroid(y);
real H = n[0]*(x[0]-y[0])+n[1]*(x[1]-y[1])+n[2]*(x[2]-y[2]);
if( f.FrontCell().isValid() )
{
f.FrontCell().Centroid(y);
H += n[0]*(y[0]-x[0])+n[1]*(y[1]-x[1])+n[2]*(y[2]-x[2]);
}
real err = tag_P[f] - p_min - tag_refP[f];
real vol = f.Area()*H;
if( C < fabs(err) ) C = fabs(err);
L2 += err*err*vol;
volume += vol;
tag_E[f] = err;
real V[3], nV;
V[0] = refU(x,visc,T).GetValue();
V[1] = refV(x,visc,T).GetValue();
V[2] = refW(x,visc,T).GetValue();
nV = n[0]*V[0]+n[1]*V[1]+n[2]*V[2];
err = tag_nU[f] - nV;
if( Cnu < fabs(err) ) Cnu = fabs(err);
L2nu += err*err*vol;
volume += vol;
tag_EnU[f] = err;
}
volume = m.Integrate(volume);
L2 = sqrt(m.Integrate(L2)/volume);
L2nu = sqrt(m.Integrate(L2nu)/volume);
C = m.AggregateMax(C);
Cnu = m.AggregateMax(Cnu);
if( m.GetProcessorRank() == 0 )
{
if( !two_point ) std::cout << "Pressure error on faces, C-norm " << C << " L2-norm " << L2 << std::endl;
std::cout << "Normal velocity error on faces, C-norm " << Cnu << " L2-norm " << L2nu << std::endl;
}
real div_max = -1.0e20, div_min = 1.0e20, div_int = 0;
for(integer it = 0; it < m.CellLastLocalID(); ++it) if( m.isValidCell(it) ) //loop over cells
{
Cell c = m.CellByLocalID(it);
if( c.GetStatus() == Element::Ghost ) continue;
ElementArray<Face> cfaces = c.getFaces(); //obtain faces of the cell
real div = 0, a, s;
for(integer kt = 0; kt < cfaces.size(); ++kt)
{
a = cfaces[kt].Area();
s = cfaces[kt].FaceOrientedOutside(c)?1.0:-1.0;
div += tag_nU[cfaces[kt]]*a*s;
}
if( div > div_max ) div_max = div;
if( div < div_min ) div_min = div;
div_int += div;
}
div_max = m.AggregateMax(div_max);
div_min = m.AggregateMin(div_min);
div_int = m.Integrate(div_int);
if( m.GetProcessorRank() == 0 )
{
std::cout << "Divergence of velocity: " << div_min << ":" << div_max << " integral " << div_int << std::endl;
}
}
}
int main(int argc,char ** argv)
{
Solver::Initialize(&argc,&argv,""); // Initialize the solver and MPI activity
#if defined(USE_PARTITIONER)
Partitioner::Initialize(&argc,&argv); // Initialize the partitioner activity
#endif
if( argc > 1 )
{
//parameters
real dT = 1.0/50.0;
real T = 1.0/10.0;
real visc = 0.0;
if( argc > 2 ) dT = atof(argv[2]);
if( argc > 3 ) T = atof(argv[3]);
if( argc > 4 ) visc = atof(argv[4]);
if( visc > 0 )
{
std::cout << "Viscosity is not calculated right now." << std::endl;
visc = 0;
}
bool repartition = false; // Is it required to redistribute the mesh?
Mesh m; // Create an empty mesh
{ // Load the mesh
m.SetCommunicator(INMOST_MPI_COMM_WORLD); // Set the MPI communicator for the mesh
if( m.GetProcessorRank() == 0 ) std::cout << "Processors: " << m.GetProcessorsNumber() << std::endl;
if( m.isParallelFileFormat(argv[1]) ) //The format is
{
m.Load(argv[1]); // Load mesh from the parallel file format
repartition = true; // Ask to repartition the mesh
}
else if( m.GetProcessorRank() == 0 ) m.Load(argv[1]); // Load mesh from the serial file format
}
#if defined(USE_PARTITIONER)
if (m.GetProcessorsNumber() > 1 )//&& !repartition) // Currently only non-distributed meshes are supported by Inner_RCM partitioner
{
{ // Compute mesh partitioning
Partitioner p(&m); //Create Partitioning object
p.SetMethod(Partitioner::Inner_RCM,repartition ? Partitioner::Repartition : Partitioner::Partition); // Specify the partitioner
p.Evaluate(); // Compute the partitioner and store new processor ID in the mesh
}
{ //Distribute the mesh
m.Redistribute(); // Redistribute the mesh data
m.ReorderEmpty(CELL|FACE|EDGE|NODE); // Clean the data after reordring
}
}
#endif
//map grid into [-1,1] box
real xmax[3] = {-1.0e20,-1.0e20,-1.0e20};
real xmin[3] = {+1.0e20,+1.0e20,+1.0e20};
for(integer it = 0; it < m.NodeLastLocalID(); ++it) if( m.isValidNode(it) )
{
real_array c = m.NodeByLocalID(it).Coords();
for(int k = 0; k < 3; ++k)
{
if( xmin[k] > c[k] ) xmin[k] = c[k];
if( xmax[k] < c[k] ) xmax[k] = c[k];
}
}
m.AggregateMax(xmax,3);
m.AggregateMin(xmin,3);
for(integer it = 0; it < m.NodeLastLocalID(); ++it) if( m.isValidNode(it) )
{
real_array c = m.NodeByLocalID(it).Coords();
for(int k = 0; k < 3; ++k)
c[k] = 2.0*(c[k]-xmin[k])/(xmax[k]-xmin[k])-1.0;
}
{ // prepare geometrical data on the mesh
Mesh::GeomParam table;
table[CENTROID] = CELL | FACE; //Compute averaged center of mass
table[NORMAL] = FACE; //Compute normals
table[ORIENTATION] = FACE; //Check and fix normal orientation
table[MEASURE] = CELL | FACE; //Compute volumes and areas
//table[BARYCENTER] = CELL | FACE; //Compute volumetric center of mass
m.PrepareGeometricData(table); //Ask to precompute the data
}
// data tags for
TagReal tag_nU; // Normal component of velocity
TagReal tag_nUold; // on previous step
TagReal tag_P; // Pressure
TagReal tag_Q; // Pressure update
TagReal tag_F; // Forcing term
TagRealArray tag_W; // Gradient matrix
TagInteger tag_i; // row number for current cell in back cell matrix
if( m.GetProcessorsNumber() > 1 ) //skip for one processor job
{ // Exchange ghost cells
m.ExchangeGhost(1,FACE); // Produce layer of ghost cells
}
{ //initialize data
tag_P = m.CreateTag("PRESSURE",DATA_REAL,CELL|FACE,NONE,1); // Create a new tag for the pressure
tag_W = m.CreateTag("W",DATA_REAL,CELL,NONE);
tag_i = m.CreateTag("row",DATA_INTEGER,FACE,NONE,1);
tag_F = m.CreateTag("FORCE",DATA_REAL,CELL,NONE,1);
tag_Q = m.CreateTag("PUPDATE",DATA_REAL,CELL|FACE,NONE,1);
tag_nU = m.CreateTag("NORMAL_VELOCITY",DATA_REAL,FACE,NONE,1);
tag_nUold = m.CreateTag("NORMAL_VELOCITY_OLD",DATA_REAL,FACE,NONE,1);
//Assemble gradient matrix W on cells
#if defined(USE_OMP)
#pragma omp parallel for
#endif
for(integer it = 0; it < m.CellLastLocalID(); ++it ) if( m.isValidCell(it) )
{
Cell c = m.CellByLocalID(it);
ElementArray<Face> faces = c.getFaces(); //obtain faces of the cell
tag_W[c].resize(faces.size()*faces.size());
Matrix<real,real_array> W(tag_W[c],faces.size(),faces.size());
real xP[3]; //center of the cell
real yF[3]; //center of the face
real nF[3]; //normal to the face
real aF; //area of the face
real vP = c.Volume(); //volume of the cell
c.Centroid(xP);
rMatrix N(faces.size(),3), R(faces.size(),3); //big gradient matrix, co-normals, directions
for(int k = 0; k < faces.size(); ++k) //loop over faces
{
aF = faces[k].Area();
faces[k].Centroid(yF);
faces[k].OrientedUnitNormal(c,nF);
// assemble matrix of directions
R(k,0) = (yF[0]-xP[0])*aF;
R(k,1) = (yF[1]-xP[1])*aF;
R(k,2) = (yF[2]-xP[2])*aF;
// assemble matrix of co-normals
N(k,0) = nF[0];
N(k,1) = nF[1];
N(k,2) = nF[2];
// remember position
if( faces[k].BackCell() == c )
tag_i[faces[k]] = k;
} //end of loop over faces
W = N*(N.Transpose()*R).Invert(true).first*N.Transpose(); //stability part
W += (rMatrix::Unit(faces.size()) - R*(R.Transpose()*R).Invert(true).first*R.Transpose())*(4.0/faces.size()*W.Trace());
} //end of loop over cells
//initialize normal velocity
#if defined(USE_OMP)
#pragma omp parallel for
#endif
for(integer it = 0; it < m.FaceLastLocalID(); ++it) if( m.isValidFace(it) )
{
Face f = m.FaceByLocalID(it);
real x[3], n[3], u[3];
f.UnitNormal(n);
f.Centroid(x);
u[0] = refU(x,visc,0).GetValue();
u[1] = refV(x,visc,0).GetValue();
u[2] = refW(x,visc,0).GetValue();
tag_nUold[f] = tag_nU[f] = u[0]*n[0]+u[1]*n[1]+u[2]*n[2];
//pressure
tag_P[f] = refP(x,visc,0);
}
//initialize pressure
#if defined(USE_OMP)
#pragma omp parallel for
#endif
for(integer it = 0; it < m.CellLastLocalID(); ++it) if( m.isValidCell(it) )
{
Cell c = m.CellByLocalID(it);
real x[3];
c.Centroid(x);
tag_P[c] = refP(x,visc,0);
}
} //end of initialize data
MarkerType boundary = m.CreateMarker();
m.MarkBoundaryFaces(boundary);
if( m.GetProcessorRank() == 0 ) std::cout << "Initialization done" << std::endl;
{ //Main loop for problem solution
Automatizator aut; // declare class to help manage unknowns
dynamic_variable Q(aut,aut.RegisterTag(tag_Q,two_point ? CELL : (CELL|FACE))); //register pressure as primary unknown
aut.EnumerateTags(); //enumerate all primary variables
std::cout << "Enumeration done, size " << aut.GetLastIndex() - aut.GetFirstIndex() << std::endl;
Residual R("",aut.GetFirstIndex(),aut.GetLastIndex());
Sparse::LockService Locks(aut.GetFirstIndex(),aut.GetLastIndex());
Sparse::AnnotationService Text(aut.GetFirstIndex(),aut.GetLastIndex());
Sparse::Vector Update ("",aut.GetFirstIndex(),aut.GetLastIndex()); //vector for update
{//Annotate matrix
for( int q = 0; q < m.CellLastLocalID(); ++q ) if( m.isValidCell(q) )
{
Cell cell = m.CellByLocalID(q);
if( cell.GetStatus() != Element::Ghost )
Text.SetAnnotation(Q.Index(cell),"Cell-centered pressure value");
}
if( !two_point ) for( int q = 0; q < m.FaceLastLocalID(); ++q ) if( m.isValidFace(q) )
{
Face face = m.FaceByLocalID(q);
if( face.GetStatus() != Element::Ghost )
{
//if( tag_BC.isValid() && face.HaveData(tag_BC) )
// Text.SetAnnotation(P.Index(face),"Pressure guided by boundary condition");
//else
Text.SetAnnotation(Q.Index(face),"Interface pressure");
}
}
}
if( m.GetProcessorRank() == 0 ) std::cout << "Matrix was annotated" << std::endl;
//checks
/*
for(integer it = 0; it < m.FaceLastLocalID(); ++it) if( m.isValidFace(it) )
{
Face f = m.FaceByLocalID(it);
Cell c = f.BackCell();
ElementArray<Face> cfaces = c.getFaces();
Matrix<real,real_array> W(tag_W[c],cfaces.size(),cfaces.size());
real ngrad[3] = {0,0,0};
real n[3], xf[3], xc[3];
f.UnitNormal(n);
c.Centroid(xc);
for(integer kt = 0; kt < cfaces.size(); ++kt)
{
f.Centroid(xf);
ngrad[0] += W(tag_i[f],kt)*(xf[0] - xc[0])*cfaces[kt].Area();
ngrad[1] += W(tag_i[f],kt)*(xf[1] - xc[1])*cfaces[kt].Area();
ngrad[2] += W(tag_i[f],kt)*(xf[2] - xc[2])*cfaces[kt].Area();
}
std::cout << "1/area: " << 1.0/f.Area() << std::endl;
std::cout << "normal: " << n[0] << " " << n[1] << " " << n[2] << std::endl;
std::cout << "ngrad: " << ngrad[0] << " " << ngrad[1] << " " << ngrad[2] << std::endl;
}
*/
int step = 0;
real alpha = 1, beta = -1, gamma = 0;
for(real t = 0; t < T-1.0e-9; t+=dT)
{
//if( step == 1 ) break;
if( t > T ) t = T;
if( report_residual )
{
if( m.GetProcessorRank() == 0 ) std::cout << "Before velocity update " << std::endl;
CheckResidual(m,tag_P,tag_nU,visc,t);
}
#if defined(USE_OMP)
#pragma omp parallel for
#endif
for(integer it = 0; it < m.FaceLastLocalID(); ++it) if( m.isValidFace(it) )
{
Face f = m.FaceByLocalID(it);
real x[3], n[3], V[3], nV, U, Uold;
f.UnitNormal(n);
f.Centroid(x);
U = tag_nU[f];
Uold = tag_nUold[f];
if( f.GetMarker(boundary) ) //Boundary condition
{
f.Centroid(x);
V[0] = refU(x,visc,t+dT).GetValue();
V[1] = refV(x,visc,t+dT).GetValue();
V[2] = refW(x,visc,t+dT).GetValue();
nV = n[0]*V[0]+n[1]*V[1]+n[2]*V[2];
tag_nU[f] = nV;
tag_nUold[f] = U;
}
else //Update with advection and pressure gradient
{
real ngradP = 0;
if( two_point )
{
if( f.FrontCell().isValid() )
{
real xf[3],xb[3],n[3],d;
f.UnitNormal(n);
f.FrontCell().Centroid(xf);
f.BackCell().Centroid(xb);
d = n[0]*(xf[0]-xb[0]) + n[1]*(xf[1]-xb[1]) + n[2]*(xf[2]-xb[2]);
ngradP = (tag_P[f.FrontCell()] - tag_P[f.BackCell()])/d;
}
}
else
{
Cell c = f.BackCell();
ElementArray<Face> cfaces = c.getFaces();
Matrix<real,real_array> W(tag_W[c],cfaces.size(),cfaces.size());
for(integer kt = 0; kt < cfaces.size(); ++kt)
ngradP += W(tag_i[f],kt)*(tag_P[cfaces[kt]] - tag_P[c])*cfaces[kt].Area();
}
variable u = refU(x,visc,t+dT);
variable v = refV(x,visc,t+dT);
variable w = refW(x,visc,t+dT);
V[0] = u.GetValue()*u.GetRow()[0] + v.GetValue()*u.GetRow()[1] + w.GetValue()*u.GetRow()[2];
V[1] = u.GetValue()*v.GetRow()[0] + v.GetValue()*v.GetRow()[1] + w.GetValue()*v.GetRow()[2];
V[2] = u.GetValue()*w.GetRow()[0] + v.GetValue()*w.GetRow()[1] + w.GetValue()*w.GetRow()[2];
nV = n[0]*V[0]+n[1]*V[1]+n[2]*V[2];
tag_nU[f] = -(beta*U + gamma*Uold + (nV + ngradP)*dT)/alpha;
tag_nUold[f] = U;
}
}
m.ExchangeData(tag_nU,FACE,0);
if( report_residual )
{
if( m.GetProcessorRank() == 0 ) std::cout << "After velocity update " << std::endl;
CheckResidual(m,tag_P,tag_nU,visc,t+dT);
}
R.Clear(); //clean up the residual
//clean up Q just in case, check does this affect the result
//for(ElementType etype = FACE; etype <= CELL; etype = NextElementType(etype) )
// for(integer it = 0; it < m.LastLocalID(etype); ++it) if( m.isValidElement(etype,it) )
// tag_Q[m.ElementByLocalID(etype,it)] = 0;
Automatizator::MakeCurrent(&aut);
#if defined(USE_OMP)
#pragma omp parallel
#endif
for(integer it = 0; it < m.CellLastLocalID(); ++it) if( m.isValidCell(it) ) //loop over cells
{
Cell c = m.CellByLocalID(it);
ElementArray<Face> cfaces = c.getFaces(); //obtain faces of the cell
vMatrix FLUX(cfaces.size(),1); //computed flux on faces
if( two_point )
{
for(integer kt = 0; kt < cfaces.size(); ++kt)
{
Face f = cfaces[kt];
if( f.FrontCell().isValid() )
{
real xf[3],xb[3],n[3],d,s;
f.UnitNormal(n);
f.FrontCell().Centroid(xf);
f.BackCell().Centroid(xb);