1 #include <petsc/private/pcmgimpl.h> /*I "petscksp.h" I*/ 2 #include <petscdm.h> 3 4 static PetscErrorCode xfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 5 { 6 PetscInt k = *((PetscInt *) ctx), c; 7 8 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[0], k); 9 return 0; 10 } 11 static PetscErrorCode yfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 12 { 13 PetscInt k = *((PetscInt *) ctx), c; 14 15 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[1], k); 16 return 0; 17 } 18 static PetscErrorCode zfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 19 { 20 PetscInt k = *((PetscInt *) ctx), c; 21 22 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[2], k); 23 return 0; 24 } 25 static PetscErrorCode xsin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 26 { 27 PetscInt k = *((PetscInt *) ctx), c; 28 29 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI*(k+1)*coords[0]); 30 return 0; 31 } 32 static PetscErrorCode ysin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 33 { 34 PetscInt k = *((PetscInt *) ctx), c; 35 36 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI*(k+1)*coords[1]); 37 return 0; 38 } 39 static PetscErrorCode zsin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 40 { 41 PetscInt k = *((PetscInt *) ctx), c; 42 43 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI*(k+1)*coords[2]); 44 return 0; 45 } 46 47 PetscErrorCode DMSetBasisFunction_Internal(PetscInt Nf, PetscBool usePoly, PetscInt dir, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) 48 { 49 PetscInt f; 50 51 PetscFunctionBeginUser; 52 for (f = 0; f < Nf; ++f) { 53 if (usePoly) { 54 switch (dir) { 55 case 0: funcs[f] = xfunc;break; 56 case 1: funcs[f] = yfunc;break; 57 case 2: funcs[f] = zfunc;break; 58 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No function for direction %D", dir); 59 } 60 } else { 61 switch (dir) { 62 case 0: funcs[f] = xsin;break; 63 case 1: funcs[f] = ysin;break; 64 case 2: funcs[f] = zsin;break; 65 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No function for direction %D", dir); 66 } 67 } 68 } 69 PetscFunctionReturn(0); 70 } 71 72 static PetscErrorCode PCMGCreateCoarseSpaceDefault_Private(PC pc, PetscInt level, PCMGCoarseSpaceType cstype, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 73 { 74 PetscBool poly = cstype == PCMG_POLYNOMIAL ? PETSC_TRUE : PETSC_FALSE; 75 PetscErrorCode (**funcs)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar*,void*); 76 void **ctxs; 77 PetscInt dim, d, Nf, f, k; 78 79 PetscFunctionBegin; 80 PetscCall(DMGetCoordinateDim(dm, &dim)); 81 PetscCall(DMGetNumFields(dm, &Nf)); 82 PetscCheckFalse(Nc % dim,PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "The number of coarse vectors %D must be divisible by the dimension %D", Nc, dim); 83 PetscCall(PetscMalloc2(Nf, &funcs, Nf, &ctxs)); 84 if (!*coarseSpace) PetscCall(PetscCalloc1(Nc, coarseSpace)); 85 for (k = 0; k < Nc/dim; ++k) { 86 for (f = 0; f < Nf; ++f) {ctxs[f] = &k;} 87 for (d = 0; d < dim; ++d) { 88 if (!(*coarseSpace)[k*dim+d]) PetscCall(DMCreateGlobalVector(dm, &(*coarseSpace)[k*dim+d])); 89 PetscCall(DMSetBasisFunction_Internal(Nf, poly, d, funcs)); 90 PetscCall(DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, (*coarseSpace)[k*dim+d])); 91 } 92 } 93 PetscCall(PetscFree2(funcs, ctxs)); 94 PetscFunctionReturn(0); 95 } 96 97 static PetscErrorCode PCMGCreateCoarseSpace_Polynomial(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 98 { 99 PetscFunctionBegin; 100 PetscCall(PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_POLYNOMIAL, dm, ksp, Nc, initialGuess, coarseSpace)); 101 PetscFunctionReturn(0); 102 } 103 104 PetscErrorCode PCMGCreateCoarseSpace_Harmonic(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 105 { 106 PetscFunctionBegin; 107 PetscCall(PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_HARMONIC, dm, ksp, Nc, initialGuess, coarseSpace)); 108 PetscFunctionReturn(0); 109 } 110 111 /* 112 PCMGComputeCoarseSpace_Internal - Compute vectors on level l that must be accurately interpolated. 113 114 Input Parameters: 115 + pc - The PCMG 116 . l - The level 117 . Nc - The size of the space (number of vectors) 118 - cspace - The space from level l-1, or NULL 119 120 Output Parameter: 121 . space - The space which must be accurately interpolated. 122 123 Level: developer 124 125 Note: This space is normally used to adapt the interpolator. 126 127 .seealso: PCMGAdaptInterpolator_Private() 128 */ 129 PetscErrorCode PCMGComputeCoarseSpace_Internal(PC pc, PetscInt l, PCMGCoarseSpaceType cstype, PetscInt Nc, const Vec cspace[], Vec *space[]) 130 { 131 PetscErrorCode (*coarseConstructor)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec*[]); 132 DM dm; 133 KSP smooth; 134 135 PetscFunctionBegin; 136 switch (cstype) { 137 case PCMG_POLYNOMIAL: coarseConstructor = &PCMGCreateCoarseSpace_Polynomial;break; 138 case PCMG_HARMONIC: coarseConstructor = &PCMGCreateCoarseSpace_Harmonic;break; 139 case PCMG_EIGENVECTOR: 140 if (l > 0) PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_MEV", &coarseConstructor)); 141 else PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_EV", &coarseConstructor)); 142 break; 143 case PCMG_GENERALIZED_EIGENVECTOR: 144 if (l > 0) PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_MGEV", &coarseConstructor)); 145 else PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_GEV", &coarseConstructor)); 146 break; 147 default: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle coarse space type %D", cstype); 148 } 149 PetscCall(PCMGGetSmoother(pc, l, &smooth)); 150 PetscCall(KSPGetDM(smooth, &dm)); 151 PetscCall((*coarseConstructor)(pc, l, dm, smooth, Nc, cspace, space)); 152 PetscFunctionReturn(0); 153 } 154 155 /* 156 PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level 1 157 158 Input Parameters: 159 + pc - The PCMG 160 . l - The level l 161 . csmooth - The (coarse) smoother for level l-1 162 . fsmooth - The (fine) smoother for level l 163 . Nc - The size of the subspace used for adaptation 164 . cspace - The (coarse) vectors in the subspace for level l-1 165 - fspace - The (fine) vectors in the subspace for level l 166 167 Level: developer 168 169 Note: This routine resets the interpolation and restriction for level l. 170 171 .seealso: PCMGComputeCoarseSpace_Private() 172 */ 173 PetscErrorCode PCMGAdaptInterpolator_Internal(PC pc, PetscInt l, KSP csmooth, KSP fsmooth, PetscInt Nc, Vec cspace[], Vec fspace[]) 174 { 175 PC_MG *mg = (PC_MG *) pc->data; 176 DM dm, cdm; 177 Mat Interp, InterpAdapt; 178 179 PetscFunctionBegin; 180 /* There is no interpolator for the coarse level */ 181 if (!l) PetscFunctionReturn(0); 182 PetscCall(KSPGetDM(csmooth, &cdm)); 183 PetscCall(KSPGetDM(fsmooth, &dm)); 184 PetscCall(PCMGGetInterpolation(pc, l, &Interp)); 185 186 PetscCall(DMAdaptInterpolator(cdm, dm, Interp, fsmooth, Nc, fspace, cspace, &InterpAdapt, pc)); 187 if (mg->mespMonitor) PetscCall(DMCheckInterpolator(dm, InterpAdapt, Nc, cspace, fspace, 0.5/* PETSC_SMALL */)); 188 PetscCall(PCMGSetInterpolation(pc, l, InterpAdapt)); 189 PetscCall(PCMGSetRestriction(pc, l, InterpAdapt)); 190 PetscCall(MatDestroy(&InterpAdapt)); 191 PetscFunctionReturn(0); 192 } 193 194 /* 195 PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted 196 197 Input Parameters: 198 + pc - The PCMG 199 - l - The level l 200 201 Level: developer 202 203 Note: This routine recomputes the Galerkin triple product for the operator on level l. 204 */ 205 PetscErrorCode PCMGRecomputeLevelOperators_Internal(PC pc, PetscInt l) 206 { 207 Mat fA, fB; /* The system and preconditioning operators on level l+1 */ 208 Mat A, B; /* The system and preconditioning operators on level l */ 209 Mat Interp, Restrc; /* The interpolation operator from level l to l+1, and restriction operator from level l+1 to l */ 210 KSP smooth, fsmooth; /* The smoothers on levels l and l+1 */ 211 PCMGGalerkinType galerkin; /* The Galerkin projection flag */ 212 MatReuse reuse = MAT_REUSE_MATRIX; /* The matrices are always assumed to be present already */ 213 PetscBool doA = PETSC_FALSE; /* Updates the system operator */ 214 PetscBool doB = PETSC_FALSE; /* Updates the preconditioning operator (A == B, then update B) */ 215 PetscInt n; /* The number of multigrid levels */ 216 217 PetscFunctionBegin; 218 PetscCall(PCMGGetGalerkin(pc, &galerkin)); 219 if (galerkin >= PC_MG_GALERKIN_NONE) PetscFunctionReturn(0); 220 PetscCall(PCMGGetLevels(pc, &n)); 221 /* Do not recompute operator for the finest grid */ 222 if (l == n-1) PetscFunctionReturn(0); 223 PetscCall(PCMGGetSmoother(pc, l, &smooth)); 224 PetscCall(KSPGetOperators(smooth, &A, &B)); 225 PetscCall(PCMGGetSmoother(pc, l+1, &fsmooth)); 226 PetscCall(KSPGetOperators(fsmooth, &fA, &fB)); 227 PetscCall(PCMGGetInterpolation(pc, l+1, &Interp)); 228 PetscCall(PCMGGetRestriction(pc, l+1, &Restrc)); 229 if ((galerkin == PC_MG_GALERKIN_PMAT) || (galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE; 230 if ((galerkin == PC_MG_GALERKIN_MAT) || ((galerkin == PC_MG_GALERKIN_BOTH) && (fA != fB))) doA = PETSC_TRUE; 231 if (doA) PetscCall(MatGalerkin(Restrc, fA, Interp, reuse, 1.0, &A)); 232 if (doB) PetscCall(MatGalerkin(Restrc, fB, Interp, reuse, 1.0, &B)); 233 PetscFunctionReturn(0); 234 } 235