Commit dc4c4552 authored by Kirill Terekhov's avatar Kirill Terekhov

Tests into TestSuite

cavity and backward facing step
parent 592ec9b4
Pipeline #138 passed with stages
in 10 minutes and 24 seconds
......@@ -26,10 +26,43 @@ add_executable(test_transform_grid transform_grid.cpp)
add_executable(test_two_wells_3d two_wells_3d.cpp)
add_executable(test_edwards_3dtest4 edwards_3dtest4.cpp)
add_executable(test_discontinuous_4zones discontinuous_4zones.cpp)
add_executable(test_linear_sol linear_sol.cpp)
add_executable(test_linear_sol linear_sol.cpp)
add_executable(test_ns_linear ns_linear.cpp)
add_executable(test_ns_cavity ns_cavity.cpp)
add_executable(test_ns_backward_step ns_backward_step.cpp)
target_link_libraries(test_ns_linear inmost)
if(USE_MPI)
message("linking test_ns_linear with MPI")
target_link_libraries(test_ns_linear ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(test_ns_linear PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS test_ns_linear EXPORT projects-targets RUNTIME DESTINATION bin/TestSuite)
target_link_libraries(test_ns_cavity inmost)
if(USE_MPI)
message("linking test_ns_cavity with MPI")
target_link_libraries(test_ns_cavity ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(test_ns_cavity PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS test_ns_cavity EXPORT projects-targets RUNTIME DESTINATION bin/TestSuite)
target_link_libraries(test_ns_backward_step inmost)
if(USE_MPI)
message("linking test_ns_backward_step with MPI")
target_link_libraries(test_ns_backward_step ${MPI_LIBRARIES})
if(MPI_LINK_FLAGS)
set_target_properties(test_ns_backward_step PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}")
endif()
endif(USE_MPI)
install(TARGETS test_ns_backward_step EXPORT projects-targets RUNTIME DESTINATION bin/TestSuite)
target_link_libraries(test_single_well inmost)
......
#include "inmost.h"
using namespace INMOST;
//setup on mesh:
//
// FORCE - 3 entries
//
// BOUNDARY_CONDITION_VELOCITY - 7 entries
// r = (a5,a6,a7)
// n.(a1*u + a2*t) = n.r
// (I - nn.)(a3*u + a4*t) = (I-nn.)r
//
// BOUNDARY_CONDITION_PRESSURE - 1 entry
//
// REFERENCE_VELOCITY - 3 entries
//
// REFERENCE_PRESSURE - 1 entry
//
typedef Storage::real real;
int main(int argc, char ** argv)
{
if( argc < 2 )
{
std::cout << "Usage: " << argv[0] << " mesh is2D=1 height=0.6 length=0.2 mesh_out=grid_out.pmf" << std::endl;
return 0;
}
Mesh * m = new Mesh;
m->SetFileOption("VERBOSITY","2");
try
{
m->Load(argv[1]);
}
catch(...)
{
std::cout << "Cannot load the mesh " << argv[1] << std::endl;
return -1;
}
std::string fout = "grid_out.pmf";
int is2D = 1;
double height = 0.6;
double length = 0.2;
if( argc > 2 ) is2D = atoi(argv[2]);
if( argc > 3 ) height = atof(argv[3]);
if( argc > 4 ) length = atof(argv[4]);
if( argc > 5 ) fout = std::string(argv[5]);
//TagRealArray force = m->CreateTag("FORCE",DATA_REAL,CELL,NONE,3);
TagRealArray bc = m->CreateTag("BOUNDARY_CONDITION_VELOCITY",DATA_REAL,FACE,FACE,7);
TagReal bcp = m->CreateTag("BOUNDARY_CONDITION_PRESSURE",DATA_REAL,FACE,FACE,1);
TagRealArray uvw = m->CreateTag("UVW",DATA_REAL,CELL,NONE,3);
TagReal p = m->CreateTag("P",DATA_REAL,CELL,NONE,1);
//cut part of the grid
rMatrix n(3,1), x(3,1);
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
{
it->Centroid(x.data());
if( x(0,0) < length && x(1,0) < height )
it->Delete();
}
//this should not be needed?
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
it->FixNormalOrientation();
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
{
uvw(*it,3,1).Zero();
p[*it] = 0;
}
if( is2D ) std::cout << "2D setup" << std::endl;
else std::cout << "3D setup" << std::endl;
std::cout << "height=" << height << " length=" << length << std::endl;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
{
it->Centroid(x.data());
it->UnitNormal(n.data());
if( is2D && (fabs(n(2,0)-1) < 1.0-3 || fabs(n(2,0)+1) < 1.0e-3 ) ) //roller bc
{
bc[*it][0] = 1;
bc[*it][1] = 0;
bc[*it][2] = 0;
bc[*it][3] = 1;
bc[*it][4] = 0;
bc[*it][5] = 0;
bc[*it][6] = 0;
}
else if( fabs(n(0,0)+1) < 1.0e-3 && x(0,0) < 1.0e-3 ) //inlet
{
std::cout << "hi" << std::endl;
bc[*it][0] = 1;
bc[*it][1] = 0;
bc[*it][2] = 1;
bc[*it][3] = 0;
bc[*it][4] = 1;
bc[*it][5] = 0;
bc[*it][6] = 0;
}
else if( fabs(n(0,0)-1) < 1.0e-3 ) //outlet
bcp[*it] = 0;
else //no-slip
{
bc[*it][0] = 1;
bc[*it][1] = 0;
bc[*it][2] = 1;
bc[*it][3] = 0;
bc[*it][4] = 0;
bc[*it][5] = 0;
bc[*it][6] = 0;
}
}
std::cout << "Saving output to " << fout << std::endl;
m->Save(fout);
return 0;
}
#include "inmost.h"
using namespace INMOST;
//setup on mesh:
//
// FORCE - 3 entries
//
// BOUNDARY_CONDITION_VELOCITY - 7 entries
// r = (a5,a6,a7)
// n.(a1*u + a2*t) = n.r
// (I - nn.)(a3*u + a4*t) = (I-nn.)r
//
// BOUNDARY_CONDITION_PRESSURE - 1 entry
//
// REFERENCE_VELOCITY - 3 entries
//
// REFERENCE_PRESSURE - 1 entry
//
typedef Storage::real real;
int main(int argc, char ** argv)
{
if( argc < 2 )
{
std::cout << "Usage: " << argv[0] << " mesh is2D=1 mesh_out=grid_out.pmf" << std::endl;
return 0;
}
Mesh * m = new Mesh;
m->SetFileOption("VERBOSITY","2");
try
{
m->Load(argv[1]);
}
catch(...)
{
std::cout << "Cannot load the mesh " << argv[1] << std::endl;
return -1;
}
std::string fout = "grid_out.pmf";
int is2D = 1;
if( argc > 2 ) is2D = atoi(argv[2]);
if( argc > 3 ) fout = std::string(argv[3]);
//TagRealArray force = m->CreateTag("FORCE",DATA_REAL,CELL,NONE,3);
TagRealArray bc = m->CreateTag("BOUNDARY_CONDITION_VELOCITY",DATA_REAL,FACE,FACE,7);
//TagReal bcp = m->CreateTag("BOUNDARY_CONDITION_PRESSURE",DATA_REAL,FACE,FACE,1);
TagRealArray uvw = m->CreateTag("UVW",DATA_REAL,CELL,NONE,3);
TagReal p = m->CreateTag("P",DATA_REAL,CELL,NONE,1);
rMatrix n(3,1);
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
{
uvw(*it,3,1).Zero();
p[*it] = 0;
}
if( is2D ) std::cout << "2D setup" << std::endl;
else std::cout << "3D setup" << std::endl;
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it) if( it->Boundary() )
{
it->UnitNormal(n.data());
if( is2D && (fabs(n(2,0)-1) < 1.0-3 || fabs(n(2,0)+1) < 1.0e-3 ) ) //roller bc
{
bc[*it][0] = 1;
bc[*it][1] = 0;
bc[*it][2] = 0;
bc[*it][3] = 1;
bc[*it][4] = 0;
bc[*it][5] = 0;
bc[*it][6] = 0;
}
else
{
bc[*it][0] = 1;
bc[*it][1] = 0;
bc[*it][2] = 1;
bc[*it][3] = 0;
if( fabs(n(1,0)-1) < 1.0e-3 ) //top
{
bc[*it][4] = 1;
bc[*it][5] = 0;
bc[*it][6] = 0;
}
else
{
bc[*it][4] = 0;
bc[*it][5] = 0;
bc[*it][6] = 0;
}
}
}
std::cout << "Saving output to " << fout << std::endl;
m->Save(fout);
return 0;
}
#include "inmost.h"
using namespace INMOST;
//setup on mesh:
//
// FORCE - 3 entries
//
// BOUNDARY_CONDITION_VELOCITY - 7 entries
// r = (a5,a6,a7)
// n.(a1*u + a2*t) = n.r
// (I - nn.)(a3*u + a4*t) = (I-nn.)r
//
// BOUNDARY_CONDITION_PRESSURE - 1 entry
//
// REFERENCE_VELOCITY - 3 entries
//
// REFERENCE_PRESSURE - 1 entry
//
typedef Storage::real real;
vMatrix refUVW(real vx, real vy, real vz, real t, real nu)
{
unknown x(vx,0);
unknown y(vy,1);
unknown z(vz,2);
vMatrix UVW(3,1);
UVW(0,0) = 1*x+2*y+3*z;
UVW(1,0) = 4*x+5*y+6*z;
UVW(2,0) = 7*x+8*y+9*z;
return UVW;
}
variable refP(real x, real y, real z, real t, real nu)
{
return 0;
}
static real refq(real vx, real vy, real vz, real t, real nu, const rMatrix & n)
{
rMatrix UVW = refUVW(vx,vy,vz,t,nu);
return n.DotProduct(UVW);
}
static rMatrix reft(real vx, real vy, real vz, real t, real nu, real rho, const rMatrix & n)
{
static const rMatrix I = rMatrix::Unit(3);
rMatrix N;
//if( fullstress)
// N = 0.5*(I.Kronecker(n.Transpose())+n.Transpose().Kronecker(I));
//else
N = I.Kronecker(n.Transpose());
vMatrix UVW = refUVW(vx,vy,vz,t,nu);
rMatrix vUVW = UVW;
rMatrix G(9,1);
//real p = (addpressure || addpressure_analytic) ? get_value(refP(vx,vy,vz,t,nu)) : 0.0;
real p = get_value(refP(vx,vy,vz,t,nu));
for(int k = 0; k < 3; ++k)
for(int l = 0; l < 3; ++l)
G(l+3*k,0) = UVW(k,0).GetRow()[l];
rMatrix ret = rho*UVW*UVW.Transpose()*n - nu*N*G + p*n;
return ret;
}
int main(int argc, char ** argv)
{
if( argc < 2 )
{
std::cout << "Usage: " << argv[0] << " mesh time=1 nu=0.1 rho=1 mesh_out=grid_out.pmf" << std::endl;
return 0;
}
Mesh * m = new Mesh;
m->SetFileOption("VERBOSITY","2");
try
{
m->Load(argv[1]);
}
catch(...)
{
std::cout << "Cannot load the mesh " << argv[1] << std::endl;
return -1;
}
std::string fout = "grid_out.pmf";
double time = 1;
double nu = 0.1;
double rho = 1;
if( argc > 2 ) time = atof(argv[2]);
if( argc > 3 ) nu = atof(argv[3]);
if( argc > 4 ) rho = atof(argv[4]);
if( argc > 5 ) fout = std::string(argv[5]);
TagRealArray force = m->CreateTag("FORCE",DATA_REAL,CELL,NONE,3);
TagRealArray bc = m->CreateTag("BOUNDARY_CONDITION_VELOCITY",DATA_REAL,FACE,FACE,7);
//TagReal bcp = m->CreateTag("BOUNDARY_CONDITION_PRESSURE",DATA_REAL,FACE,FACE,1);
TagRealArray uvw = m->CreateTag("REFERENCE_VELOCITY",DATA_REAL,CELL,NONE,3);
TagReal p = m->CreateTag("REFERENCE_PRESSURE",DATA_REAL,CELL,NONE,1);
rMatrix x(3,1), n(3,1);
for(Mesh::iteratorCell it = m->BeginCell(); it != m->EndCell(); ++it)
{
it->Centroid(x.data());
uvw(*it,3,1) = refUVW(x(0,0),x(1,0),x(2,0),time,nu);
p[*it] = get_value(refP(x(0,0),x(1,0),x(2,0),time,nu));
ElementArray<Face> faces = it->getFaces();
force(*it,3,1).Zero();
for(ElementArray<Face>::iterator f = faces.begin(); f != faces.end(); ++f)
{
f->Barycenter(x.data());
f->OrientedUnitNormal(it->self(),n.data());
force(*it,3,1) += reft(x(0,0),x(1,0),x(2,0),time,nu,rho,n)*f->Area();
}
force(*it,3,1) /= it->Volume();
}
for(Mesh::iteratorFace it = m->BeginFace(); it != m->EndFace(); ++it)
{
it->Centroid(x.data());
it->UnitNormal(n.data());
if( it->Boundary() )
{
bc[*it][0] = 1;
bc[*it][1] = 0;
bc[*it][2] = 1;
bc[*it][3] = 0;
bc(*it,7,1)(4,7,0,1) = refUVW(x(0,0),x(1,0),x(2,0),time,nu);
}
}
std::cout << "Saving output to " << fout << std::endl;
m->Save(fout);
return 0;
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment