1 /// This is the incompressible Stokes solver benchmark based on this project:
2 /// http://www.featflow.de/en/benchmarks/cfdbenchmarking.html
3 ///
4 /// Author: A.Shamanskiy (2016 - ...., TU Kaiserslautern)
5 #include <gismo.h>
6 #include <gsElasticity/gsNsAssembler.h>
7 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
8 
9 using namespace gismo;
10 
refineBoundaryLayer(gsMultiBasis<> & velocity,gsMultiBasis<> & pressure)11 void refineBoundaryLayer(gsMultiBasis<> & velocity, gsMultiBasis<> & pressure)
12 {
13     gsMatrix<> boxSouth(2,2);
14     boxSouth << 0.,0.,0.,0.2;
15     for (index_t p = 0; p < 4; ++p)
16     {
17         velocity.refine(p,boxSouth);
18         pressure.refine(p,boxSouth);
19     }
20 }
21 
main(int argc,char * argv[])22 int main(int argc, char* argv[]){
23     gsInfo << "Benchmark 2D-0: steady-state flow of an incompressible super-viscous fluid.\n";
24 
25     //=====================================//
26                     // Input //
27     //=====================================//
28 
29     std::string filename = ELAST_DATA_DIR"/flow_around_cylinder.xml";
30     real_t viscosity = 0.001; // kinematic viscosity
31     real_t meanVelocity = 0.2; // inflow velocity
32     real_t density = 1.;
33     // space discretization
34     index_t numUniRef = 1;
35     index_t numDegElev = 0;
36     index_t numBLRef = 1;
37     bool subgridOrTaylorHood = false;
38     // output
39     index_t numPlotPoints = 10000;
40     bool plotMesh = false;
41 
42     // minimalistic user interface for terminal
43     gsCmdLine cmd("Testing the Stokes solver in 2D.");
44     cmd.addInt("r","refine","Number of uniform refinement applications",numUniRef);
45     cmd.addInt("d","degelev","Number of degree elevations",numDegElev);
46     cmd.addInt("l","blayer","Number of additional boundary layer refinements for the fluid",numBLRef);
47     cmd.addSwitch("e","element","Mixed element: false = subgrid (default), true = Taylor-Hood",subgridOrTaylorHood);
48     cmd.addInt("p","points","Number of points to plot to Paraview",numPlotPoints);
49     cmd.addSwitch("m","mesh","Plot computational mesh",plotMesh);
50     try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
51     gsInfo << "Using " << (subgridOrTaylorHood ? "Taylor-Hood " : "subgrid ") << "mixed elements.\n";
52 
53     //=============================================//
54         // Scanning geometry and creating bases //
55     //=============================================//
56 
57     // scanning geometry
58     gsMultiPatch<> geometry;
59     gsReadFile<>(filename, geometry);
60 
61     // creating bases
62     gsMultiBasis<> basisVelocity(geometry);
63     gsMultiBasis<> basisPressure(geometry);
64     for (index_t i = 0; i < numDegElev; ++i)
65     {
66         basisVelocity.degreeElevate();
67         basisPressure.degreeElevate();
68     }
69     for (index_t i = 0; i < numUniRef; ++i)
70     {
71         basisVelocity.uniformRefine();
72         basisPressure.uniformRefine();
73     }
74     // additional refinement of the boundary layer around the cylinder
75     for (index_t i = 0; i < numBLRef; ++i)
76         refineBoundaryLayer(basisVelocity,basisPressure);
77     // additional velocity refinement for stable mixed FEM
78     if (!subgridOrTaylorHood) // subgrid
79         basisVelocity.uniformRefine();
80     else // Taylor-Hood
81         basisVelocity.degreeElevate();
82 
83     //=============================================//
84         // Setting loads and boundary conditions //
85     //=============================================//
86 
87     // inflow velocity profile U(y) = U_max*y*(H-y)/(H/2)^2; channel height H = 0.41
88     gsFunctionExpr<> inflow(util::to_string(meanVelocity) + "*6*y*(0.41-y)/0.41^2",2);
89 
90     // boundary conditions
91     gsBoundaryConditions<> bcInfo;
92     bcInfo.addCondition(0,boundary::north,condition_type::dirichlet,&inflow,0);
93     bcInfo.addCondition(0,boundary::north,condition_type::dirichlet,0,1);
94     for (index_t d = 0; d < 2; ++d)
95     {   // no slip conditions
96         bcInfo.addCondition(0,boundary::south,condition_type::dirichlet,0,d);
97         bcInfo.addCondition(1,boundary::south,condition_type::dirichlet,0,d);
98         bcInfo.addCondition(1,boundary::north,condition_type::dirichlet,0,d);
99         bcInfo.addCondition(2,boundary::south,condition_type::dirichlet,0,d);
100         bcInfo.addCondition(3,boundary::south,condition_type::dirichlet,0,d);
101         bcInfo.addCondition(3,boundary::north,condition_type::dirichlet,0,d);
102         bcInfo.addCondition(4,boundary::south,condition_type::dirichlet,0,d);
103         bcInfo.addCondition(4,boundary::north,condition_type::dirichlet,0,d);
104     }
105 
106     // source function, rhs
107     gsConstantFunction<> g(0.,0.,2);
108 
109     //=============================================//
110               // Assembling & solving //
111     //=============================================//
112 
113     // creating assembler
114     gsNsAssembler<real_t> assembler(geometry,basisVelocity,basisPressure,bcInfo,g);
115     assembler.options().setReal("Viscosity",viscosity);
116     assembler.options().setReal("Density",density);
117     assembler.options().setInt("DirichletValues",dirichlet::interpolation);
118     gsInfo<<"Assembling...\n";
119     gsStopwatch clock;
120     clock.restart();
121     assembler.assemble();
122     gsInfo << "Assembled a system (matrix and load vector) with "
123            << assembler.numDofs() << " dofs in " << clock.stop() << "s.\n";
124 
125     gsInfo << "Solving...\n";
126     clock.restart();
127 
128 #ifdef GISMO_WITH_PARDISO
129     gsSparseSolver<>::PardisoLDLT solver(assembler.matrix());
130     gsVector<> solVector = solver.solve(assembler.rhs());
131     gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() <<"s.\n";
132 #else
133     gsSparseSolver<>::SimplicialLDLT solver(assembler.matrix());
134     gsVector<> solVector = solver.solve(assembler.rhs());
135     gsInfo << "Solved the system with EigenLDLT solver in " << clock.stop() <<"s.\n";
136 #endif
137 
138     //=============================================//
139                       // Output //
140     //=============================================//
141 
142     // constructing solution as an IGA function
143     gsMultiPatch<> velocity, pressure;
144     assembler.constructSolution(solVector,assembler.allFixedDofs(),velocity,pressure);
145 
146     if (numPlotPoints > 0) // visualization
147     {
148         // constructing isogeometric field (geometry + solution)
149         gsField<> velocityField(assembler.patches(),velocity);
150         gsField<> pressureField(assembler.patches(),pressure);
151         // creating a container to plot all fields to one Paraview file
152         std::map<std::string,const gsField<> *> fields;
153         fields["Velocity"] = &velocityField;
154         fields["Pressure"] = &pressureField;
155         gsWriteParaviewMultiPhysics(fields,"aroundCylinder",numPlotPoints,plotMesh);
156         gsInfo << "Open \"aroundCylinder.pvd\" in Paraview for visualization.\n";
157     }
158 
159     // computing forces acting on the surface of the solid body
160     std::vector<std::pair<index_t, boxSide> > bdrySides;
161     bdrySides.push_back(std::pair<index_t,index_t>(0,boxSide(boundary::south)));
162     bdrySides.push_back(std::pair<index_t,index_t>(1,boxSide(boundary::south)));
163     bdrySides.push_back(std::pair<index_t,index_t>(2,boxSide(boundary::south)));
164     bdrySides.push_back(std::pair<index_t,index_t>(3,boxSide(boundary::south)));
165     gsMatrix<> force = assembler.computeForce(velocity,pressure,bdrySides);
166 
167     // evaluating pressure difference at the far front and the far rear points of the cylinder
168     gsMatrix<> point(2,1);
169     point << 0.5, 0;
170     real_t pressureFront = pressure.patch(0).eval(point)(0,0);
171     real_t pressureBack = pressure.patch(2).eval(point)(0,0);
172     real_t L = 0.1; // characteristic length
173 
174     gsInfo << "Drag coefficient: " << 2.*force.at(0)/L/pow(meanVelocity,2) << std::endl;
175     gsInfo << "Lift coefficient: " << 2.*force.at(1)/L/pow(meanVelocity,2) << std::endl;
176     gsInfo << "Pressure difference: " << pressureFront - pressureBack << std::endl;
177 
178     return 0;
179 }
180