33 const int n_pressure_bases,
35 const std::vector<basis::ElementBases> &bases,
36 const std::vector<basis::ElementBases> &gbases,
39 const std::vector<int> &boundary_nodes,
40 const bool use_avg_pressure,
41 const int problem_dim,
43 const double beta_dt,
const Eigen::VectorXd &prev_sol,
46 const Eigen::MatrixXd &rhs, Eigen::VectorXd &
x)
48 assert(velocity_assembler.
name() ==
"NavierStokes");
51 logger().debug(
"\tinternal solver {}", solver->name());
53 const int precond_num = problem_dim * n_bases;
62 Eigen::VectorXd prev_sol_mass(rhs.size());
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)
69 velocity_stiffness + velocity_mass, mixed_stiffness, pressure_stiffness,
73 logger().debug(
"\tStokes matrix assembly time {}s", time.getElapsedTimeInSec());
77 Eigen::VectorXd b = rhs + prev_sol_mass;
83 dirichlet_solve(*solver, stoke_stiffness, b, boundary_nodes,
x, precond_num,
"",
false,
true, use_avg_pressure);
87 logger().debug(
"\tStokes solve time {}s", time.getElapsedTimeInSec());
88 logger().debug(
"\tStokes solver error: {}", (stoke_stiffness *
x - b).norm());
91 std::vector<bool> zero_col(stoke_stiffness.cols(),
true);
92 for (
int k = 0; k < stoke_stiffness.outerSize(); ++k)
94 for (StiffnessMatrix::InnerIterator it(stoke_stiffness, k); it; ++it)
96 if (fabs(it.value()) > 1e-12)
97 zero_col[it.col()] =
false;
100 std::vector<int> skipping;
101 for (
int i = 0; i < zero_col.size(); ++i)
105 skipping.push_back(i);
113 double nlres_norm = 0;
114 b = rhs + prev_sol_mass;
116 if (use_avg_pressure)
132 velocity_stiffness, mixed_stiffness, pressure_stiffness, velocity_mass, b, 1e-3, solver, nlres_norm,
x);
145 velocity_stiffness, mixed_stiffness, pressure_stiffness, velocity_mass, b,
gradNorm, solver, nlres_norm,
x);
158 logger().info(
"finished with niter: {}, ||g||_2 = {}", it, nlres_norm);
162 const bool is_picard,
163 const std::vector<int> &skipping,
165 const int n_pressure_bases,
167 const std::vector<basis::ElementBases> &bases,
168 const std::vector<basis::ElementBases> &gbases,
171 const std::vector<int> &boundary_nodes,
172 const bool use_avg_pressure,
173 const int problem_dim,
174 const bool is_volume,
177 const Eigen::VectorXd &rhs,
const double grad_norm,
178 std::unique_ptr<linear::Solver> &solver,
double &nlres_norm,
182 const int precond_num = problem_dim * n_bases;
190 velocity_assembler.
assemble_hessian(is_volume, n_bases,
false, bases, gbases, ass_vals_cache, t, 0,
x, Eigen::MatrixXd(), mat_cache, nl_matrix);
192 (velocity_stiffness + nl_matrix) + velocity_mass, mixed_stiffness, pressure_stiffness,
196 logger().debug(
"\tNavier Stokes assembly time {}s", time.getElapsedTimeInSec());
198 Eigen::VectorXd nlres = -(total_matrix *
x) + rhs;
199 for (
int i : boundary_nodes)
201 for (
int i : skipping)
204 nlres_norm = nlres.norm();
205 logger().debug(
"\tInitial residula norm {}", nlres_norm);
209 while (nlres_norm > grad_norm && it <
iterations)
217 velocity_assembler.
assemble_hessian(is_volume, n_bases,
false, bases, gbases, ass_vals_cache, t, 0,
x, Eigen::MatrixXd(), mat_cache, nl_matrix);
219 (velocity_stiffness + nl_matrix) + velocity_mass, mixed_stiffness, pressure_stiffness,
222 dirichlet_solve(*solver, total_matrix, nlres, boundary_nodes, dx, precond_num,
"",
false,
true, use_avg_pressure);
227 logger().debug(
"\tinverting time {}s", time.getElapsedTimeInSec());
234 velocity_assembler.
assemble_hessian(is_volume, n_bases,
false, bases, gbases, ass_vals_cache, t, 0,
x, Eigen::MatrixXd(), mat_cache, nl_matrix);
236 (velocity_stiffness + nl_matrix) + velocity_mass, mixed_stiffness, pressure_stiffness,
239 logger().debug(
"\tassembly time {}s", time.getElapsedTimeInSec());
242 nlres = -(total_matrix *
x) + rhs;
243 for (
int i : boundary_nodes)
245 for (
int i : skipping)
250 nlres_norm = nlres.norm();
252 logger().debug(
"\titer: {}, ||g||_2 = {}, ||step|| = {}\n",
253 it, nlres_norm, dx.norm());
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,...
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)