1 #include <petsc-private/snesimpl.h> /*I "petscsnes.h" I*/ 2 #include <petscdmshell.h> 3 4 typedef struct { 5 PetscInt n; /* local subdomains */ 6 SNES *subsnes; /* nonlinear solvers for each subdomain */ 7 8 Vec *r; /* function vectors */ 9 Vec *x; /* solution vectors */ 10 Vec *b; /* rhs vectors */ 11 12 IS *ois; 13 IS *iis; 14 PetscBool usesdm; /* use the outer DM's communication facilities rather than ISes */ 15 } SNES_NASM; 16 17 #undef __FUNCT__ 18 #define __FUNCT__ "SNESReset_NASM" 19 PetscErrorCode SNESReset_NASM(SNES snes) 20 { 21 SNES_NASM *nasm = (SNES_NASM *)snes->data; 22 PetscErrorCode ierr; 23 PetscInt i; 24 PetscFunctionBegin; 25 for (i=0;i<nasm->n;i++) { 26 if (nasm->x) { ierr = VecDestroy(&nasm->x[i]);CHKERRQ(ierr); } 27 if (nasm->r) { ierr = VecDestroy(&nasm->r[i]);CHKERRQ(ierr); } 28 if (nasm->b) { ierr = VecDestroy(&nasm->b[i]);CHKERRQ(ierr); } 29 30 if (nasm->subsnes) { ierr = SNESDestroy(&nasm->subsnes[i]);CHKERRQ(ierr); } 31 if (nasm->ois) { ierr = ISDestroy(&nasm->ois[i]);CHKERRQ(ierr); } 32 if (nasm->iis) { ierr = ISDestroy(&nasm->iis[i]);CHKERRQ(ierr); } 33 } 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "SNESDestroy_NASM" 39 PetscErrorCode SNESDestroy_NASM(SNES snes) 40 { 41 PetscErrorCode ierr; 42 PetscFunctionBegin; 43 ierr = SNESReset_NASM(snes);CHKERRQ(ierr); 44 ierr = PetscFree(snes->data); 45 PetscFunctionReturn(0); 46 } 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "SNESSetUp_NASM" 50 PetscErrorCode SNESSetUp_NASM(SNES snes) 51 { 52 SNES_NASM *nasm = (SNES_NASM *)snes->data; 53 PetscErrorCode ierr; 54 DM dm; 55 DM subdm; 56 PetscErrorCode (*f)(SNES,Vec,Vec,void*); 57 58 Mat A; 59 PetscErrorCode (*fj)(SNES,Vec,Mat*,Mat*,MatStructure*,void*); 60 61 PetscInt n; 62 void *ctx; 63 const char *optionsprefix; 64 65 PetscFunctionBegin; 66 67 if (!nasm->subsnes) { 68 if (snes->dm) { 69 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 70 nasm->n = 1; 71 nasm->usesdm = PETSC_TRUE; 72 /* create a single solver */ 73 ierr = PetscMalloc(sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 74 ierr = PetscMalloc(sizeof(Vec),&nasm->r);CHKERRQ(ierr); 75 ierr = PetscMalloc(sizeof(Vec),&nasm->x);CHKERRQ(ierr); 76 ierr = PetscMalloc(sizeof(Vec),&nasm->b);CHKERRQ(ierr); 77 ierr = SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[0]);CHKERRQ(ierr); 78 79 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 80 ierr = SNESAppendOptionsPrefix(nasm->subsnes[0],optionsprefix);CHKERRQ(ierr); 81 ierr = SNESAppendOptionsPrefix(nasm->subsnes[0],"sub_");CHKERRQ(ierr); 82 83 ierr = SNESGetDM(nasm->subsnes[0],&subdm);CHKERRQ(ierr); 84 85 /* set up the local function */ 86 ierr = DMSNESGetBlockFunction(dm,&f,&ctx);CHKERRQ(ierr); 87 ierr = DMCreateLocalVector(dm,&nasm->r[0]);CHKERRQ(ierr); 88 ierr = DMShellSetGlobalVector(dm,nasm->r[0]);CHKERRQ(ierr); 89 ierr = DMCreateLocalVector(dm,&nasm->b[0]);CHKERRQ(ierr); 90 ierr = DMCreateLocalVector(dm,&nasm->x[0]);CHKERRQ(ierr); 91 92 if (!f) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to provide a block solve!");CHKERRQ(ierr); 93 94 ierr = SNESSetFunction(nasm->subsnes[0],nasm->r[0],f,ctx);CHKERRQ(ierr); 95 96 /* set up the local jacobian -- TODO: do this correctly */ 97 ierr = DMSNESGetBlockJacobian(dm,&fj,&ctx);CHKERRQ(ierr); 98 ierr = VecGetSize(nasm->r[0],&n);CHKERRQ(ierr); 99 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,PETSC_DEFAULT,PETSC_NULL,&A);CHKERRQ(ierr); 100 ierr = DMShellSetMatrix(subdm,A);CHKERRQ(ierr); 101 ierr = SNESSetJacobian(nasm->subsnes[0],A,A,fj,ctx);CHKERRQ(ierr); 102 103 ierr = SNESSetFromOptions(nasm->subsnes[0]);CHKERRQ(ierr); 104 } else { 105 SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!"); 106 } 107 } else { 108 SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!"); 109 } 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "SNESSetFromOptions_NASM" 115 PetscErrorCode SNESSetFromOptions_NASM(SNES snes) 116 { 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = PetscOptionsHead("SNES NASM options");CHKERRQ(ierr); 121 ierr = PetscOptionsTail();CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "SNESView_NASM" 127 PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer) 128 { 129 PetscFunctionBegin; 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "SNESNASMSetSubdomains" 135 PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],IS iis[],IS ois[]) { 136 PetscErrorCode ierr; 137 PetscErrorCode (*f)(SNES,PetscInt,SNES*,IS*,IS*); 138 PetscFunctionBegin; 139 ierr = PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",(void (**)(void))&f);CHKERRQ(ierr); 140 ierr = (f)(snes,n,subsnes,iis,ois);CHKERRQ(ierr); 141 PetscFunctionReturn(0); 142 } 143 144 EXTERN_C_BEGIN 145 #undef __FUNCT__ 146 #define __FUNCT__ "SNESNASMSetSubdomains_NASM" 147 PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],IS iis[],IS ois[]) { 148 PetscInt i; 149 PetscErrorCode ierr; 150 SNES_NASM *nasm = (SNES_NASM *)snes->data; 151 PetscFunctionBegin; 152 if (snes->setupcalled)SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp()."); 153 154 if (!snes->setupcalled) { 155 nasm->n = n; 156 nasm->ois = 0; 157 nasm->iis = 0; 158 if (ois) { 159 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)ois[i]);CHKERRQ(ierr);} 160 } 161 if (iis) { 162 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)iis[i]);CHKERRQ(ierr);} 163 } 164 if (ois) { 165 ierr = PetscMalloc(n*sizeof(IS),&nasm->ois);CHKERRQ(ierr); 166 for (i=0; i<n; i++) { 167 nasm->ois[i] = ois[i]; 168 } 169 if (!iis) { 170 ierr = PetscMalloc(n*sizeof(IS),&nasm->iis);CHKERRQ(ierr); 171 for (i=0; i<n; i++) { 172 for (i=0; i<n; i++) {ierr = PetscObjectReference((PetscObject)ois[i]);CHKERRQ(ierr);} 173 nasm->iis[i] = ois[i]; 174 } 175 } 176 } 177 if (iis) { 178 ierr = PetscMalloc(n*sizeof(IS),&nasm->iis);CHKERRQ(ierr); 179 for (i=0; i<n; i++) { 180 nasm->iis[i] = iis[i]; 181 } 182 if (!ois) { 183 ierr = PetscMalloc(n*sizeof(IS),&nasm->ois);CHKERRQ(ierr); 184 for (i=0; i<n; i++) { 185 for (i=0; i<n; i++) { 186 ierr = PetscObjectReference((PetscObject)iis[i]);CHKERRQ(ierr); 187 nasm->ois[i] = iis[i]; 188 } 189 } 190 } 191 } 192 } 193 if (subsnes) { 194 ierr = PetscMalloc(n*sizeof(SNES),&nasm->subsnes);CHKERRQ(ierr); 195 for (i=0; i<n; i++) { 196 nasm->subsnes[i] = subsnes[i]; 197 } 198 } 199 PetscFunctionReturn(0); 200 } 201 EXTERN_C_END 202 203 #undef __FUNCT__ 204 #define __FUNCT__ "SNESNASMSolveLocal_Private" 205 PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec X) { 206 SNES_NASM *nasm = (SNES_NASM *)snes->data; 207 PetscInt i; 208 PetscErrorCode ierr; 209 Vec Xl,Bl; 210 DM dm; 211 PetscFunctionBegin; 212 /* restrict to the local system */ 213 if (nasm->usesdm) { 214 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 215 ierr = DMGlobalToLocalBegin(dm,X,INSERT_VALUES,nasm->x[0]);CHKERRQ(ierr); 216 ierr = DMGlobalToLocalEnd(dm,X,INSERT_VALUES,nasm->x[0]);CHKERRQ(ierr); 217 if (B) { 218 ierr = DMGlobalToLocalBegin(dm,B,INSERT_VALUES,nasm->b[0]);CHKERRQ(ierr); 219 ierr = DMGlobalToLocalEnd(dm,B,INSERT_VALUES,nasm->b[0]);CHKERRQ(ierr); 220 } 221 } else { 222 /* 223 for (i = 0;i < nasm->n;i++) { 224 ierr = VecScatterBegin(nasm->gorestriction,X,nasm->gx,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 225 ierr = VecScatterEnd(nasm->gorestriction,X,nasm->gx,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 226 if (B) { 227 ierr = VecScatterBegin(nasm->gorestriction,B,nasm->gb,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 228 ierr = VecScatterEnd(nasm->gorestriction,B,nasm->gb,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 229 } else { 230 } 231 } 232 */ 233 } 234 for(i=0;i<nasm->n;i++) { 235 Xl = nasm->x[i]; 236 if (B) { 237 Bl = nasm->b[i]; 238 } else { 239 Bl = PETSC_NULL; 240 } 241 ierr = SNESSolve(nasm->subsnes[i],Bl,Xl);CHKERRQ(ierr); 242 } 243 244 if (nasm->usesdm) { 245 ierr = DMLocalToGlobalBegin(dm,Xl,INSERT_VALUES,X);CHKERRQ(ierr); 246 ierr = DMLocalToGlobalEnd(dm,Xl,INSERT_VALUES,X);CHKERRQ(ierr); 247 } else { 248 /* 249 ierr = VecScatterBegin(nasm->girestriction,nasm->gx,X,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 250 ierr = VecScatterEnd(nasm->girestriction,nasm->gx,X,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 251 */ 252 } 253 PetscFunctionReturn(0); 254 } 255 256 #undef __FUNCT__ 257 #define __FUNCT__ "SNESSolve_NASM" 258 PetscErrorCode SNESSolve_NASM(SNES snes) 259 { 260 Vec F; 261 Vec X; 262 Vec B; 263 PetscInt i; 264 PetscReal fnorm; 265 PetscErrorCode ierr; 266 SNESNormType normtype; 267 268 PetscFunctionBegin; 269 270 X = snes->vec_sol; 271 F = snes->vec_func; 272 B = snes->vec_rhs; 273 274 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 275 snes->iter = 0; 276 snes->norm = 0.; 277 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 278 snes->reason = SNES_CONVERGED_ITERATING; 279 ierr = SNESGetNormType(snes, &normtype);CHKERRQ(ierr); 280 if (normtype == SNES_NORM_FUNCTION || normtype == SNES_NORM_INITIAL_ONLY || normtype == SNES_NORM_INITIAL_FINAL_ONLY) { 281 /* compute the initial function and preconditioned update delX */ 282 if (!snes->vec_func_init_set) { 283 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 284 if (snes->domainerror) { 285 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 286 PetscFunctionReturn(0); 287 } 288 } else { 289 snes->vec_func_init_set = PETSC_FALSE; 290 } 291 292 /* convergence test */ 293 if (!snes->norm_init_set) { 294 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 295 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 296 } else { 297 fnorm = snes->norm_init; 298 snes->norm_init_set = PETSC_FALSE; 299 } 300 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 301 snes->iter = 0; 302 snes->norm = fnorm; 303 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 304 SNESLogConvHistory(snes,snes->norm,0); 305 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 306 307 /* set parameter for default relative tolerance convergence test */ 308 snes->ttol = fnorm*snes->rtol; 309 310 /* test convergence */ 311 ierr = (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 312 if (snes->reason) PetscFunctionReturn(0); 313 } else { 314 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 315 SNESLogConvHistory(snes,snes->norm,0); 316 ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr); 317 } 318 319 /* Call general purpose update function */ 320 if (snes->ops->update) { 321 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 322 } 323 324 for (i = 0; i < snes->max_its; i++) { 325 ierr = SNESNASMSolveLocal_Private(snes,B,X);CHKERRQ(ierr); 326 if (normtype == SNES_NORM_FUNCTION || ((i == snes->max_its - 1) && (normtype == SNES_NORM_INITIAL_FINAL_ONLY || normtype == SNES_NORM_FINAL_ONLY))) { 327 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); 328 if (snes->domainerror) { 329 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 330 PetscFunctionReturn(0); 331 } 332 ierr = VecNorm(F, NORM_2, &fnorm);CHKERRQ(ierr); /* fnorm <- ||F|| */ 333 if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm"); 334 } 335 /* Monitor convergence */ 336 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 337 snes->iter = i+1; 338 snes->norm = fnorm; 339 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 340 SNESLogConvHistory(snes,snes->norm,0); 341 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 342 /* Test for convergence */ 343 if (normtype == SNES_NORM_FUNCTION) ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 344 if (snes->reason) PetscFunctionReturn(0); 345 /* Call general purpose update function */ 346 if (snes->ops->update) { 347 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 348 } 349 } 350 if (normtype == SNES_NORM_FUNCTION) { 351 if (i == snes->max_its) { 352 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);CHKERRQ(ierr); 353 if(!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT; 354 } 355 } else { 356 if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */ 357 } 358 PetscFunctionReturn(0); 359 } 360 361 /*MC 362 SNESNASM - Nonlinear Additive Schwartz 363 364 Level: advanced 365 366 .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types) 367 M*/ 368 369 EXTERN_C_BEGIN 370 #undef __FUNCT__ 371 #define __FUNCT__ "SNESCreate_NASM" 372 PetscErrorCode SNESCreate_NASM(SNES snes) 373 { 374 SNES_NASM *nasm; 375 PetscErrorCode ierr; 376 377 PetscFunctionBegin; 378 379 380 ierr = PetscNewLog(snes, SNES_NASM, &nasm);CHKERRQ(ierr); 381 snes->data = (void*)nasm; 382 383 nasm->n = PETSC_DECIDE; 384 nasm->subsnes = 0; 385 nasm->x = 0; 386 nasm->b = 0; 387 nasm->ois = 0; 388 nasm->iis = 0; 389 390 snes->ops->destroy = SNESDestroy_NASM; 391 snes->ops->setup = SNESSetUp_NASM; 392 snes->ops->setfromoptions = SNESSetFromOptions_NASM; 393 snes->ops->view = SNESView_NASM; 394 snes->ops->solve = SNESSolve_NASM; 395 snes->ops->reset = SNESReset_NASM; 396 397 snes->usesksp = PETSC_FALSE; 398 snes->usespc = PETSC_FALSE; 399 400 if (!snes->tolerancesset) { 401 snes->max_its = 10000; 402 snes->max_funcs = 10000; 403 } 404 405 ierr = PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNASMSetSubdomains_C","SNESNASMSetSubdomains_NASM", 406 SNESNASMSetSubdomains_NASM);CHKERRQ(ierr); 407 408 PetscFunctionReturn(0); 409 } 410 EXTERN_C_END 411