PetscSolverFeti.cc 40.9 KB
Newer Older
1001

1002
    // === Start global assembly procedure. ===
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016

    MatAssemblyBegin(mat_b_b, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_b_b, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_primal_primal, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_primal_primal, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_b_primal, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_b_primal, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_primal_b, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_primal_b, MAT_FINAL_ASSEMBLY);
	  

1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
    MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_bound_bound, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_bound_bound, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_interior_bound, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_interior_bound, MAT_FINAL_ASSEMBLY);

    MatAssemblyBegin(mat_bound_interior, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(mat_bound_interior, MAT_FINAL_ASSEMBLY);


1030
    // === Create and fill PETSc's right hand side vectors. ===
1031

1032
1033
1034
    VecCreate(PETSC_COMM_WORLD, &f_b);
    VecSetSizes(f_b, nRankB * nComponents, nOverallB * nComponents);
    VecSetType(f_b, VECMPI);
1035

1036
1037
    VecCreate(PETSC_COMM_WORLD, &f_primal);
    VecSetSizes(f_primal, nRankPrimals * nComponents, 
1038
		nOverallPrimals * nComponents);
1039
    VecSetType(f_primal, VECMPI);
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
    
    for (int i = 0; i < nComponents; i++) {
      DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
      for (dofIt.reset(); !dofIt.end(); ++dofIt) {
	int index = dofIt.getDOFIndex();
	if (primals.count(index)) {
	  TEST_EXIT_DBG(globalPrimalIndex.count(index))
	    ("Should not happen!\n");

	  index = globalPrimalIndex[index] * nComponents + i;
	  double value = *dofIt;
1051
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
1052
1053
1054
1055
1056
1057
	} else {
	  TEST_EXIT_DBG(globalIndexB.count(index))
	    ("Should not happen!\n");

	  index = globalIndexB[index] * nComponents + i;
	  double value = *dofIt;
1058
	  VecSetValues(f_b, 1, &index, &value, ADD_VALUES);
1059
1060
1061
1062
	}      
      }
    }

1063
1064
    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);
1065

1066
1067
    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);
1068
1069


1070
    // === Create and fill PETSc matrix for Lagrange constraints. ===
1071

1072
    createMatLagrange();
1073
1074

    
1075
1076
1077
1078
1079
1080
    // === Create PETSc solver for the Schur complement on primal variables. ===
    
    createSchurPrimalKsp();


    // === Create PETSc solver for the FETI-DP operator. ===
1081
1082

    createFetiKsp();
1083
1084
1085
  }


1086
  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
1087
  {
1088
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
1089

1090
1091
1092
    // Create transpose of Lagrange matrix.
    Mat mat_lagrange_transpose;
    MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
1093
1094


1095
    // === Create nested matrix which will contain the overall FETI system. ===
1096

1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
    Mat A;
    Mat nestedA[3][3];
    nestedA[0][0] = mat_b_b;
    nestedA[0][1] = mat_b_primal;
    nestedA[0][2] = mat_lagrange_transpose;
    nestedA[1][0] = mat_primal_b;
    nestedA[1][1] = mat_primal_primal;
    nestedA[1][2] = PETSC_NULL;
    nestedA[2][0] = mat_lagrange;
    nestedA[2][1] = PETSC_NULL;
    nestedA[2][2] = PETSC_NULL;
1108

1109
    MatCreateNest(PETSC_COMM_WORLD, 3, PETSC_NULL, 3, PETSC_NULL, &(nestedA[0][0]), &A);
1110

1111
1112
1113
    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  
1114
1115


1116
1117
1118
    int nRankNest = (nRankB + nRankPrimals + nRankLagrange) * nComponents;
    int nOverallNest = (nOverallB + nOverallPrimals + nOverallLagrange) * nComponents;
    int rStartNest = (rStartB + rStartPrimals + rStartLagrange) * nComponents;
1119

1120
1121
    {
      // === Test some matrix sizes. ===
1122

1123
1124
1125
1126
1127
      int matRow, matCol;
      MatGetLocalSize(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == nRankNest)("Should not happen!\n");
      mpi::globalAdd(matRow);
      TEST_EXIT_DBG(matRow == nOverallNest)("Should not happen!\n");
1128

1129
1130
1131
      MatGetOwnershipRange(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == rStartNest)("Should not happen!\n");
    }
1132

1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
    // === Create rhs and solution vectors for the overall FETI system. ===

    Vec f;
    VecCreate(PETSC_COMM_WORLD, &f);
    VecSetSizes(f, nRankNest, nOverallNest);
    VecSetType(f, VECMPI);

    Vec b;
    VecDuplicate(f, &b);

    
    // === Fill rhs vector by coping the primal and non primal PETSc vectors. ===
1145

1146
1147
    PetscScalar *local_f_b;
    VecGetArray(f_b, &local_f_b);
1148

1149
1150
    PetscScalar *local_f_primal;
    VecGetArray(f_primal, &local_f_primal);
1151

1152
1153
1154
1155
1156
1157
1158
    {
      int size;
      VecGetLocalSize(f_b, &size);
      TEST_EXIT_DBG(size == nRankB * nComponents)("Should not happen!\n");
      VecGetLocalSize(f_primal, &size);
      TEST_EXIT_DBG(size == nRankPrimals * nComponents)("Should not happen!\n");
    }
1159

1160
1161
    PetscScalar *local_f;
    VecGetArray(f, &local_f);
1162

1163
1164
1165
1166
1167
    int index = 0;
    for (int i = 0; i < nRankB * nComponents; i++)
      local_f[index++] = local_f_b[i];
    for (int i = 0; i < nRankPrimals * nComponents; i++)
      local_f[index++] = local_f_primal[i];
1168

1169
1170
1171
    VecRestoreArray(f, &local_f);  
    VecRestoreArray(f_b, &local_f_b);
    VecRestoreArray(f_primal, &local_f_primal);
1172

1173
1174
    
    // === Create solver and solve the overall FETI system. ===
1175

1176
1177
1178
1179
    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, A, A, SAME_NONZERO_PATTERN);
    KSPSetFromOptions(ksp);
1180
1181


1182
    KSPSolve(ksp, f, b);
1183
1184


1185
1186
1187
1188
1189
    // === Reconstruct FETI solution vectors. ===
    
    Vec u_b, u_primal;
    VecDuplicate(f_b, &u_b);
    VecDuplicate(f_primal, &u_primal);
1190
1191
    

1192
1193
    PetscScalar *local_b;
    VecGetArray(b, &local_b);
1194

1195
1196
    PetscScalar *local_u_b;
    VecGetArray(u_b, &local_u_b);
1197

1198
1199
    PetscScalar *local_u_primal;
    VecGetArray(u_primal, &local_u_primal);
1200

1201
1202
1203
1204
1205
    index = 0;
    for (int i = 0; i < nRankB * nComponents; i++)
      local_u_b[i] = local_b[index++];
    for (int i = 0; i < nRankPrimals * nComponents; i++)
      local_u_primal[i] = local_b[index++];
1206

1207
1208
1209
    VecRestoreArray(f, &local_b);
    VecRestoreArray(u_b, &local_u_b);
    VecRestoreArray(u_primal, &local_u_primal);
1210

1211
    recoverSolution(u_b, u_primal, vec);
1212

1213
1214
1215
1216
    VecDestroy(&u_b);
    VecDestroy(&u_primal);
    VecDestroy(&b);
    VecDestroy(&f);
1217

1218
    KSPDestroy(&ksp);
1219
  }
1220
1221


1222
1223
1224
  void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
1225

1226
    // === Create solver for the non primal (thus local) variables. ===
1227

1228
1229
1230
1231
1232
1233
1234
    KSPCreate(PETSC_COMM_WORLD, &ksp_b);
    KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_b, "solver_b_");
    KSPSetFromOptions(ksp_b);

    // RHS and solution vector.
    Vec vec_rhs;
1235

1236
1237
1238
1239
1240
1241
1242
1243
    // Some temporary vectors.
    Vec tmp_b0, tmp_b1, tmp_lagrange0, tmp_primal0, tmp_primal1;
    MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0);
    MatGetVecs(mat_lagrange, PETSC_NULL, &vec_rhs);
    MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b0);
    MatGetVecs(mat_b_b, PETSC_NULL, &tmp_b1);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal0);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal1);
1244
1245


1246
    // === Create new rhs ===
1247

1248
1249
1250
    // vec_rhs = L * inv(K_BB) * f_b
    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, vec_rhs);
1251

1252
1253
    // tmp_primal0 = M_PiB * inv(K_BB) * f_b
    MatMult(mat_primal_b, tmp_b0, tmp_primal0);
1254

1255
1256
    // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_b
    VecAXPBY(tmp_primal0, -1.0, 1.0, f_primal);
1257

1258
1259
    // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_b)
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1260

1261
1262
1263
1264
    //
    MatMult(mat_b_primal, tmp_primal0, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);
1265

1266
1267
    //
    VecAXPBY(vec_rhs, 1.0, 1.0, tmp_lagrange0);
1268
1269


1270
    // === Solve with FETI-DP operator. ===
1271

1272
    MSG("START FETI SOLVE!\n");
1273
    KSPSolve(ksp_feti, vec_rhs, vec_rhs);
1274
1275

   
1276
    // === Solve for u_primals. ===
1277

1278
    VecCopy(f_primal, tmp_primal0);
1279

1280
1281
    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1282

1283
    VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1);
1284

1285
1286
1287
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1288

1289
1290
    VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1291
1292

    
1293
    // === Solve for u_b. ===
1294

1295
1296
1297
    VecCopy(f_b, tmp_b0);
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1298

1299
1300
    MatMult(mat_b_primal, tmp_primal0, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1301

1302
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
1303
1304


1305
    // === And recover AMDiS solution vectors. ===
1306
    
1307
    recoverSolution(tmp_b0, tmp_primal0, vec);
1308
1309


1310
    // === Destroy all data structures. ===
1311
    
1312
1313
1314
1315
1316
1317
    VecDestroy(&vec_rhs);
    VecDestroy(&tmp_b0);
    VecDestroy(&tmp_b1);
    VecDestroy(&tmp_lagrange0);
    VecDestroy(&tmp_primal0);
    VecDestroy(&tmp_primal1);
1318
1319
	    

1320
    KSPDestroy(&ksp_b);
1321

1322
1323
1324
1325
1326
    MatDestroy(&mat_b_b);
    MatDestroy(&mat_primal_primal);
    MatDestroy(&mat_b_primal);
    MatDestroy(&mat_primal_b);
    MatDestroy(&mat_lagrange);
1327

1328
1329
    VecDestroy(&f_b);
    VecDestroy(&f_primal);
1330

1331
    destroySchurPrimalKsp();
1332

1333
1334
1335
1336
1337
    destroyFetiKsp();

    
    // === Destroy preconditioner data structures. ===

1338
1339
1340
1341
    MatDestroy(&mat_interior_interior);
    MatDestroy(&mat_bound_bound);
    MatDestroy(&mat_interior_bound);
    MatDestroy(&mat_bound_interior);
1342
  }
1343

1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357

  void PetscSolverFeti::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverFeti::solvePetscMatrix()");

    int debug = 0;
    Parameters::get("parallel->debug feti", debug);

    if (debug) {
      WARNING("FETI matrix is solved globally, thus without reducing to the lagrange multipliers!\n");

      solveFetiMatrix(vec);
    } else {
      solveReducedFetiMatrix(vec);
1358
    }      
1359
  }
1360

1361
1362
1363
#endif

}
For faster browsing, not all history is shown. View entire blame