PetscSolverFeti.cc 43.7 KB
Newer Older
1001
	    if (!rowPrimal && !colPrimal) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1002
1003
	      int rowIndex = localDofMap[feSpace][*cursor].local;
	      int colIndex = localDofMap[feSpace][col(*icursor)].local;
1004
1005
1006
		
	      if (rowIndex < nLocalInterior) {
		if (colIndex < nLocalInterior) {
1007
		  int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
1008
1009
1010
		  colsLocal.push_back(colIndex2);
		  valuesLocal.push_back(value(*icursor));
		} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1011
1012
		  int colIndex2 = 
		    (localDofMap[feSpace][col(*icursor)].local - nLocalInterior) *
Thomas Witkowski's avatar
Thomas Witkowski committed
1013
		    nComponents + j;
1014
1015
1016
1017
1018
1019

		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		}
	      } else {
		if (colIndex < nLocalInterior) {
1020
		  int colIndex2 = localDofMap.mapLocal(col(*icursor), j);
1021
1022
1023
		  colsLocalOther.push_back(colIndex2);
		  valuesLocalOther.push_back(value(*icursor));
		} else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1024
		  int colIndex2 = (localDofMap[feSpace][col(*icursor)].local - nLocalInterior) *
Thomas Witkowski's avatar
Thomas Witkowski committed
1025
		    nComponents + j;
1026
1027
1028
1029
1030
1031
1032
1033

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


1034
	  }  // for each nnz in row
1035

1036
1037
1038

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

1039
	  if (rowPrimal) {
1040
	    int rowIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
1041
	      primalDofMap[feSpace][*cursor].local * nComponents + i;
1042
1043

	    for (unsigned int k = 0; k < cols.size(); k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
1044
	      cols[k] = primalDofMap[feSpace][cols[k]].local * nComponents + j;
1045

1046
1047
	    MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);
1048

1049
1050
1051
	    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
1052
1053
 	      
	      MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
1054
1055
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
1056
	  } else {
1057
1058
	    int localRowIndex = localDofMap.mapLocal(*cursor, i);

Thomas Witkowski's avatar
Thomas Witkowski committed
1059
	    for (unsigned int k = 0; k < cols.size(); k++)
1060
	      cols[k] = localDofMap.mapLocal(cols[k], j);
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
1061
	    
1062
1063
  	    MatSetValues(mat_b_b, 1, &localRowIndex, cols.size(),
  			 &(cols[0]), &(values[0]), ADD_VALUES);
1064

1065
	    if (colsOther.size()) {
1066
1067
1068
1069
	      int globalRowIndex = localDofMap.mapGlobal(*cursor, i);

	      for (unsigned int k = 0; k < colsOther.size(); k++)
		colsOther[k] = 
Thomas Witkowski's avatar
Thomas Witkowski committed
1070
		  primalDofMap[feSpace][colsOther[k]].local * nComponents + j;
1071
1072
1073

  	      MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(),
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
1074
	    }
1075
	  }
1076

1077
	  // === Set matrix values for preconditioner ===
1078
1079

	  if (!rowPrimal) {
1080
1081
1082
	    switch (fetiPreconditioner) {
	    case FETI_DIRICHLET:
	      {
Thomas Witkowski's avatar
Thomas Witkowski committed
1083
		int rowIndex = localDofMap[feSpace][*cursor].local;
1084
1085
		
		if (rowIndex < nLocalInterior) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1086
		  int rowIndex2 = localDofMap[feSpace][*cursor].local * nComponents + i;
1087
1088
1089
1090
1091
1092
1093
1094
1095
		  
		  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 = 
Thomas Witkowski's avatar
Thomas Witkowski committed
1096
		    (localDofMap[feSpace][*cursor].local - nLocalInterior) * nComponents + i;
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
		  
		  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;
1108

1109
1110
	    case FETI_LUMPED:
	      {
Thomas Witkowski's avatar
Thomas Witkowski committed
1111
		int rowIndex = localDofMap[feSpace][*cursor].local;
1112
1113
1114
		
		if (rowIndex >= nLocalInterior) {
		  int rowIndex2 = 
Thomas Witkowski's avatar
Thomas Witkowski committed
1115
		    (localDofMap[feSpace][*cursor].local - nLocalInterior) * nComponents + i;
1116
1117
1118
1119
1120
1121
		  
		  MatSetValues(mat_duals_duals, 1, &rowIndex2, colsLocal.size(),
			       &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);		
		}
	      }		
	      break;
1122
1123
	    default:
	      break;
1124
	    }	  
1125
	  }
1126
1127
1128
	} 
      }
    }
1129

1130
    // === Start global assembly procedure. ===
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143

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

1144

1145
    // === Start global assembly procedure for preconditioner matrices. ===
1146

1147
1148
1149
1150
    if (fetiPreconditioner != FETI_NONE) {
      MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY); 
    }
1151

1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
    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);
    }
1162
1163


1164
1165
    // === Create and fill PETSc matrix for Lagrange constraints. ===

1166
    createMatLagrange(feSpaces);
1167

1168
1169
1170
1171
1172
1173
1174
1175

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

1176
1177
1178
    
    // === Create PETSc solver for the Schur complement on primal variables. ===
    
1179
    createSchurPrimalKsp(feSpaces);
1180
1181
1182
1183


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

1184
    createFetiKsp(feSpaces);
1185
1186
  }

1187

1188
1189
1190
1191
  void PetscSolverFeti::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverFeti::fillPetscRhs()");

1192
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(vec);
1193
    const FiniteElemSpace *feSpace = meshDistributor->getFeSpace(0);
1194
    int nComponents = vec->getSize();
1195

1196
    VecCreateMPI(PETSC_COMM_WORLD, 
1197
		 localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b);
1198
    VecCreateMPI(PETSC_COMM_WORLD, 
1199
1200
		 primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
		 &f_primal);
1201
1202
1203
1204
1205
    
    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();
1206
	if (isPrimal(feSpace, index)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1207
	  index = primalDofMap[feSpace][index].local * nComponents + i;
1208
	  double value = *dofIt;
1209
	  VecSetValues(f_primal, 1, &index, &value, ADD_VALUES);
1210
	} else {
1211
	  index = localDofMap.mapGlobal(index, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
1212
	  VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
1213
1214
1215
1216
	}      
      }
    }

1217
1218
    VecAssemblyBegin(f_b);
    VecAssemblyEnd(f_b);
1219

1220
1221
    VecAssemblyBegin(f_primal);
    VecAssemblyEnd(f_primal);
1222
1223
  }

1224
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1225
1226
1227
  void PetscSolverFeti::solveLocalProblem(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverFeti::solveLocalProblem()");
1228

Thomas Witkowski's avatar
Thomas Witkowski committed
1229
    Vec tmp;
1230
    VecCreateSeq(PETSC_COMM_SELF, localDofMap.getRankDofs(), &tmp);
1231

1232
1233
    PetscScalar *tmpValues, *rhsValues;
    VecGetArray(tmp, &tmpValues);
Thomas Witkowski's avatar
Thomas Witkowski committed
1234
    VecGetArray(rhs, &rhsValues);
1235

1236
    for (int i = 0; i < localDofMap.getRankDofs(); i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
1237
      tmpValues[i] = rhsValues[i];
1238

Thomas Witkowski's avatar
Thomas Witkowski committed
1239
1240
    VecRestoreArray(rhs, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);
1241

Thomas Witkowski's avatar
Thomas Witkowski committed
1242
    KSPSolve(ksp_b, tmp, tmp);
1243

Thomas Witkowski's avatar
Thomas Witkowski committed
1244
    VecGetArray(tmp, &tmpValues);
1245
    VecGetArray(sol, &rhsValues);
1246

1247
1248
    for (int i = 0; i < localDofMap.getRankDofs(); i++) 
      rhsValues[i] = tmpValues[i];
1249

1250
    VecRestoreArray(sol, &rhsValues);
Thomas Witkowski's avatar
Thomas Witkowski committed
1251
    VecRestoreArray(tmp, &tmpValues);
1252

Thomas Witkowski's avatar
Thomas Witkowski committed
1253
    VecDestroy(&tmp);
1254
  }
1255
1256


Thomas Witkowski's avatar
Thomas Witkowski committed
1257
1258
1259
1260
  void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveFetiMatrix()");

1261
1262
1263
    // 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
1264
1265
1266
  }


1267
1268
1269
  void PetscSolverFeti::solveReducedFetiMatrix(SystemVector &vec)
  {
    FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
1270

1271
    // RHS vector.
1272
    Vec vec_rhs;
1273

1274
1275
1276
1277
    // 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
1278
1279
    VecDuplicate(f_b, &tmp_b0);
    VecDuplicate(f_b, &tmp_b1);
1280
1281
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal0);
    MatGetVecs(mat_primal_primal, PETSC_NULL, &tmp_primal1);
1282
1283


1284
    // === Create new rhs ===
1285

1286
1287
    //    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
1288
    solveLocalProblem(f_b, tmp_b0);
1289
    MatMult(mat_lagrange, tmp_b0, vec_rhs);
1290

1291
    // tmp_primal0 = M_PiB * inv(K_BB) * f_B
1292
    MatMult(mat_primal_b, tmp_b0, tmp_primal0);
1293

1294
1295
    // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_B
    VecAXPBY(tmp_primal0, 1.0, -1.0, f_primal);
1296

1297
    // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B)
1298
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1299

1300
1301
    //
    MatMult(mat_b_primal, tmp_primal0, tmp_b0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1302
    solveLocalProblem(tmp_b0, tmp_b0);
1303
    MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);
1304

1305
    //
1306
    VecAXPBY(vec_rhs, -1.0, 1.0, tmp_lagrange0);
1307
1308


1309
1310
    // === Solve with FETI-DP operator. ===
    KSPSolve(ksp_feti, vec_rhs, vec_rhs);
1311

1312
    // === Solve for u_primals. ===
1313

1314
    VecCopy(f_primal, tmp_primal0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1315
    solveLocalProblem(f_b, tmp_b0);
1316
1317
1318
    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
1319
    solveLocalProblem(tmp_b0, tmp_b0);
1320
    MatMult(mat_primal_b, tmp_b0, tmp_primal1);
1321

1322
1323
    VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
    KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
1324
1325

    
1326
    // === Solve for u_b. ===
1327

1328
1329
1330
    VecCopy(f_b, tmp_b0);
    MatMultTranspose(mat_lagrange, vec_rhs, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
1331

1332
1333
    MatMult(mat_b_primal, tmp_primal0, tmp_b1);
    VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
Thomas Witkowski's avatar
Thomas Witkowski committed
1334
    solveLocalProblem(tmp_b0, tmp_b0);
1335

1336
    // === And recover AMDiS solution vectors. ===
1337
    
1338
    recoverSolution(tmp_b0, tmp_primal0, vec);
1339
1340


1341
    // === Destroy all data structures. ===
1342
    
1343
1344
1345
1346
1347
1348
    VecDestroy(&vec_rhs);
    VecDestroy(&tmp_b0);
    VecDestroy(&tmp_b1);
    VecDestroy(&tmp_lagrange0);
    VecDestroy(&tmp_primal0);
    VecDestroy(&tmp_primal1);
1349
	    
1350
1351
1352
    VecDestroy(&f_b);
    VecDestroy(&f_primal);
  }
1353

1354
1355
1356
1357

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

1359
1360
1361
1362
1363
    MatDestroy(&mat_b_b);
    MatDestroy(&mat_primal_primal);
    MatDestroy(&mat_b_primal);
    MatDestroy(&mat_primal_b);
    MatDestroy(&mat_lagrange);
1364

1365
1366
    // === Destroy preconditioner data structures. ===

1367
1368
1369
1370
1371
1372
1373
1374
    if (fetiPreconditioner != FETI_NONE)
      MatDestroy(&mat_duals_duals);

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

1376
1377
1378
1379
1380
1381
    destroySchurPrimalKsp();

    destroyFetiKsp();

    KSPDestroy(&ksp_b);
  }
1382

Thomas Witkowski's avatar
Thomas Witkowski committed
1383

1384
1385
1386
1387
1388
1389
1390
  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
1391
1392
1393
1394
1395
1396
    if (debug) {
      solveFetiMatrix(vec);
    } else {
      solveReducedFetiMatrix(vec);
    } 

1397
    MeshDistributor::globalMeshDistributor->synchVector(vec);
1398
  }
1399

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