PetscSolverFeti.cc 44.7 KB
Newer Older
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	if (!(*mat)[i][j])
	  continue;

	traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
	
	// Traverse all rows.
	for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
1012

1013
	  bool rowPrimal = primalDofMap[feSpace].isSet(*cursor);
1014
  
1015
1016
	  cols.clear();
	  colsOther.clear();
1017
	  values.clear();	  
1018
	  valuesOther.clear();
1019
1020
1021
1022
1023
1024

	  colsLocal.clear();
	  colsLocalOther.clear();
	  valuesLocal.clear();
	  valuesLocalOther.clear();

1025
1026
1027
1028
1029
	  
	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

1030
	    bool colPrimal = primalDofMap[feSpace].isSet(col(*icursor));
1031
1032

	    if (colPrimal) {
1033
	      if (rowPrimal) {
1034
		cols.push_back(col(*icursor));
1035
		values.push_back(value(*icursor));
1036
	      } else {
1037
		colsOther.push_back(col(*icursor));
1038
1039
1040
1041
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      if (rowPrimal) {
1042
		colsOther.push_back(col(*icursor));
1043
1044
		valuesOther.push_back(value(*icursor));
	      } else {
1045
		cols.push_back(col(*icursor));
1046
		values.push_back(value(*icursor));
1047
1048
	      }
	    }
1049
1050
1051
1052
1053


	    // === For preconditioner ===

	    if (!rowPrimal && !colPrimal) {
1054
1055
	      int rowIndex = localDofMap[feSpace][*cursor];
	      int colIndex = localDofMap[feSpace][col(*icursor)];
1056
1057
1058
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
1059
		  int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
1060
1061
1062
		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {
1063
		  int colIndex2 = (localDofMap[feSpace][col(*icursor)] - nLocalInterior) *
Thomas Witkowski's avatar
Thomas Witkowski committed
1064
		    nComponents + j;
1065
1066
1067
1068
1069
1070

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (colIndex < nLocalInterior) {
1071
		  int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
1072
1073
1074
		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
1075
		  int colIndex2 = (localDofMap[feSpace][col(*icursor)] - nLocalInterior) *
Thomas Witkowski's avatar
Thomas Witkowski committed
1076
		    nComponents + j;
1077
1078
1079
1080
1081
1082
1083
1084

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


1085
	  }  // for each nnz in row
1086

1087
1088
1089

	  // === Set matrix values. ===

1090
	  if (rowPrimal) {
1091
1092
	    int rowIndex = 
	      primalDofMap[feSpace][*cursor] * nComponents + i;
1093
1094
1095
1096

	    for (unsigned int k = 0; k < cols.size(); k++)
	      cols[k] = primalDofMap[feSpace][cols[k]] * nComponents + j;

1097
1098
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1099

1100
1101
1102
	    if (colsOther.size()) {
	      for (unsigned int k = 0; k < colsOther.size(); k++)
		colsOther[k] = localDofMap.mapGlobal(colsOther[k], j);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
1103
1104
 	      
	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
1105
1106
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
1107
	  } else {
1108
1109
	    int localRowIndex = localDofMap.mapLocal(*cursor, i);

Thomas Witkowski's avatar
Thomas Witkowski committed
1110
	    for (unsigned int k = 0; k < cols.size(); k++)
1111
	      cols[k] = localDofMap.mapLocal(cols[k], j);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
1112
	    
1113
1114
  	    MatSetValues(mat_b_b, 1, &localRowIndex, cols.size(),
  			 &(cols[0]), &(values[0]), ADD_VALUES);
1115

1116
	    if (colsOther.size()) {
1117
1118
1119
1120
1121
1122
1123
1124
	      int globalRowIndex = localDofMap.mapGlobal(*cursor, i);

	      for (unsigned int k = 0; k < colsOther.size(); k++)
		colsOther[k] = 
		  primalDofMap[feSpace][colsOther[k]] * nComponents + j;

  	      MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(),
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
1125
	    }
1126
	  }
1127

1128
	  // === Set matrix values for preconditioner ===
1129
1130

	  if (!rowPrimal) {
1131
1132
1133
	    switch (fetiPreconditioner) {
	    case FETI_DIRICHLET:
	      {
1134
		int rowIndex = localDofMap[feSpace][*cursor];
1135
1136
		
		if (rowIndex < nLocalInterior) {
1137
		  int rowIndex2 = localDofMap[feSpace][*cursor] * nComponents + i;
1138
1139
1140
1141
1142
1143
1144
1145
1146
		  
		  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 = 
1147
		    (localDofMap[feSpace][*cursor] - nLocalInterior) * nComponents + i;
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
		  
		  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;
1159

1160
1161
	    case FETI_LUMPED:
	      {
1162
		int rowIndex = localDofMap[feSpace][*cursor];
1163
1164
1165
		
		if (rowIndex >= nLocalInterior) {
		  int rowIndex2 = 
1166
		    (localDofMap[feSpace][*cursor] - nLocalInterior) * nComponents + i;
1167
1168
1169
1170
1171
1172
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);		
		}
	      }		
	      break;
1173
1174
	    default:
	      break;
1175
	    }	  
1176
	  }
1177
1178
1179
	} 
      }
    }
1180

1181
    // === Start global assembly procedure. ===
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194

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

1195

1196
    // === Start global assembly procedure for preconditioner matrices. ===
1197

1198
1199
1200
1201
    if (fetiPreconditioner != FETI_NONE) {
      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
    }
1202

1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
    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);
    }
1213
1214


1215
1216
    // === Create and fill PETSc matrix for Lagrange constraints. ===

1217
    createMatLagrange(feSpaces);
1218

1219
1220
1221
1222
1223
1224
1225
1226

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

    KSPCreate(PETSC_COMM_SELF, &ksp_b);
    KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(ksp_b, "solver_b_");
    KSPSetFromOptions(ksp_b);

1227
1228
1229
    
    // === Create PETSc solver for the Schur complement on primal variables. ===
    
1230
    createSchurPrimalKsp(feSpaces);
1231
1232
1233
1234


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

1235
    createFetiKsp(feSpaces);
1236
1237
  }

1238

1239
1240
1241
1242
  void PetscSolverFeti::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverFeti::fillPetscRhs()");

1243
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(vec);
1244
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
1245
    int nComponents = vec->getSize();
1246

1247
    VecCreateMPI(PETSC_COMM_WORLD, 
1248
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b);
1249
    VecCreateMPI(PETSC_COMM_WORLD, 
1250
1251
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
		 &f_primal);
1252
1253
1254
1255
1256
    
    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();
1257
	if (isPrimal(feSpace, index)) {
1258
	  index = primalDofMap[feSpace][index] * nComponents + i;
1259
	  double value = *dofIt;
1260
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
1261
	} else {
1262
	  index = localDofMap.mapGlobal(index, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
1263
	  VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
1264
1265
1266
1267
	}      
      }
    }

1268
1269
    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);
1270

1271
1272
    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);
1273
1274
  }

1275
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1276
1277
1278
  void PetscSolverFeti::solveLocalProblem(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverFeti::solveLocalProblem()");
1279

Thomas Witkowski's avatar
Thomas Witkowski committed
1280
    Vec tmp;
1281
    VecCreateSeq(PETSC_COMM_SELF, localDofMap.getRankDofs(), &tmp);
1282

1283
1284
    PetscScalar *tmpValues, *rhsValues;
    VecGetArray(tmp, &tmpValues);
Thomas Witkowski's avatar
Thomas Witkowski committed
1285
    VecGetArray(rhs, &rhsValues);
1286

1287
    for (int i = 0; i < localDofMap.getRankDofs(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
1288
      tmpValues[i] = rhsValues[i];
1289

Thomas Witkowski's avatar
Thomas Witkowski committed
1290
1291
    VecRestoreArray(rhs, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);
1292

Thomas Witkowski's avatar
Thomas Witkowski committed
1293
    KSPSolve(ksp_b, tmp, tmp);
1294

Thomas Witkowski's avatar
Thomas Witkowski committed
1295
    VecGetArray(tmp, &tmpValues);
1296
    VecGetArray(sol, &rhsValues);
1297

1298
1299
    for (int i = 0; i < localDofMap.getRankDofs(); i++) 
      rhsValues[i] = tmpValues[i];
1300

1301
    VecRestoreArray(sol, &rhsValues);
Thomas Witkowski's avatar
Thomas Witkowski committed
1302
    VecRestoreArray(tmp, &tmpValues);
1303

Thomas Witkowski's avatar
Thomas Witkowski committed
1304
    VecDestroy(&tmp);
1305
  }
1306
1307


Thomas Witkowski's avatar
Thomas Witkowski committed
1308
1309
1310
1311
  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");

1312
1313
1314
    // The function was removed when FETI-DP was reformulated for mixed finite
    // elements. To reconstruct this function, check for revision number 2024.
    ERROR_EXIT("This function was removed!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1315
1316
1317
  }


1318
1319
1320
  void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
1321

1322
    // RHS vector.
1323
    Vec vec_rhs;
1324

1325
1326
1327
1328
    // 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);
Thomas Witkowski's avatar
Thomas Witkowski committed
1329
1330
    VecDuplicate(f_b, &tmp_b0);
    VecDuplicate(f_b, &tmp_b1);
1331
1332
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal0);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal1);
1333
1334


1335
    // === Create new rhs ===
1336

1337
1338
    //    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
Thomas Witkowski's avatar
Thomas Witkowski committed
1339
    solveLocalProblem(f_b, tmp_b0);
1340
    MatMult(mat_lagrange, tmp_b0, vec_rhs);
1341

1342
    // tmp_primal0 = M_PiB * inv(K_BB) * f_B
1343
    MatMult(mat_primal_b, tmp_b0, tmp_primal0);
1344

1345
1346
    // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_B
    VecAXPBY(tmp_primal0, 1.0, -1.0, f_primal);
1347

1348
    // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B)
1349
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1350

1351
1352
    //
    MatMult(mat_b_primal, tmp_primal0, tmp_b0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1353
    solveLocalProblem(tmp_b0, tmp_b0);
1354
    MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);
1355

1356
    //
1357
    VecAXPBY(vec_rhs, -1.0, 1.0, tmp_lagrange0);
1358
1359


1360
1361
    // === Solve with FETI-DP operator. ===
    KSPSolve(ksp_feti, vec_rhs, vec_rhs);
1362

1363
    // === Solve for u_primals. ===
1364

1365
    VecCopy(f_primal, tmp_primal0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1366
    solveLocalProblem(f_b, tmp_b0);
1367
1368
1369
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
    VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1);
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1370
    solveLocalProblem(tmp_b0, tmp_b0);
1371
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1372

1373
1374
    VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1375
1376

    
1377
    // === Solve for u_b. ===
1378

1379
1380
1381
    VecCopy(f_b, tmp_b0);
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1382

1383
1384
    MatMult(mat_b_primal, tmp_primal0, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
Thomas Witkowski's avatar
Thomas Witkowski committed
1385
    solveLocalProblem(tmp_b0, tmp_b0);
1386

1387
    // === And recover AMDiS solution vectors. ===
1388
    
1389
    recoverSolution(tmp_b0, tmp_primal0, vec);
1390
1391


1392
    // === Destroy all data structures. ===
1393
    
1394
1395
1396
1397
1398
1399
    VecDestroy(&vec_rhs);
    VecDestroy(&tmp_b0);
    VecDestroy(&tmp_b1);
    VecDestroy(&tmp_lagrange0);
    VecDestroy(&tmp_primal0);
    VecDestroy(&tmp_primal1);
1400
	    
1401
1402
1403
    VecDestroy(&f_b);
    VecDestroy(&f_primal);
  }
1404

1405
1406
1407
1408

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

1410
1411
1412
1413
1414
    MatDestroy(&mat_b_b);
    MatDestroy(&mat_primal_primal);
    MatDestroy(&mat_b_primal);
    MatDestroy(&mat_primal_b);
    MatDestroy(&mat_lagrange);
1415

1416
1417
    // === Destroy preconditioner data structures. ===

1418
1419
1420
1421
1422
1423
1424
1425
    if (fetiPreconditioner != FETI_NONE)
      MatDestroy(&mat_duals_duals);

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

1427
1428
1429
1430
1431
1432
    destroySchurPrimalKsp();

    destroyFetiKsp();

    KSPDestroy(&ksp_b);
  }
1433

Thomas Witkowski's avatar
Thomas Witkowski committed
1434

1435
1436
1437
1438
1439
1440
1441
  void PetscSolverFeti::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverFeti::solvePetscMatrix()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1442
1443
1444
1445
1446
1447
    if (debug) {
      solveFetiMatrix(vec);
    } else {
      solveReducedFetiMatrix(vec);
    } 

1448
    MeshDistributor::globalMeshDistributor->synchVector(vec);
1449
  }
1450

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