39 namespace PressureAdvectionDiffusionValidation
52 wind[0] = sin(6.0 * x[1]);
53 wind[1] = cos(6.0 * x[0]);
69 u[2] = x[0] * x[0] * pow(1.0 - x[0], 2.0) * x[1] * x[1] *
74 u[2] = 0.1E1 -
Peclet * x[0] * (0.1E1 - 0.5 * x[0]);
83 u = x[0] * x[0] * pow(1.0 - x[0], 2.0) * x[1] * x[1] *
88 u = 0.1E1 -
Peclet * x[0] * (0.1E1 - 0.5 * x[0]);
106 (sin(0.6E1 * x[1]) * (2.0 * x[0] * pow(1.0 - x[0], 2.0) * x[1] *
107 x[1] * pow(1.0 - x[1], 2.0) -
108 2.0 * x[0] * x[0] * (1.0 - x[0]) * x[1] *
109 x[1] * pow(1.0 - x[1], 2.0)) +
110 cos(0.6E1 * x[0]) * (2.0 * x[0] * x[0] * pow(1.0 - x[0], 2.0) *
111 x[1] * pow(1.0 - x[1], 2.0) -
112 2.0 * x[0] * x[0] * pow(1.0 - x[0], 2.0) *
113 x[1] * x[1] * (1.0 - x[1]))) -
114 2.0 * pow(1.0 - x[0], 2.0) * x[1] * x[1] * pow(1.0 - x[1], 2.0) +
115 8.0 * x[0] * (1.0 - x[0]) * x[1] * x[1] * pow(1.0 - x[1], 2.0) -
116 2.0 * x[0] * x[0] * x[1] * x[1] * pow(1.0 - x[1], 2.0) -
117 2.0 * x[0] * x[0] * pow(1.0 - x[0], 2.0) * pow(1.0 - x[1], 2.0) +
118 8.0 * x[0] * x[0] * pow(1.0 - x[0], 2.0) * x[1] * (1.0 - x[1]) -
119 2.0 * x[0] * x[0] * pow(1.0 - x[0], 2.0) * x[1] * x[1];
123 source =
Peclet * (-0.1E1 *
Peclet * (0.1E1 - 0.5 * x[0]) +
150 bool doc_block_matrices =
false;
153 bool raytime_flag =
false;
165 double clean_up_memory_time =
166 t_clean_up_memory_end - t_clean_up_memory_start;
169 oomph_info <<
"LSC: clean_up_memory_time: " << clean_up_memory_time
178 std::ostringstream error_message;
179 error_message <<
"The navier stokes elements mesh pointer must be set.\n"
180 <<
"Use method set_navier_stokes_mesh(...)";
182 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
202 if (cr_matrix_pt == 0)
204 std::ostringstream error_message;
206 <<
"NavierStokesSchurComplementPreconditioner only works with "
207 <<
"CRDoubleMatrix matrices" << std::endl;
209 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
214 if (doc_block_matrices)
216 std::stringstream junk;
217 junk <<
"j_matrix" <<
comm_pt()->my_rank() <<
".dat";
218 oomph_info <<
"About to output " << junk.str() << std::endl;
220 oomph_info <<
"Done output of " << junk.str() << std::endl;
251 double block_setup_time = t_block_setup_finish - t_block_setup_start;
254 oomph_info <<
"Time for block_setup(...) [sec]: " << block_setup_time
260 oomph_info <<
"LSC: block_setup: " << block_setup_time << std::endl;
269 if (F_block_preconditioner_pt == 0)
280 double get_B_time = t_get_B_finish - t_get_B_start;
283 oomph_info <<
"Time to get B [sec]: " << get_B_time <<
"\n";
288 oomph_info <<
"LSC: get block B get_B_time: " << get_B_time << std::endl;
291 if (doc_block_matrices)
293 std::stringstream junk;
294 junk <<
"b_matrix" <<
comm_pt()->my_rank() <<
".dat";
296 oomph_info <<
"Done output of " << junk.str() << std::endl;
309 inv_p_mass_pt, inv_v_mass_pt,
false);
315 inv_p_mass_pt, inv_v_mass_pt,
true);
320 double ivmm_assembly_time = ivmm_assembly_finish_t - ivmm_assembly_start_t;
323 oomph_info <<
"Time to assemble inverse diagonal velocity and pressure"
324 <<
"mass matrices) [sec]: " << ivmm_assembly_time <<
"\n";
328 oomph_info <<
"LSC: ivmm_assembly_time: " << ivmm_assembly_time
333 if (doc_block_matrices)
335 std::stringstream junk;
336 junk <<
"inv_v_mass_matrix" <<
comm_pt()->my_rank() <<
".dat";
338 oomph_info <<
"Done output of " << junk.str() << std::endl;
348 double t_get_Bt_time = t_get_Bt_finish - t_get_Bt_start;
351 oomph_info <<
"Time to get Bt [sec]: " << t_get_Bt_time << std::endl;
355 oomph_info <<
"LSC: get block Bt: " << t_get_Bt_time << std::endl;
358 if (doc_block_matrices)
360 std::stringstream junk;
361 junk <<
"bt_matrix" <<
comm_pt()->my_rank() <<
".dat";
363 oomph_info <<
"Done output of " << junk.str() << std::endl;
373 inv_v_mass_pt->
multiply(*bt_pt, *qbt_pt);
381 double t_QBt_time = t_QBt_matrix_finish - t_QBt_matrix_start;
384 oomph_info <<
"Time to generate QBt [sec]: " << t_QBt_time << std::endl;
386 delete inv_v_mass_pt;
390 oomph_info <<
"LSC: t_QBt_time (matrix multiplicaton): " << t_QBt_time
397 b_pt->
multiply(*bt_pt, *p_matrix_pt);
400 double t_p_time = t_p_matrix_finish - t_p_matrix_start;
403 oomph_info <<
"Time to generate P [sec]: " << t_p_time << std::endl;
411 oomph_info <<
"LSC: t_p_time (matrix multiplication): " << t_p_time
422 double t_p_time2 = t_QBt_MV_finish - t_QBt_MV_start;
425 oomph_info <<
"Time to build QBt matrix vector operator [sec]: "
426 << t_p_time2 << std::endl;
436 oomph_info <<
"LSC: QBt (setup MV product): " << t_p_time2 << std::endl;
454 double t_get_Fp_time = t_get_Fp_finish - t_get_Fp_start;
455 oomph_info <<
"Time to get Fp [sec]: " << t_get_Fp_time << std::endl;
461 fp_matrix_pt->
multiply(*inv_p_mass_pt, *fp_qp_inv_pt);
470 double t_p_time = t_Fp_Qp_inv_MV_finish - t_Fp_Qp_inv_MV_start;
471 oomph_info <<
"Time to build Fp Qp^{-1} matrix vector operator [sec]: "
472 << t_p_time << std::endl;
475 delete inv_p_mass_pt;
488 double t_get_F_time = t_get_F_finish - t_get_F_start;
491 oomph_info <<
"Time to get F [sec]: " << t_get_F_time << std::endl;
495 oomph_info <<
"LSC: get_block t_get_F_time: " << t_get_F_time
505 double t_F_MV_time = t_F_MV_finish - t_F_MV_start;
508 oomph_info <<
"Time to build F Matrix Vector Operator [sec]: "
509 << t_F_MV_time << std::endl;
513 oomph_info <<
"LSC: MV product setup t_F_MV_time: " << t_F_MV_time
531 double t_get_Bt_time2 = t_get_Bt_finish - t_get_Bt_start;
534 oomph_info <<
"Time to get Bt [sec]: " << t_get_Bt_time2 << std::endl;
538 oomph_info <<
"LSC: get_block t_get_Bt_time2: " << t_get_Bt_time2
559 double t_Bt_MV_time = t_Bt_MV_finish - t_Bt_MV_start;
562 oomph_info <<
"LSC: MV product setup t_Bt_MV_time: " << t_Bt_MV_time
576 if (doc_block_matrices)
578 std::stringstream junk;
579 junk <<
"p_matrix" <<
comm_pt()->my_rank() <<
".dat";
581 oomph_info <<
"Done output of " << junk.str() << std::endl;
589 double t_p_prec_time = t_p_prec_finish - t_p_prec_start;
592 oomph_info <<
"P sub-preconditioner setup time [sec]: " << t_p_prec_time
597 oomph_info <<
"LSC: p_prec setup time: " << t_p_prec_time << std::endl;
615 unsigned nvelocity_dof_types =
619 for (
unsigned i = 0;
i < nvelocity_dof_types;
i++)
637 double t_f_prec_time = t_f_prec_finish - t_f_prec_start;
640 oomph_info <<
"F sub-preconditioner setup time [sec]: " << t_f_prec_time
645 oomph_info <<
"LSC: f_prec setup time: " << t_f_prec_time << std::endl;
664 std::ostringstream error_message;
665 error_message <<
"setup must be called before using preconditioner_solve";
667 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
673 std::ostringstream error_message;
674 error_message <<
"The vectors z and r must have the same number of "
677 OOMPH_CURRENT_FUNCTION,
678 OOMPH_EXCEPTION_LOCATION);
712 std::ostringstream error_message;
713 error_message <<
"P_preconditioner_pt has not been set.";
715 OOMPH_CURRENT_FUNCTION,
716 OOMPH_EXCEPTION_LOCATION);
729 another_temp_vec.
clear();
738 another_temp_vec.
clear();
760 std::ostringstream error_message;
761 error_message <<
"P_preconditioner_pt has not been set.";
763 OOMPH_CURRENT_FUNCTION,
764 OOMPH_EXCEPTION_LOCATION);
769 another_temp_vec.
clear();
782 temp_vec -= another_temp_vec;
800 another_temp_vec.
clear();
805 another_temp_vec += temp_vec;
814 std::ostringstream error_message;
815 error_message <<
"F_preconditioner_pt has not been set.";
817 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
857 double* v_values =
new double[v_nrow_local];
858 for (
unsigned i = 0;
i < v_nrow_local;
i++)
865 unsigned p_first_row = 0;
866 unsigned p_nrow_local = 0;
868 double* p_values = 0;
877 p_values =
new double[p_nrow_local];
878 for (
unsigned i = 0;
i < p_nrow_local;
i++)
900 unsigned nproc =
comm_pt()->nproc();
903 unsigned my_rank =
comm_pt()->my_rank();
910 unsigned first_lookup_row = 0;
911 unsigned last_lookup_row = 0;
924 unsigned max_block = 0;
926 for (
unsigned block_index = 0; block_index <= max_block; block_index++)
929 unsigned v_or_p_first_row = v_first_row;
930 double* v_or_p_values = v_values;
932 if (block_index == 1)
934 v_or_p_first_row = p_first_row;
935 v_or_p_values = p_values;
961 for (
unsigned e = 0;
e < n_el;
e++)
972 unsigned el_dof = el_pt->
ndof();
979 unsigned which_one = 2;
980 if (block_index == 1) which_one = 1;
988 el_pmm_diagonal, el_vmm_diagonal, which_one);
992 for (
unsigned i = 0;
i < el_dof;
i++)
998 if ((eqn_number >= first_lookup_row) &&
999 (eqn_number <= last_lookup_row))
1002 if (this->
block_number(eqn_number) == int(block_index))
1008 for (
unsigned p = 0; p < nproc; p++)
1010 if ((index >= velocity_or_press_dist_pt->
first_row(p)) &&
1011 (index < (velocity_or_press_dist_pt->
first_row(p) +
1018 if (block_index == 0)
1020 v_or_p_values[index - v_or_p_first_row] +=
1023 else if (block_index == 1)
1025 v_or_p_values[index - v_or_p_first_row] +=
1032 if (block_index == 0)
1034 classified_contributions_send[p].push_back(
1035 el_vmm_diagonal[
i]);
1036 classified_indices_send[p].push_back(index);
1038 else if (block_index == 1)
1040 classified_contributions_send[p].push_back(
1041 el_pmm_diagonal[
i]);
1042 classified_indices_send[p].push_back(index);
1068 if (block_index == 0)
1070 unclassified_contributions_send[p].push_back(
1071 el_vmm_diagonal[
i]);
1072 unclassified_indices_send[p].push_back(eqn_number);
1074 else if (block_index == 1)
1076 unclassified_contributions_send[p].push_back(
1077 el_pmm_diagonal[
i]);
1078 unclassified_indices_send[p].push_back(eqn_number);
1090 unsigned* n_unclassified_send =
new unsigned[nproc];
1091 for (
unsigned p = 0; p < nproc; p++)
1095 n_unclassified_send[p] = 0;
1099 n_unclassified_send[p] = unclassified_contributions_send[p].size();
1104 unsigned* n_unclassified_recv =
new unsigned[nproc];
1105 MPI_Alltoall(n_unclassified_send,
1108 n_unclassified_recv,
1114 MPI_Aint base_displacement;
1115 MPI_Get_address(v_or_p_values, &base_displacement);
1124 for (
unsigned p = 0; p < nproc; p++)
1129 if (n_unclassified_recv[p] > 0)
1131 unclassified_contributions_recv[p] =
1132 new double[n_unclassified_recv[p]];
1133 unclassified_indices_recv[p] =
1134 new unsigned[n_unclassified_recv[p]];
1137 MPI_Datatype recv_types[2];
1138 MPI_Aint recv_displacements[2];
1142 MPI_Type_contiguous(
1143 n_unclassified_recv[p], MPI_DOUBLE, &recv_types[0]);
1144 MPI_Type_commit(&recv_types[0]);
1145 MPI_Get_address(unclassified_contributions_recv[p],
1146 &recv_displacements[0]);
1147 recv_displacements[0] -= base_displacement;
1151 MPI_Type_contiguous(
1152 n_unclassified_recv[p], MPI_UNSIGNED, &recv_types[1]);
1153 MPI_Type_commit(&recv_types[1]);
1154 MPI_Get_address(unclassified_indices_recv[p],
1155 &recv_displacements[1]);
1156 recv_displacements[1] -= base_displacement;
1160 MPI_Datatype final_recv_type;
1161 MPI_Type_create_struct(
1162 2, recv_sz, recv_displacements, recv_types, &final_recv_type);
1163 MPI_Type_commit(&final_recv_type);
1167 MPI_Irecv(v_or_p_values,
1174 unclassified_recv_requests.push_back(req);
1175 unclassified_recv_proc.push_back(p);
1176 MPI_Type_free(&recv_types[0]);
1177 MPI_Type_free(&recv_types[1]);
1178 MPI_Type_free(&final_recv_type);
1182 if (n_unclassified_send[p] > 0)
1185 MPI_Datatype send_types[2];
1186 MPI_Aint send_displacements[2];
1190 MPI_Type_contiguous(
1191 n_unclassified_send[p], MPI_DOUBLE, &send_types[0]);
1192 MPI_Type_commit(&send_types[0]);
1193 MPI_Get_address(&unclassified_contributions_send[p][0],
1194 &send_displacements[0]);
1195 send_displacements[0] -= base_displacement;
1199 MPI_Type_contiguous(
1200 n_unclassified_send[p], MPI_UNSIGNED, &send_types[1]);
1201 MPI_Type_commit(&send_types[1]);
1202 MPI_Get_address(&unclassified_indices_send[p][0],
1203 &send_displacements[1]);
1204 send_displacements[1] -= base_displacement;
1208 MPI_Datatype final_send_type;
1209 MPI_Type_create_struct(
1210 2, send_sz, send_displacements, send_types, &final_send_type);
1211 MPI_Type_commit(&final_send_type);
1215 MPI_Isend(v_or_p_values,
1222 unclassified_send_requests.push_back(req);
1223 MPI_Type_free(&send_types[0]);
1224 MPI_Type_free(&send_types[1]);
1225 MPI_Type_free(&final_send_type);
1231 unsigned n_unclassified_recv_req = unclassified_recv_requests.size();
1232 while (n_unclassified_recv_req > 0)
1237 MPI_Waitany(n_unclassified_recv_req,
1238 &unclassified_recv_requests[0],
1241 unsigned p = unclassified_recv_proc[req_num];
1242 unclassified_recv_requests.erase(unclassified_recv_requests.begin() +
1244 unclassified_recv_proc.erase(unclassified_recv_proc.begin() +
1246 n_unclassified_recv_req--;
1250 unsigned n_recv = n_unclassified_recv[p];
1251 for (
unsigned i = 0;
i < n_recv;
i++)
1253 unsigned eqn_number = unclassified_indices_recv[p][
i];
1255 if (this->
block_number(eqn_number) ==
int(block_index))
1261 for (
unsigned pp = 0; pp < nproc; pp++)
1263 if ((index >= velocity_or_press_dist_pt->
first_row(pp)) &&
1264 (index < (velocity_or_press_dist_pt->
first_row(pp) +
1271 v_or_p_values[index - v_or_p_first_row] +=
1272 unclassified_contributions_recv[p][
i];
1277 double v = unclassified_contributions_recv[p][
i];
1278 classified_contributions_send[pp].push_back(v);
1279 classified_indices_send[pp].push_back(index);
1287 delete[] unclassified_contributions_recv[p];
1288 delete[] unclassified_indices_recv[p];
1290 delete[] n_unclassified_recv;
1299 unsigned* n_classified_send =
new unsigned[nproc];
1300 for (
unsigned p = 0; p < nproc; p++)
1304 n_classified_send[p] = 0;
1308 n_classified_send[p] = classified_contributions_send[p].size();
1313 unsigned* n_classified_recv =
new unsigned[nproc];
1314 MPI_Alltoall(n_classified_send,
1329 for (
unsigned p = 0; p < nproc; p++)
1334 if (n_classified_recv[p] > 0)
1336 classified_contributions_recv[p] =
1337 new double[n_classified_recv[p]];
1338 classified_indices_recv[p] =
new unsigned[n_classified_recv[p]];
1341 MPI_Datatype recv_types[2];
1342 MPI_Aint recv_displacements[2];
1346 MPI_Type_contiguous(
1347 n_classified_recv[p], MPI_DOUBLE, &recv_types[0]);
1348 MPI_Type_commit(&recv_types[0]);
1349 MPI_Get_address(classified_contributions_recv[p],
1350 &recv_displacements[0]);
1351 recv_displacements[0] -= base_displacement;
1355 MPI_Type_contiguous(
1356 n_classified_recv[p], MPI_UNSIGNED, &recv_types[1]);
1357 MPI_Type_commit(&recv_types[1]);
1358 MPI_Get_address(classified_indices_recv[p],
1359 &recv_displacements[1]);
1360 recv_displacements[1] -= base_displacement;
1364 MPI_Datatype final_recv_type;
1365 MPI_Type_create_struct(
1366 2, recv_sz, recv_displacements, recv_types, &final_recv_type);
1367 MPI_Type_commit(&final_recv_type);
1371 MPI_Irecv(v_or_p_values,
1378 classified_recv_requests.push_back(req);
1379 classified_recv_proc.push_back(p);
1380 MPI_Type_free(&recv_types[0]);
1381 MPI_Type_free(&recv_types[1]);
1382 MPI_Type_free(&final_recv_type);
1386 if (n_classified_send[p] > 0)
1389 MPI_Datatype send_types[2];
1390 MPI_Aint send_displacements[2];
1394 MPI_Type_contiguous(
1395 n_classified_send[p], MPI_DOUBLE, &send_types[0]);
1396 MPI_Type_commit(&send_types[0]);
1397 MPI_Get_address(&classified_contributions_send[p][0],
1398 &send_displacements[0]);
1399 send_displacements[0] -= base_displacement;
1403 MPI_Type_contiguous(
1404 n_classified_send[p], MPI_UNSIGNED, &send_types[1]);
1405 MPI_Type_commit(&send_types[1]);
1406 MPI_Get_address(&classified_indices_send[p][0],
1407 &send_displacements[1]);
1408 send_displacements[1] -= base_displacement;
1412 MPI_Datatype final_send_type;
1413 MPI_Type_create_struct(
1414 2, send_sz, send_displacements, send_types, &final_send_type);
1415 MPI_Type_commit(&final_send_type);
1419 MPI_Isend(v_or_p_values,
1426 classified_send_requests.push_back(req);
1427 MPI_Type_free(&send_types[0]);
1428 MPI_Type_free(&send_types[1]);
1429 MPI_Type_free(&final_send_type);
1435 unsigned n_classified_recv_req = classified_recv_requests.size();
1436 while (n_classified_recv_req > 0)
1441 MPI_Waitany(n_classified_recv_req,
1442 &classified_recv_requests[0],
1445 unsigned p = classified_recv_proc[req_num];
1446 classified_recv_requests.erase(classified_recv_requests.begin() +
1448 classified_recv_proc.erase(classified_recv_proc.begin() + req_num);
1449 n_classified_recv_req--;
1453 unsigned n_recv = n_classified_recv[p];
1454 for (
unsigned i = 0;
i < n_recv;
i++)
1456 v_or_p_values[classified_indices_recv[p][
i] - v_or_p_first_row] +=
1457 classified_contributions_recv[p][
i];
1461 delete[] classified_contributions_recv[p];
1462 delete[] classified_indices_recv[p];
1466 unsigned n_unclassified_send_req = unclassified_send_requests.size();
1467 if (n_unclassified_send_req > 0)
1469 MPI_Waitall(n_unclassified_send_req,
1470 &unclassified_send_requests[0],
1473 delete[] unclassified_contributions_send;
1474 delete[] unclassified_indices_send;
1475 delete[] n_unclassified_send;
1478 unsigned n_classified_send_req = classified_send_requests.size();
1479 if (n_classified_send_req > 0)
1481 MPI_Waitall(n_classified_send_req,
1482 &classified_send_requests[0],
1485 delete[] classified_indices_send;
1486 delete[] classified_contributions_send;
1487 delete[] n_classified_recv;
1488 delete[] n_classified_send;
1491 if (block_index == 0)
1493 v_values = v_or_p_values;
1495 else if (block_index == 1)
1497 p_values = v_or_p_values;
1510 unsigned which_one = 0;
1514 for (
unsigned e = 0;
e < n_el;
e++)
1522 unsigned el_dof = el_pt->
ndof();
1532 if (cast_el_pt != 0)
1535 el_pmm_diagonal, el_vmm_diagonal, which_one);
1540 std::ostringstream error_message;
1542 <<
"Navier-Stokes mesh contains element that is not of type \n"
1543 <<
"NavierStokesElementWithDiagonalMassMatrices. \n"
1544 <<
"The element is in fact of type " <<
typeid(*el_pt).name()
1545 <<
"\nWe'll assume that it does not make a used contribution \n"
1546 <<
"to the inverse diagonal mass matrix used in the "
1548 <<
"and (to avoid divisions by zero) fill in dummy unit entries.\n"
1549 <<
"[This case currently arises with flux control problems\n"
1550 <<
"where for simplicity the NetFluxControlElement has been added "
1552 <<
"to the Navier Stokes mesh -- this should probably be changed "
1554 <<
"some point -- if you get this warning in any other context\n"
1555 <<
"you should check your code very carefully]\n";
1557 "NavierStokesSchurComplementPreconditioner::assemble_"
1558 "inv_press_and_veloc_mass_matrix_diagonal()",
1559 OOMPH_EXCEPTION_LOCATION);
1563 for (
unsigned j = 0; j < el_dof; j++)
1565 el_vmm_diagonal[j] = 1.0;
1566 el_pmm_diagonal[j] = 1.0;
1571 for (
unsigned i = 0;
i < el_dof;
i++)
1583 if ((index >= v_first_row) &&
1584 (index < (v_first_row + v_nrow_local)))
1586 v_values[index - v_first_row] += el_vmm_diagonal[
i];
1598 if ((index >= p_first_row) &&
1599 (index < (p_first_row + p_nrow_local)))
1601 p_values[index - p_first_row] += el_pmm_diagonal[
i];
1610 int* v_column_index =
new int[v_nrow_local];
1611 int* v_row_start =
new int[v_nrow_local + 1];
1612 for (
unsigned i = 0;
i < v_nrow_local;
i++)
1615 if (v_values[
i] == 0.0)
1617 std::ostringstream error_message;
1618 error_message <<
"Zero entry in diagonal of velocity mass matrix\n"
1619 <<
"Index: " <<
i << std::endl;
1621 OOMPH_CURRENT_FUNCTION,
1622 OOMPH_EXCEPTION_LOCATION);
1625 v_values[
i] = 1.0 / v_values[
i];
1626 v_column_index[
i] = v_first_row +
i;
1629 v_row_start[v_nrow_local] = v_nrow_local;
1634 v_nrow, v_nrow_local, v_values, v_column_index, v_row_start);
1640 int* p_column_index =
new int[p_nrow_local];
1641 int* p_row_start =
new int[p_nrow_local + 1];
1642 for (
unsigned i = 0;
i < p_nrow_local;
i++)
1645 if (p_values[
i] == 0.0)
1647 std::ostringstream error_message;
1648 error_message <<
"Zero entry in diagonal of pressure mass matrix\n"
1649 <<
"Index: " <<
i << std::endl;
1651 OOMPH_CURRENT_FUNCTION,
1652 OOMPH_EXCEPTION_LOCATION);
1655 p_values[
i] = 1.0 / p_values[
i];
1657 p_column_index[
i] = p_first_row +
i;
1660 p_row_start[p_nrow_local] = p_nrow_local;
1665 p_nrow, p_nrow_local, p_values, p_column_index, p_row_start);
void return_block_vector(const unsigned &n, const DoubleVector &b, DoubleVector &v) const
Takes the n-th block ordered vector, b, and copies its entries to the appropriate entries in the natu...
unsigned ndof_types_in_mesh(const unsigned &i) const
Return the number of DOF types in mesh i. WARNING: This should only be used by the upper-most master ...
const LinearAlgebraDistribution * master_distribution_pt() const
Access function to the distribution of the master preconditioner. If this preconditioner does not hav...
void get_block(const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
int index_in_block(const unsigned &i_dof) const
Given a global dof number, returns the index in the block it belongs to. This is the overall index,...
void get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
Takes the naturally ordered vector, v and returns the n-th block vector, b. Here n is the block numbe...
int block_number(const unsigned &i_dof) const
Return the block number corresponding to a global index i_dof.
bool is_subsidiary_block_preconditioner() const
Return true if this preconditioner is a subsidiary preconditioner.
void get_block_other_matrix(const unsigned &i, const unsigned &j, CRDoubleMatrix *in_matrix_pt, CRDoubleMatrix &output_matrix)
Get a block from a different matrix using the blocking scheme that has already been set up.
unsigned ndof_types() const
Return the total number of DOF types.
void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< MATRIX > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
Function to turn this preconditioner into a subsidiary preconditioner that operates within a bigger "...
void setup_matrix_vector_product(MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const Vector< unsigned > &block_col_indices)
Setup a matrix vector product. matvec_prod_pt is a pointer to the MatrixVectorProduct,...
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
CRDoubleMatrix * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*,...
const LinearAlgebraDistribution * block_distribution_pt(const unsigned &b) const
Access function to the block distributions (const version).
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Set the i-th mesh for this block preconditioner. Note: The method set_nmesh(...) must be called befor...
A class for compressed row matrices. This is a distributable object.
void sparse_indexed_output_with_offset(std::string filename)
Indexed output function to print a matrix to a file as i,j,a(i,j) for a(i,j)!=0 only....
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
void build_without_copy(const unsigned &ncol, const unsigned &nnz, double *value, int *column_index, int *row_start)
keeps the existing distribution and just matrix that is stored without copying the matrix data
bool distributed() const
distribution is serial or distributed
unsigned nrow() const
access function to the number of global rows.
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
A vector in the mathematical sense, initially developed for linear algebra type applications....
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void clear()
wipes the DoubleVector
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
A Generalised Element class.
bool is_halo() const
Is this element a halo?
unsigned ndof() const
Return the number of equations/dofs in the element.
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
void multiply_transpose(const DoubleVector &x, DoubleVector &y) const
Apply the transpose of the operator to the vector x and return the result in the vector y.
void multiply(const DoubleVector &x, DoubleVector &y) const
Apply the operator to the vector x and return the result in the vector y.
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
unsigned long nelement() const
Return number of elements in the mesh.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)=0
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed,...
void setup()
Setup the preconditioner.
bool Doc_time
Set Doc_time to true for outputting results of timings.
Preconditioner * F_preconditioner_pt
Pointer to the 'preconditioner' for the F matrix.
bool Preconditioner_has_been_setup
Control flag is true if the preconditioner has been setup (used so we can wipe the data when the prec...
void clean_up_memory()
Helper function to delete preconditioner data.
Mesh * Navier_stokes_mesh_pt
the pointer to the mesh of block preconditionable Navier Stokes elements.
MatrixVectorProduct * F_mat_vec_pt
MatrixVectorProduct operator for F.
MatrixVectorProduct * E_mat_vec_pt
MatrixVectorProduct operator for E = Fp Qp^{-1} (only for Fp variant)
bool Use_LSC
Boolean to indicate use of LSC (true) or Fp (false) variant.
bool Using_default_p_preconditioner
flag indicating whether the default P preconditioner is used
MatrixVectorProduct * QBt_mat_vec_pt
MatrixVectorProduct operator for Qv^{-1} Bt.
bool Using_default_f_preconditioner
flag indicating whether the default F preconditioner is used
bool Allow_multiple_element_type_in_navier_stokes_mesh
Flag to indicate if there are multiple element types in the Navier-Stokes mesh.
Preconditioner * P_preconditioner_pt
Pointer to the 'preconditioner' for the pressure matrix.
Problem * problem_pt() const
void assemble_inv_press_and_veloc_mass_matrix_diagonal(CRDoubleMatrix *&inv_p_mass_pt, CRDoubleMatrix *&inv_v_mass_pt, const bool &do_both)
Helper function to assemble the inverse diagonals of the pressure and velocity mass matrices from the...
void get_pressure_advection_diffusion_matrix(CRDoubleMatrix &fp_matrix)
Get the pressure advection diffusion matrix.
bool F_preconditioner_is_block_preconditioner
Boolean indicating whether the momentum system preconditioner is a block preconditioner.
MatrixVectorProduct * Bt_mat_vec_pt
MatrixVectorProduct operator for Bt.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to Vector r.
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
void setup(DoubleMatrixBase *matrix_pt)
Setup the preconditioner: store the matrix pointer and the communicator pointer then call preconditio...
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
virtual void preconditioner_solve(const DoubleVector &r, DoubleVector &z)=0
Apply the preconditioner. Pure virtual generic interface function. This method should apply the preco...
An interface to allow SuperLU to be used as an (exact) Preconditioner.
void get_exact_u(const Vector< double > &x, Vector< double > &u)
Exact solution as a Vector.
double source_function(const Vector< double > &x_vect)
Source function required to make the solution above an exact solution.
void wind_function(const Vector< double > &x, Vector< double > &wind)
Wind.
unsigned Flag
Flag for solution.
double Peclet
Peclet number – overwrite with actual Reynolds number.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...