PolyFEM
Loading...
Searching...
No Matches
TransientNavierStokesSolver.cpp
Go to the documentation of this file.
2
4#include <polysolve/linear/FEMSolver.hpp>
6
8
9#include <igl/Timer.h>
10
11#include <unsupported/Eigen/SparseExtra>
12
13#include <cmath>
14
15namespace polyfem
16{
17 using namespace polysolve;
18 using namespace assembler;
19 using namespace utils;
20
21 namespace solver
22 {
23
25 : solver_param(solver_param)
26 {
27 gradNorm = solver_param["nonlinear"]["grad_norm"];
28 iterations = solver_param["nonlinear"]["max_iterations"];
29 }
30
32 const int n_bases,
33 const int n_pressure_bases,
34 const double t,
35 const std::vector<basis::ElementBases> &bases,
36 const std::vector<basis::ElementBases> &gbases,
37 assembler::NavierStokesVelocity &velocity_assembler,
38 const assembler::AssemblyValsCache &ass_vals_cache,
39 const std::vector<int> &boundary_nodes,
40 const bool use_avg_pressure,
41 const int problem_dim,
42 const bool is_volume,
43 const double beta_dt, const Eigen::VectorXd &prev_sol,
44 const StiffnessMatrix &velocity_stiffness, const StiffnessMatrix &mixed_stiffness, const StiffnessMatrix &pressure_stiffness,
45 const StiffnessMatrix &velocity_mass1,
46 const Eigen::MatrixXd &rhs, Eigen::VectorXd &x)
47 {
48 assert(velocity_assembler.name() == "NavierStokes");
49
50 auto solver = linear::Solver::create(solver_param["linear"], logger());
51 logger().debug("\tinternal solver {}", solver->name());
52
53 const int precond_num = problem_dim * n_bases;
54
55 StiffnessMatrix velocity_mass = velocity_mass1 / beta_dt;
56 // velocity_mass.setZero();
57
58 igl::Timer time;
59
60 time.start();
61 StiffnessMatrix stoke_stiffness;
62 Eigen::VectorXd prev_sol_mass(rhs.size()); // prev_sol_mass=prev_sol
63 prev_sol_mass.setZero();
64 prev_sol_mass.block(0, 0, velocity_mass.rows(), 1) = velocity_mass * prev_sol.block(0, 0, velocity_mass.rows(), 1);
65 for (int i : boundary_nodes)
66 prev_sol_mass[i] = 0;
67
68 AssemblerUtils::merge_mixed_matrices(n_bases, n_pressure_bases, problem_dim, use_avg_pressure,
69 velocity_stiffness + velocity_mass, mixed_stiffness, pressure_stiffness,
70 stoke_stiffness);
71 time.stop();
72 stokes_matrix_time = time.getElapsedTimeInSec();
73 logger().debug("\tStokes matrix assembly time {}s", time.getElapsedTimeInSec());
74
75 time.start();
76
77 Eigen::VectorXd b = rhs + prev_sol_mass;
78
79 if (use_avg_pressure)
80 {
81 b[b.size() - 1] = 0;
82 }
83 dirichlet_solve(*solver, stoke_stiffness, b, boundary_nodes, x, precond_num, "", false, true, use_avg_pressure);
84 // solver->get_info(solver_info);
85 time.stop();
86 stokes_solve_time = time.getElapsedTimeInSec();
87 logger().debug("\tStokes solve time {}s", time.getElapsedTimeInSec());
88 logger().debug("\tStokes solver error: {}", (stoke_stiffness * x - b).norm());
89 // return;
90
91 std::vector<bool> zero_col(stoke_stiffness.cols(), true);
92 for (int k = 0; k < stoke_stiffness.outerSize(); ++k)
93 {
94 for (StiffnessMatrix::InnerIterator it(stoke_stiffness, k); it; ++it)
95 {
96 if (fabs(it.value()) > 1e-12)
97 zero_col[it.col()] = false;
98 }
99 }
100 std::vector<int> skipping;
101 for (int i = 0; i < zero_col.size(); ++i)
102 {
103 if (zero_col[i])
104 {
105 skipping.push_back(i);
106 }
107 }
108
109 assembly_time = 0;
110 inverting_time = 0;
111
112 int it = 0;
113 double nlres_norm = 0;
114 b = rhs + prev_sol_mass;
115
116 if (use_avg_pressure)
117 {
118 b[b.size() - 1] = 0;
119 }
120 it += minimize_aux(true, skipping,
121 n_bases,
122 n_pressure_bases,
123 t,
124 bases,
125 gbases,
126 velocity_assembler,
127 ass_vals_cache,
128 boundary_nodes,
129 use_avg_pressure,
130 problem_dim,
131 is_volume,
132 velocity_stiffness, mixed_stiffness, pressure_stiffness, velocity_mass, b, 1e-3, solver, nlres_norm, x);
133 it += minimize_aux(false, skipping,
134 n_bases,
135 n_pressure_bases,
136 t,
137 bases,
138 gbases,
139 velocity_assembler,
140 ass_vals_cache,
141 boundary_nodes,
142 use_avg_pressure,
143 problem_dim,
144 is_volume,
145 velocity_stiffness, mixed_stiffness, pressure_stiffness, velocity_mass, b, gradNorm, solver, nlres_norm, x);
146
147 solver_info["iterations"] = it;
148 solver_info["gradNorm"] = nlres_norm;
149
150 assembly_time /= it;
151 inverting_time /= it;
152
153 solver_info["time_assembly"] = assembly_time;
154 solver_info["time_inverting"] = inverting_time;
155 solver_info["time_stokes_assembly"] = stokes_matrix_time;
156 solver_info["time_stokes_solve"] = stokes_solve_time;
157
158 logger().info("finished with niter: {}, ||g||_2 = {}", it, nlres_norm);
159 }
160
162 const bool is_picard,
163 const std::vector<int> &skipping,
164 const int n_bases,
165 const int n_pressure_bases,
166 const double t,
167 const std::vector<basis::ElementBases> &bases,
168 const std::vector<basis::ElementBases> &gbases,
169 assembler::NavierStokesVelocity &velocity_assembler,
170 const assembler::AssemblyValsCache &ass_vals_cache,
171 const std::vector<int> &boundary_nodes,
172 const bool use_avg_pressure,
173 const int problem_dim,
174 const bool is_volume,
175 const StiffnessMatrix &velocity_stiffness, const StiffnessMatrix &mixed_stiffness, const StiffnessMatrix &pressure_stiffness,
176 const StiffnessMatrix &velocity_mass,
177 const Eigen::VectorXd &rhs, const double grad_norm,
178 std::unique_ptr<linear::Solver> &solver, double &nlres_norm,
179 Eigen::VectorXd &x)
180 {
181 igl::Timer time;
182 const int precond_num = problem_dim * n_bases;
183
184 StiffnessMatrix nl_matrix;
185 StiffnessMatrix total_matrix;
186 SparseMatrixCache mat_cache;
187
188 time.start();
189 velocity_assembler.set_picard(true);
190 velocity_assembler.assemble_hessian(is_volume, n_bases, false, bases, gbases, ass_vals_cache, t, 0, x, Eigen::MatrixXd(), mat_cache, nl_matrix);
191 AssemblerUtils::merge_mixed_matrices(n_bases, n_pressure_bases, problem_dim, use_avg_pressure,
192 (velocity_stiffness + nl_matrix) + velocity_mass, mixed_stiffness, pressure_stiffness,
193 total_matrix);
194 time.stop();
195 assembly_time = time.getElapsedTimeInSec();
196 logger().debug("\tNavier Stokes assembly time {}s", time.getElapsedTimeInSec());
197
198 Eigen::VectorXd nlres = -(total_matrix * x) + rhs;
199 for (int i : boundary_nodes)
200 nlres[i] = 0;
201 for (int i : skipping)
202 nlres[i] = 0;
203 Eigen::VectorXd dx;
204 nlres_norm = nlres.norm();
205 logger().debug("\tInitial residula norm {}", nlres_norm);
206
207 int it = 0;
208
209 while (nlres_norm > grad_norm && it < iterations)
210 {
211 ++it;
212
213 time.start();
214 if (!is_picard)
215 {
216 velocity_assembler.set_picard(false);
217 velocity_assembler.assemble_hessian(is_volume, n_bases, false, bases, gbases, ass_vals_cache, t, 0, x, Eigen::MatrixXd(), mat_cache, nl_matrix);
218 AssemblerUtils::merge_mixed_matrices(n_bases, n_pressure_bases, problem_dim, use_avg_pressure,
219 (velocity_stiffness + nl_matrix) + velocity_mass, mixed_stiffness, pressure_stiffness,
220 total_matrix);
221 }
222 dirichlet_solve(*solver, total_matrix, nlres, boundary_nodes, dx, precond_num, "", false, true, use_avg_pressure);
223 // for (int i : boundary_nodes)
224 // dx[i] = 0;
225 time.stop();
226 inverting_time += time.getElapsedTimeInSec();
227 logger().debug("\tinverting time {}s", time.getElapsedTimeInSec());
228
229 x += dx;
230 // TODO check for nans
231
232 time.start();
233 velocity_assembler.set_picard(true);
234 velocity_assembler.assemble_hessian(is_volume, n_bases, false, bases, gbases, ass_vals_cache, t, 0, x, Eigen::MatrixXd(), mat_cache, nl_matrix);
235 AssemblerUtils::merge_mixed_matrices(n_bases, n_pressure_bases, problem_dim, use_avg_pressure,
236 (velocity_stiffness + nl_matrix) + velocity_mass, mixed_stiffness, pressure_stiffness,
237 total_matrix);
238 time.stop();
239 logger().debug("\tassembly time {}s", time.getElapsedTimeInSec());
240 assembly_time += time.getElapsedTimeInSec();
241
242 nlres = -(total_matrix * x) + rhs;
243 for (int i : boundary_nodes)
244 nlres[i] = 0;
245 for (int i : skipping)
246 nlres[i] = 0;
247
248 // if (use_avg_pressure)
249 // nlres[nlres.size() - 1] = 0;
250 nlres_norm = nlres.norm();
251
252 logger().debug("\titer: {}, ||g||_2 = {}, ||step|| = {}\n",
253 it, nlres_norm, dx.norm());
254 }
255
256 if (it >= iterations)
257 {
258 log_and_throw_error("Reaching the max number of iterations!");
259 }
260
261 // solver_info["internal_solver"] = internal_solver;
262 // solver_info["internal_solver_first"] = internal_solver.front();
263 // solver_info["status"] = this->status();
264
265 return it;
266 }
267 } // namespace solver
268} // namespace polyfem
int x
static void merge_mixed_matrices(const int n_bases, const int n_pressure_bases, const int problem_dim, const bool add_average, const StiffnessMatrix &velocity_stiffness, const StiffnessMatrix &mixed_stiffness, const StiffnessMatrix &pressure_stiffness, StiffnessMatrix &stiffness)
utility to merge 3 blocks of mixed matrices, A=velocity_stiffness, B=mixed_stiffness,...
Caches basis evaluation and geometric mapping at every element.
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
std::string name() const override
void minimize(const int n_bases, const int n_pressure_bases, const double t, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, assembler::NavierStokesVelocity &velocity_assembler, const assembler::AssemblyValsCache &ass_vals_cache, const std::vector< int > &boundary_nodes, const bool use_avg_pressure, const int problem_dim, const bool is_volume, const double beta_dt, const Eigen::VectorXd &prev_sol, const StiffnessMatrix &velocity_stiffness, const StiffnessMatrix &mixed_stiffness, const StiffnessMatrix &pressure_stiffness, const StiffnessMatrix &velocity_mass1, const Eigen::MatrixXd &rhs, Eigen::VectorXd &x)
int minimize_aux(const bool is_picard, const std::vector< int > &skipping, const int n_bases, const int n_pressure_bases, const double t, const std::vector< basis::ElementBases > &bases, const std::vector< basis::ElementBases > &gbases, assembler::NavierStokesVelocity &velocity_assembler, const assembler::AssemblyValsCache &ass_vals_cache, const std::vector< int > &boundary_nodes, const bool use_avg_pressure, const int problem_dim, const bool is_volume, const StiffnessMatrix &velocity_stiffness, const StiffnessMatrix &mixed_stiffness, const StiffnessMatrix &pressure_stiffness, const StiffnessMatrix &velocity_mass, const Eigen::VectorXd &rhs, const double grad_norm, std::unique_ptr< polysolve::linear::Solver > &solver, double &nlres_norm, Eigen::VectorXd &x)
spdlog::logger & logger()
Retrieves the current logger.
Definition Logger.cpp:42
nlohmann::json json
Definition Common.hpp:9
void log_and_throw_error(const std::string &msg)
Definition Logger.cpp:71
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix
Definition Types.hpp:22