Commits

Peter Brune committed 0e32a25 Merge

Merge branch 'prbrune/fas-gscolorsecant'

Comments (0)

Files changed (10)

include/petsc-private/snesimpl.h

 
 PetscErrorCode SNESScaleStep_Private(SNES,Vec,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
 
-PETSC_EXTERN PetscLogEvent SNES_Solve, SNES_LineSearch, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval, SNES_NPCSolve;
+PETSC_EXTERN PetscLogEvent SNES_Solve, SNES_LineSearch, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval, SNES_GSFuncEval, SNES_NPCSolve;
 
 #endif

src/snes/examples/tutorials/makefile

 	   ${DIFF} output/ex19_ngmres_fas.out ex19_ngmres_fas.tmp || printf "${PWD}\nPossible problem with ex19_ngmres_fas, diffs above\n=========================================\n"; \
            ${RM} -f ex19_ngmres_fas.tmp
 
+runex19_ngmres_fas_gssecant: #test ex19 with NGMRES preconditioned by FAS with pointwise GS smoother
+	-@${MPIEXEC} -n 1 ./ex19 -da_refine 3 -snes_monitor_short -snes_type ngmres -npc_snes_type fas \
+        -npc_fas_levels_snes_type gs -npc_fas_levels_snes_max_it 6 -npc_fas_levels_snes_gs_secant -npc_fas_levels_snes_gs_max_it 1 -npc_fas_coarse_snes_max_it 1 \
+        -lidvelocity 100 -grashof 4e4 \
+         > ex19_ngmres_fas_gssecant.tmp 2>&1; \
+	   ${DIFF} output/ex19_ngmres_fas_gssecant.out ex19_ngmres_fas_gssecant.tmp || echo  ${PWD} "\nPossible problem with ex19_ngmres_fas_gssecant, diffs above \n========================================="; \
+           ${RM} -f ex19_ngmres_fas_gssecant.tmp
+
 runex19_ngmres_fas_ms: #test ex19 with NGMRES preconditioned by FAS with multi-stage smoother
 	-@${MPIEXEC} -n 2 ./ex19 -snes_grid_sequence 2 -lidvelocity 200 -grashof 1e4 -snes_monitor_short -snes_view -snes_converged_reason -snes_type ngmres -npc_snes_type fas -npc_fas_coarse_snes_type newtonls -npc_fas_coarse_ksp_type preonly -npc_snes_max_it 1 -npc_fas_snes_type ms -npc_fas_snes_max_it 1 -npc_fas_ksp_type preonly -npc_fas_pc_type none -npc_fas_snes_ms_type m62 -npc_fas_snes_max_it 1 -npc_fas_snes_ms_damping 0.22 > ex19_ngmres_fas_ms.tmp 2>&1; \
 	   ${DIFF} output/ex19_ngmres_fas_ms.out ex19_ngmres_fas_ms.tmp || printf "${PWD}\nPossible problem with ex19_ngmres_fas_ms, diffs above\n=========================================\n"; \

src/snes/examples/tutorials/output/ex19_ngmres_fas_gssecant.out

+lid velocity = 100, prandtl # = 1, grashof # = 40000
+  0 SNES Function norm 1667.67 
+  1 SNES Function norm 227.403 
+  2 SNES Function norm 94.4607 
+  3 SNES Function norm 20.4337 
+  4 SNES Function norm 17.4825 
+  5 SNES Function norm 7.36826 
+  6 SNES Function norm 3.09886 
+  7 SNES Function norm 0.802506 
+  8 SNES Function norm 0.554519 
+  9 SNES Function norm 0.317552 
+ 10 SNES Function norm 0.112429 
+ 11 SNES Function norm 0.0653065 
+ 12 SNES Function norm 0.0269513 
+ 13 SNES Function norm 0.0119604 
+ 14 SNES Function norm 0.00281523 
+ 15 SNES Function norm 0.00123359 
+ 16 SNES Function norm 0.000380515 
+ 17 SNES Function norm 0.000176141 
+ 18 SNES Function norm 4.13561e-05 
+ 19 SNES Function norm 7.31869e-06 
+Number of SNES iterations = 19

src/snes/impls/fas/fas.c

   ierr = MatDestroy(&fas->interpolate);CHKERRQ(ierr);
   ierr = MatDestroy(&fas->restrct);CHKERRQ(ierr);
   ierr = VecDestroy(&fas->rscale);CHKERRQ(ierr);
-  if (fas->next) ierr = SNESReset(fas->next);CHKERRQ(ierr);
+  if (fas->galerkin) {
+    ierr = VecDestroy(&fas->Xg);CHKERRQ(ierr);
+    ierr = VecDestroy(&fas->Fg);CHKERRQ(ierr);
+  }
+  if (fas->next) {ierr = SNESReset(fas->next);CHKERRQ(ierr);}
   PetscFunctionReturn(0);
 }
 

src/snes/impls/gs/gsimpl.h

+#if !defined(__GSIMPL)
+#define __GSIMPL
+
+#include <petsc-private/snesimpl.h>      /*I "petscsnes.h"  I*/
+#include <petscdm.h>
+
+typedef struct {
+  PetscInt  sweeps;     /* number of sweeps through the local subdomain before neighbor communication */
+  PetscInt  max_its;    /* maximum iterations of the inner pointblock solver */
+  PetscReal rtol;       /* relative tolerance of the inner pointblock solver */
+  PetscReal abstol;     /* absolute tolerance of the inner pointblock solver */
+  PetscReal stol;       /* step tolerance of the inner pointblock solver */
+  PetscReal h;          /* differencing for secant variants */
+  PetscBool secant_mat; /* use the Jacobian to get the coloring for the secant */
+} SNES_GS;
+
+PETSC_EXTERN PetscErrorCode SNESComputeGSDefaultSecant(SNES,Vec,Vec,void *);
+
+#endif

src/snes/impls/gs/gssecant.c

+#include <../src/snes/impls/gs/gsimpl.h>
+
+#undef __FUNCT__
+#define __FUNCT__ "GSDestroy_Private"
+PetscErrorCode GSDestroy_Private(ISColoring coloring)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "SNESComputeGSDefaultSecant"
+PETSC_EXTERN PetscErrorCode SNESComputeGSDefaultSecant(SNES snes,Vec X,Vec B,void *ctx)
+{
+  PetscErrorCode ierr;
+  SNES_GS *gs = (SNES_GS*)snes->data;
+  PetscInt       i,j,k,ncolors;
+  DM             dm;
+  PetscBool      flg;
+  ISColoring     coloring;
+  MatColoring    mc;
+  Vec            W,G,F;
+  PetscScalar    h=gs->h;
+  IS             *coloris;
+  PetscScalar    f,g,x,w,d;
+  PetscReal      dxt,xt,ft,ft1=0;
+  const PetscInt *idx;
+  PetscInt       size,s;
+  PetscReal      atol,rtol,stol;
+  PetscInt       its;
+  PetscErrorCode (*func)(SNES,Vec,Vec,void*);
+  void           *fctx;
+  PetscContainer colorcontainer;
+  PetscBool      mat = gs->secant_mat,equal,isdone,alldone;
+  PetscScalar    *xa,*fa,*wa,*ga;
+
+  PetscFunctionBegin;
+  if (snes->nwork < 3) {
+    ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr);
+  }
+  W = snes->work[0];
+  G = snes->work[1];
+  F = snes->work[2];
+  ierr = VecGetOwnershipRange(X,&s,NULL);CHKERRQ(ierr);
+  ierr = SNESGSGetTolerances(snes,&atol,&rtol,&stol,&its);CHKERRQ(ierr);
+  ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr);
+  ierr = SNESGetFunction(snes,NULL,&func,&fctx);CHKERRQ(ierr);
+  ierr = PetscObjectQuery((PetscObject)snes,"SNESGSColoring",(PetscObject*)&colorcontainer);CHKERRQ(ierr);
+  if (!colorcontainer) {
+    /* create the coloring */
+    ierr = DMHasColoring(dm,&flg);CHKERRQ(ierr);
+    if (flg && !mat) {
+      ierr = DMCreateColoring(dm,IS_COLORING_GLOBAL,&coloring);CHKERRQ(ierr);
+    } else {
+      if (!snes->jacobian) {ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr);}
+      ierr = MatColoringCreate(snes->jacobian,&mc);CHKERRQ(ierr);
+      ierr = MatColoringSetDistance(mc,1);CHKERRQ(ierr);
+      ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
+      ierr = MatColoringApply(mc,&coloring);CHKERRQ(ierr);
+      ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
+    }
+    ierr = PetscContainerCreate(PetscObjectComm((PetscObject)snes),&colorcontainer);CHKERRQ(ierr);
+    ierr = PetscContainerSetPointer(colorcontainer,(void *)coloring);CHKERRQ(ierr);
+    ierr = PetscContainerSetUserDestroy(colorcontainer,(PetscErrorCode (*)(void *))GSDestroy_Private);CHKERRQ(ierr);
+    ierr = PetscObjectCompose((PetscObject)snes,"SNESGSColoring",(PetscObject)colorcontainer);CHKERRQ(ierr);
+    ierr = PetscContainerDestroy(&colorcontainer);CHKERRQ(ierr);
+  } else {
+    ierr = PetscContainerGetPointer(colorcontainer,(void **)&coloring);CHKERRQ(ierr);
+  }
+  ierr = ISColoringGetIS(coloring,&ncolors,&coloris);CHKERRQ(ierr);
+  ierr = VecEqual(X,snes->vec_sol,&equal);CHKERRQ(ierr);
+  if (equal && snes->normschedule == SNES_NORM_ALWAYS) {
+    /* assume that the function is already computed */
+    ierr = VecCopy(snes->vec_func,F);CHKERRQ(ierr);
+  } else {
+    ierr = PetscLogEventBegin(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr);
+    ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr);
+    ierr = PetscLogEventEnd(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr);
+    if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);}
+  }
+  ierr = VecGetArray(X,&xa);CHKERRQ(ierr);
+  ierr = VecGetArray(F,&fa);CHKERRQ(ierr);
+  ierr = VecGetArray(G,&ga);CHKERRQ(ierr);
+  ierr = VecGetArray(W,&wa);CHKERRQ(ierr);
+  for (i=0;i<ncolors;i++) {
+    ierr = ISGetIndices(coloris[i],&idx);CHKERRQ(ierr);
+    ierr = ISGetLocalSize(coloris[i],&size);CHKERRQ(ierr);
+    ierr = VecCopy(X,W);CHKERRQ(ierr);
+    for (j=0;j<size;j++) {
+      wa[idx[j]-s] += h;
+    }
+    ierr = PetscLogEventBegin(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr);
+    ierr = (*func)(snes,W,G,fctx);CHKERRQ(ierr);
+    ierr = PetscLogEventEnd(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr);
+    if (B) {ierr = VecAXPY(G,-1.0,B);CHKERRQ(ierr);}
+    for (k=0;k<its;k++) {
+      dxt = 0.;
+      xt = 0.;
+      ft = 0.;
+      for (j=0;j<size;j++) {
+        f = fa[idx[j]-s];
+        x = xa[idx[j]-s];
+        g = ga[idx[j]-s];
+        w = wa[idx[j]-s];
+        if (PetscAbsScalar(g-f) > atol) {
+          d = (x*g-w*f) / PetscRealPart(g-f);
+        } else {
+          d = x;
+        }
+        dxt += PetscRealPart(PetscSqr(d-x));
+        xt += PetscRealPart(PetscSqr(x));
+        ft += PetscRealPart(PetscSqr(f));
+        xa[idx[j]-s] = d;
+      }
+
+      if (k == 0) ft1 = PetscSqrtReal(ft);
+      if (k<its-1) {
+        isdone = PETSC_FALSE;
+        if (stol*PetscSqrtReal(xt) > PetscSqrtReal(dxt)) isdone = PETSC_TRUE;
+        if (PetscSqrtReal(ft) < atol) isdone = PETSC_TRUE;
+        if (rtol*ft1 > PetscSqrtReal(ft)) isdone = PETSC_TRUE;
+        ierr = MPI_Allreduce(&isdone,&alldone,1,MPIU_BOOL,MPI_BAND,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
+        if (alldone) break;
+      }
+      if (i < ncolors-1 || k < its-1) {
+        ierr = PetscLogEventBegin(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr);
+        ierr = (*func)(snes,X,F,fctx);CHKERRQ(ierr);
+        ierr = PetscLogEventEnd(SNES_GSFuncEval,snes,X,B,0);CHKERRQ(ierr);
+        if (B) {ierr = VecAXPY(F,-1.0,B);CHKERRQ(ierr);}
+      }
+      if (k<its-1) {
+        ierr = VecSwap(X,W);CHKERRQ(ierr);
+        ierr = VecSwap(F,G);CHKERRQ(ierr);
+      }
+    }
+  }
+  ierr = VecRestoreArray(X,&xa);CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,&fa);CHKERRQ(ierr);
+  ierr = VecRestoreArray(G,&ga);CHKERRQ(ierr);
+  ierr = VecRestoreArray(W,&wa);CHKERRQ(ierr);
+  ierr = ISColoringRestoreIS(coloring,&coloris);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}

src/snes/impls/gs/makefile

 
 CFLAGS   =
 FFLAGS   =
-SOURCEC  = snesgs.c
+SOURCEC  = snesgs.c gssecant.c
 SOURCEF  =
 SOURCEH  =
 LIBBASE  = libpetscsnes

src/snes/impls/gs/snesgs.c

-#include <petsc-private/snesimpl.h>             /*I   "petscsnes.h"   I*/
-
-typedef struct {
-  PetscInt  sweeps;     /* number of sweeps through the local subdomain before neighbor communication */
-  PetscInt  max_its;    /* maximum iterations of the inner pointblock solver */
-  PetscReal rtol;       /* relative tolerance of the inner pointblock solver */
-  PetscReal abstol;     /* absolute tolerance of the inner pointblock solver */
-  PetscReal stol;       /* step tolerance of the inner pointblock solver */
-} SNES_GS;
-
+#include <../src/snes/impls/gs/gsimpl.h>      /*I "petscsnes.h"  I*/
 
 #undef __FUNCT__
 #define __FUNCT__ "SNESGSSetTolerances"
 #define __FUNCT__ "SNESSetUp_GS"
 PetscErrorCode SNESSetUp_GS(SNES snes)
 {
+  PetscErrorCode ierr;
+  PetscErrorCode (*f)(SNES,Vec,Vec,void*);
+
   PetscFunctionBegin;
+  ierr = SNESGetGS(snes,&f,NULL);CHKERRQ(ierr);
+  if (!f) {
+    ierr = SNESSetGS(snes,SNESComputeGSDefaultSecant,NULL);CHKERRQ(ierr);
+  }
   PetscFunctionReturn(0);
 }
 
   if (flg || flg1 || flg2 || flg3) {
     ierr = SNESGSSetTolerances(snes,atol,rtol,stol,max_its);CHKERRQ(ierr);
   }
+  flg  = PETSC_FALSE;
+  ierr = PetscOptionsBool("-snes_gs_secant","Use pointwise secant local Jacobian approximation","",flg,&flg,NULL);CHKERRQ(ierr);
+  if (flg) {
+    ierr = SNESSetGS(snes,SNESComputeGSDefaultSecant,NULL);CHKERRQ(ierr);
+    ierr = PetscInfo(snes,"Setting default finite difference coloring Jacobian matrix\n");CHKERRQ(ierr);
+  }
+  ierr = PetscOptionsReal("-snes_gs_secant_h","Differencing parameter for secant search","",gs->h,&gs->h,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-snes_gs_secant_mat_coloring","Use the Jacobian coloring for the secant GS","",gs->secant_mat,&gs->secant_mat,&flg);CHKERRQ(ierr);
+
   ierr = PetscOptionsTail();CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
     ierr = SNESMonitor(snes,0,snes->norm);CHKERRQ(ierr);
   }
 
-
   /* Call general purpose update function */
   if (snes->ops->update) {
     ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr);
   gs->abstol  = 1e-15;
   gs->stol    = 1e-12;
   gs->max_its = 50;
+  gs->h       = 1e-8;
 
   snes->data = (void*) gs;
   PetscFunctionReturn(0);

src/snes/interface/dlregissnes.c

   ierr = PetscLogEventRegister("SNESSolve",            SNES_CLASSID,&SNES_Solve);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("SNESFunctionEval",     SNES_CLASSID,&SNES_FunctionEval);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("SNESGSEval",           SNES_CLASSID,&SNES_GSEval);CHKERRQ(ierr);
+  ierr = PetscLogEventRegister("SNESGSFuncEval",       SNES_CLASSID,&SNES_GSFuncEval);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("SNESJacobianEval",     SNES_CLASSID,&SNES_JacobianEval);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("SNESLineSearch",       SNESLINESEARCH_CLASSID,&SNESLineSearch_Apply);CHKERRQ(ierr);
   ierr = PetscLogEventRegister("SNESNPCSolve",         SNES_CLASSID,&SNES_NPCSolve);CHKERRQ(ierr);

src/snes/interface/snes.c

 
 /* Logging support */
 PetscClassId  SNES_CLASSID, DMSNES_CLASSID;
-PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval, SNES_NPCSolve;
+PetscLogEvent SNES_Solve, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval, SNES_GSFuncEval, SNES_NPCSolve;
 
 #undef __FUNCT__
 #define __FUNCT__ "SNESSetErrorIfNotConverged"