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: SETERRQ1(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: SETERRQ1(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 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 82 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 83 if (Nc % dim) SETERRQ2(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "The number of coarse vectors %D must be divisible by the dimension %D", Nc, dim); 84 ierr = PetscMalloc2(Nf, &funcs, Nf, &ctxs);CHKERRQ(ierr); 85 if (!*coarseSpace) {ierr = PetscCalloc1(Nc, coarseSpace);CHKERRQ(ierr);} 86 for (k = 0; k < Nc/dim; ++k) { 87 for (f = 0; f < Nf; ++f) {ctxs[f] = &k;} 88 for (d = 0; d < dim; ++d) { 89 if (!(*coarseSpace)[k*dim+d]) {ierr = DMCreateGlobalVector(dm, &(*coarseSpace)[k*dim+d]);CHKERRQ(ierr);} 90 ierr = DMSetBasisFunction_Internal(Nf, poly, d, funcs);CHKERRQ(ierr); 91 ierr = DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, (*coarseSpace)[k*dim+d]);CHKERRQ(ierr); 92 } 93 } 94 ierr = PetscFree2(funcs, ctxs);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 static PetscErrorCode PCMGCreateCoarseSpace_Polynomial(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 99 { 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 ierr = PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_POLYNOMIAL, dm, ksp, Nc, initialGuess, coarseSpace);CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 107 PetscErrorCode PCMGCreateCoarseSpace_Harmonic(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 108 { 109 PetscErrorCode ierr; 110 111 PetscFunctionBegin; 112 ierr = PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_HARMONIC, dm, ksp, Nc, initialGuess, coarseSpace);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 /* 117 PCMGComputeCoarseSpace_Internal - Compute vectors on level l that must be accurately interpolated. 118 119 Input Parameters: 120 + pc - The PCMG 121 . l - The level 122 . Nc - The size of the space (number of vectors) 123 - cspace - The space from level l-1, or NULL 124 125 Output Parameter: 126 . space - The space which must be accurately interpolated. 127 128 Level: developer 129 130 Note: This space is normally used to adapt the interpolator. 131 132 .seealso: PCMGAdaptInterpolator_Private() 133 */ 134 PetscErrorCode PCMGComputeCoarseSpace_Internal(PC pc, PetscInt l, PCMGCoarseSpaceType cstype, PetscInt Nc, const Vec cspace[], Vec *space[]) 135 { 136 PetscErrorCode (*coarseConstructor)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec*[]); 137 DM dm; 138 KSP smooth; 139 PetscErrorCode ierr; 140 141 PetscFunctionBegin; 142 switch (cstype) { 143 case PCMG_POLYNOMIAL: coarseConstructor = &PCMGCreateCoarseSpace_Polynomial;break; 144 case PCMG_HARMONIC: coarseConstructor = &PCMGCreateCoarseSpace_Harmonic;break; 145 case PCMG_EIGENVECTOR: 146 if (l > 0) {ierr = PCMGGetCoarseSpaceConstructor("BAMG_MEV", &coarseConstructor);CHKERRQ(ierr);} 147 else {ierr = PCMGGetCoarseSpaceConstructor("BAMG_EV", &coarseConstructor);CHKERRQ(ierr);} 148 break; 149 case PCMG_GENERALIZED_EIGENVECTOR: 150 if (l > 0) {ierr = PCMGGetCoarseSpaceConstructor("BAMG_MGEV", &coarseConstructor);CHKERRQ(ierr);} 151 else {ierr = PCMGGetCoarseSpaceConstructor("BAMG_GEV", &coarseConstructor);CHKERRQ(ierr);} 152 break; 153 default: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle coarse space type %D", cstype); 154 } 155 ierr = PCMGGetSmoother(pc, l, &smooth);CHKERRQ(ierr); 156 ierr = KSPGetDM(smooth, &dm);CHKERRQ(ierr); 157 ierr = (*coarseConstructor)(pc, l, dm, smooth, Nc, cspace, space);CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 /* 162 PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level 1 163 164 Input Parameters: 165 + pc - The PCMG 166 . l - The level l 167 . csmooth - The (coarse) smoother for level l-1 168 . fsmooth - The (fine) smoother for level l 169 . Nc - The size of the subspace used for adaptation 170 . cspace - The (coarse) vectors in the subspace for level l-1 171 - fspace - The (fine) vectors in the subspace for level l 172 173 Level: developer 174 175 Note: This routine resets the interpolation and restriction for level l. 176 177 .seealso: PCMGComputeCoarseSpace_Private() 178 */ 179 PetscErrorCode PCMGAdaptInterpolator_Internal(PC pc, PetscInt l, KSP csmooth, KSP fsmooth, PetscInt Nc, Vec cspace[], Vec fspace[]) 180 { 181 PC_MG *mg = (PC_MG *) pc->data; 182 DM dm, cdm; 183 Mat Interp, InterpAdapt; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 /* There is no interpolator for the coarse level */ 188 if (!l) PetscFunctionReturn(0); 189 ierr = KSPGetDM(csmooth, &cdm);CHKERRQ(ierr); 190 ierr = KSPGetDM(fsmooth, &dm);CHKERRQ(ierr); 191 ierr = PCMGGetInterpolation(pc, l, &Interp);CHKERRQ(ierr); 192 193 ierr = DMAdaptInterpolator(cdm, dm, Interp, fsmooth, Nc, fspace, cspace, &InterpAdapt, pc);CHKERRQ(ierr); 194 if (mg->mespMonitor) {ierr = DMCheckInterpolator(dm, InterpAdapt, Nc, cspace, fspace, 0.5/* PETSC_SMALL */);CHKERRQ(ierr);} 195 ierr = PCMGSetInterpolation(pc, l, InterpAdapt);CHKERRQ(ierr); 196 ierr = PCMGSetRestriction(pc, l, InterpAdapt);CHKERRQ(ierr); 197 ierr = MatDestroy(&InterpAdapt);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 /* 202 PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted 203 204 Input Parameters: 205 + pc - The PCMG 206 - l - The level l 207 208 Level: developer 209 210 Note: This routine recomputes the Galerkin triple product for the operator on level l. 211 */ 212 PetscErrorCode PCMGRecomputeLevelOperators_Internal(PC pc, PetscInt l) 213 { 214 Mat fA, fB; /* The system and preconditioning operators on level l+1 */ 215 Mat A, B; /* The system and preconditioning operators on level l */ 216 Mat Interp, Restrc; /* The interpolation operator from level l to l+1, and restriction operator from level l+1 to l */ 217 KSP smooth, fsmooth; /* The smoothers on levels l and l+1 */ 218 PCMGGalerkinType galerkin; /* The Galerkin projection flag */ 219 MatReuse reuse = MAT_REUSE_MATRIX; /* The matrices are always assumed to be present already */ 220 PetscBool doA = PETSC_FALSE; /* Updates the system operator */ 221 PetscBool doB = PETSC_FALSE; /* Updates the preconditioning operator (A == B, then update B) */ 222 PetscInt n; /* The number of multigrid levels */ 223 PetscErrorCode ierr; 224 225 PetscFunctionBegin; 226 ierr = PCMGGetGalerkin(pc, &galerkin);CHKERRQ(ierr); 227 if (galerkin >= PC_MG_GALERKIN_NONE) PetscFunctionReturn(0); 228 ierr = PCMGGetLevels(pc, &n);CHKERRQ(ierr); 229 /* Do not recompute operator for the finest grid */ 230 if (l == n-1) PetscFunctionReturn(0); 231 ierr = PCMGGetSmoother(pc, l, &smooth);CHKERRQ(ierr); 232 ierr = KSPGetOperators(smooth, &A, &B);CHKERRQ(ierr); 233 ierr = PCMGGetSmoother(pc, l+1, &fsmooth);CHKERRQ(ierr); 234 ierr = KSPGetOperators(fsmooth, &fA, &fB);CHKERRQ(ierr); 235 ierr = PCMGGetInterpolation(pc, l+1, &Interp);CHKERRQ(ierr); 236 ierr = PCMGGetRestriction(pc, l+1, &Restrc);CHKERRQ(ierr); 237 if ((galerkin == PC_MG_GALERKIN_PMAT) || (galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE; 238 if ((galerkin == PC_MG_GALERKIN_MAT) || ((galerkin == PC_MG_GALERKIN_BOTH) && (fA != fB))) doA = PETSC_TRUE; 239 if (doA) {ierr = MatGalerkin(Restrc, fA, Interp, reuse, 1.0, &A);CHKERRQ(ierr);} 240 if (doB) {ierr = MatGalerkin(Restrc, fB, Interp, reuse, 1.0, &B);CHKERRQ(ierr);} 241 PetscFunctionReturn(0); 242 } 243