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