PetscSolverFeti.cc 45.7 KB
Newer Older
1001
1002
1003
1004
1005
1006
1007
1008
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      // Column is not a primal variable.

	      TEST_EXIT_DBG(globalIndexB.count(col(*icursor)))
		("No global B index for DOF %d!\n", col(*icursor));
1009
	      
1010
1011
1012
1013
1014
1015
1016
1017
	      int colIndex = globalIndexB[col(*icursor)] * nComponents + j;

	      if (rowPrimal) {
		colsOther.push_back(colIndex);
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(colIndex);
		values.push_back(value(*icursor));
1018
1019
	      }
	    }
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060



	    // === For preconditioner ===

	    if (!rowPrimal && !colPrimal) {
	      int rowIndex = globalIndexB[*cursor] - rStartB;
	      int colIndex = globalIndexB[col(*icursor)] - rStartB;
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;

		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (colIndex < nLocalInterior) {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
		  int colIndex2 = 
		    (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;

		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		}
	      }		
	    }


1061
	  }
1062

1063
1064
1065
	  if (rowPrimal) {
	    TEST_EXIT_DBG(globalPrimalIndex.count(*cursor))
	      ("Should not happen!\n");
1066

1067
1068
1069
	    int rowIndex = globalPrimalIndex[*cursor] * nComponents + i;
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1070

1071
1072
1073
1074
1075
1076
	    if (colsOther.size())
	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  } else {
	    TEST_EXIT_DBG(globalIndexB.count(*cursor))
	      ("Should not happen!\n");
1077

1078
1079
1080
	    int rowIndex = globalIndexB[*cursor] * nComponents + i;
	    MatSetValues(mat_b_b, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1081

1082
1083
1084
1085
	    if (colsOther.size())
	      MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(),
			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	  }
1086
1087
1088
1089

	  // === For preconditioner ===

	  if (!rowPrimal) {
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
	    switch (fetiPreconditioner) {
	    case FETI_DIRICHLET:
	      {
		int rowIndex = globalIndexB[*cursor] - rStartB;
		
		if (rowIndex < nLocalInterior) {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB) * nComponents + i;
		  
		  MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
		  
		  if (colsLocalOther.size()) 
		    MatSetValues(mat_interior_duals, 1, &rowIndex2, colsLocalOther.size(),
				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
		} else {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
		  
		  if (colsLocalOther.size()) 
		    MatSetValues(mat_duals_interior, 1, &rowIndex2, colsLocalOther.size(),
				 &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
		  
		}
	      }
	      break;
1119

1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
	    case FETI_LUMPED:
	      {
		int rowIndex = globalIndexB[*cursor] - rStartB;
		
		if (rowIndex >= nLocalInterior) {
		  int rowIndex2 = 
		    (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);		
		}
	      }		
	      break;
	    }	  
1134
	  }
1135
1136
1137
	} 
      }
    }
1138

1139
    // === Start global assembly procedure. ===
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152

    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);

1153

1154
    // === Start global assembly procedure for preconditioner matrices. ===
1155

1156
1157
1158
1159
    if (fetiPreconditioner != FETI_NONE) {
      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
    }
1160

1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY);
    }
1171
1172


1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
    // === Create and fill PETSc matrix for Lagrange constraints. ===

    createMatLagrange();

    
    // === Create PETSc solver for the Schur complement on primal variables. ===
    
    createSchurPrimalKsp();


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

    createFetiKsp();

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

    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);
  }

1195

1196
1197
1198
1199
1200
  void PetscSolverFeti::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverFeti::fillPetscRhs()");

    int nComponents = vec->getSize();
1201

1202
1203
1204
1205
1206
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankB * nComponents, nOverallB * nComponents, &f_b);
    VecCreateMPI(PETSC_COMM_WORLD, 
		 nRankPrimals * nComponents, 
		 nOverallPrimals * nComponents, &f_primal);
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
    
    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;
1218
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
1219
1220
1221
1222
1223
1224
	} else {
	  TEST_EXIT_DBG(globalIndexB.count(index))
	    ("Should not happen!\n");

	  index = globalIndexB[index] * nComponents + i;
	  double value = *dofIt;
1225
	  VecSetValues(f_b, 1, &index, &value, ADD_VALUES);
1226
1227
1228
1229
	}      
      }
    }

1230
1231
    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);
1232

1233
1234
    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);
1235
1236
1237
  }


1238
  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
1239
  {
1240
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
1241

1242
    Mat  mat_lagrange_transpose;
1243
    MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
1244

1245
    // === Create nested matrix which will contain the overall FETI system. ===
1246

1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
    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;
1258

1259
    MatCreateNest(PETSC_COMM_WORLD, 3, PETSC_NULL, 3, PETSC_NULL, &(nestedA[0][0]), &A);
1260

1261
1262
1263
    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  
1264
1265


1266
1267
1268
    int nRankNest = (nRankB + nRankPrimals + nRankLagrange) * nComponents;
    int nOverallNest = (nOverallB + nOverallPrimals + nOverallLagrange) * nComponents;
    int rStartNest = (rStartB + rStartPrimals + rStartLagrange) * nComponents;
1269

1270
1271
    {
      // === Test some matrix sizes. ===
1272

1273
1274
1275
1276
1277
      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");
1278

1279
1280
1281
      MatGetOwnershipRange(A, &matRow, &matCol);
      TEST_EXIT_DBG(matRow == rStartNest)("Should not happen!\n");
    }
1282

1283
1284
    // === Create rhs and solution vectors for the overall FETI system. ===

1285
1286
    Vec f, b;
    VecCreateMPI(PETSC_COMM_WORLD, nRankNest, nOverallNest, &f);
1287
1288
1289
1290
    VecDuplicate(f, &b);

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

1292
1293
    PetscScalar *local_f_b;
    VecGetArray(f_b, &local_f_b);
1294

1295
1296
    PetscScalar *local_f_primal;
    VecGetArray(f_primal, &local_f_primal);
1297

1298
1299
1300
1301
1302
1303
1304
    {
      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");
    }
1305

1306
1307
    PetscScalar *local_f;
    VecGetArray(f, &local_f);
1308

1309
1310
1311
1312
1313
    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];
1314

1315
1316
1317
    VecRestoreArray(f, &local_f);  
    VecRestoreArray(f_b, &local_f_b);
    VecRestoreArray(f_primal, &local_f_primal);
1318

1319
1320
    
    // === Create solver and solve the overall FETI system. ===
1321

1322
1323
1324
1325
    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD, &ksp);
    KSPSetOperators(ksp, A, A, SAME_NONZERO_PATTERN);
    KSPSetFromOptions(ksp);
1326
1327


1328
    KSPSolve(ksp, f, b);
1329
1330


1331
1332
1333
1334
1335
    // === Reconstruct FETI solution vectors. ===
    
    Vec u_b, u_primal;
    VecDuplicate(f_b, &u_b);
    VecDuplicate(f_primal, &u_primal);
1336
1337
    

1338
1339
    PetscScalar *local_b;
    VecGetArray(b, &local_b);
1340

1341
1342
    PetscScalar *local_u_b;
    VecGetArray(u_b, &local_u_b);
1343

1344
1345
    PetscScalar *local_u_primal;
    VecGetArray(u_primal, &local_u_primal);
1346

1347
1348
1349
1350
1351
    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++];
1352

1353
1354
1355
    VecRestoreArray(f, &local_b);
    VecRestoreArray(u_b, &local_u_b);
    VecRestoreArray(u_primal, &local_u_primal);
1356

1357
    recoverSolution(u_b, u_primal, vec);
1358

1359
1360
1361
1362
    VecDestroy(&u_b);
    VecDestroy(&u_primal);
    VecDestroy(&b);
    VecDestroy(&f);
1363

1364
    KSPDestroy(&ksp);
1365
  }
1366
1367


1368
1369
1370
  void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
1371

1372
    // RHS vector.
1373
    Vec vec_rhs;
1374

1375
1376
1377
1378
1379
1380
1381
1382
    // 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);
1383
1384


1385
    // === Create new rhs ===
1386

1387
1388
1389
    //    d = L inv(K_BB) f_B - L inv(K_BB) K_BPi inv(S_PiPi) [f_Pi - K_PiB inv(K_BB) f_B]

    // vec_rhs = L * inv(K_BB) * f_B
1390
1391
    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, vec_rhs);
1392

1393
    // tmp_primal0 = M_PiB * inv(K_BB) * f_B
1394
    MatMult(mat_primal_b, tmp_b0, tmp_primal0);
1395

1396
1397
    // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_B
    VecAXPBY(tmp_primal0, 1.0, -1.0, f_primal);
1398

1399
    // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B)
1400
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1401

1402
1403
1404
1405
    //
    MatMult(mat_b_primal, tmp_primal0, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);
1406

1407
    //
1408
    VecAXPBY(vec_rhs, -1.0, 1.0, tmp_lagrange0);
1409
1410


1411
    // === Solve with FETI-DP operator. ===
1412

1413
    KSPSolve(ksp_feti, vec_rhs, vec_rhs);
1414
   
1415
    // === Solve for u_primals. ===
1416

1417
    VecCopy(f_primal, tmp_primal0);
1418

1419
1420
    KSPSolve(ksp_b, f_b, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1421

1422
    VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1);
1423

1424
1425
1426
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0);
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1427

1428
1429
    VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1430
1431

    
1432
    // === Solve for u_b. ===
1433

1434
1435
1436
    VecCopy(f_b, tmp_b0);
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1437

1438
1439
    MatMult(mat_b_primal, tmp_primal0, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1440

1441
    KSPSolve(ksp_b, tmp_b0, tmp_b0);
1442

1443
    // === And recover AMDiS solution vectors. ===
1444
    
1445
    recoverSolution(tmp_b0, tmp_primal0, vec);
1446
1447


1448
    // === Destroy all data structures. ===
1449
    
1450
1451
1452
1453
1454
1455
    VecDestroy(&vec_rhs);
    VecDestroy(&tmp_b0);
    VecDestroy(&tmp_b1);
    VecDestroy(&tmp_lagrange0);
    VecDestroy(&tmp_primal0);
    VecDestroy(&tmp_primal1);
1456
	    
1457
1458
1459
    VecDestroy(&f_b);
    VecDestroy(&f_primal);
  }
1460

1461
1462
1463
1464

  void PetscSolverFeti::destroyMatrixData()
  {
    FUNCNAME("PetscSolverFeti::destroyMatrixData()");
1465

1466
1467
1468
1469
1470
    MatDestroy(&mat_b_b);
    MatDestroy(&mat_primal_primal);
    MatDestroy(&mat_b_primal);
    MatDestroy(&mat_primal_b);
    MatDestroy(&mat_lagrange);
1471

1472
1473
    // === Destroy preconditioner data structures. ===

1474
1475
1476
1477
1478
1479
1480
1481
    if (fetiPreconditioner != FETI_NONE)
      MatDestroy(&mat_duals_duals);

    if (fetiPreconditioner == FETI_DIRICHLET) {
      MatDestroy(&mat_interior_interior);
      MatDestroy(&mat_interior_duals);
      MatDestroy(&mat_duals_interior);
    }
1482

1483
1484
1485
1486
1487
1488
    destroySchurPrimalKsp();

    destroyFetiKsp();

    KSPDestroy(&ksp_b);
  }
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502

  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);
1503
1504
1505
    } 

    MeshDistributor::globalMeshDistributor->synchVector(vec);
1506
  }
1507

1508

1509
1510
1511
#endif

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