Commits

Lisandro Dalcin  committed df9a332

Fixes and enhancements to Bratu demo

  • Participants
  • Parent commits 2987bbd

Comments (0)

Files changed (2)

File demo/Bratu.c

 } AppCtx;
 
 EXTERN_C_BEGIN
-extern PetscErrorCode Bratu_Function(IGAPoint,const PetscScalar U[],PetscScalar F[],void *);
-extern PetscErrorCode Bratu_Jacobian(IGAPoint,const PetscScalar U[],PetscScalar J[],void *);
+extern PetscErrorCode Bratu_Function(IGAPoint,const PetscScalar U[],PetscScalar F[],void *ctx);
+extern PetscErrorCode Bratu_Jacobian(IGAPoint,const PetscScalar U[],PetscScalar J[],void *ctx);
 EXTERN_C_END
 
 EXTERN_C_BEGIN
                                       PetscScalar *F,void *ctx);
 EXTERN_C_END
 
+#if PETSC_VERSION_LT(3,4,0)
+#define TSSolve(ts,x) TSSolve(ts,x,NULL)
+#endif
+
 #undef __FUNCT__
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
   PetscBool steady = PETSC_TRUE;
   PetscReal lambda = 6.80;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Bratu Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-steady","Steady problem",__FILE__,steady,&steady,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsReal("-lambda","Bratu parameter",__FILE__,lambda,&lambda,NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-steady","Steady problem",__FILE__,steady,&steady,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  PetscInt dir,side;
+  PetscInt dim,dir,side;
   for (dir=0; dir<3; dir++) {
     for (side=0; side<2; side++) {
       IGABoundary bnd;
     }
   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  PetscInt dim;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   if (dim < 1) {ierr = IGASetDim(iga,dim=2);CHKERRQ(ierr);}
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   Vec x;
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+
   if (steady) {
     SNES snes;
     ierr = IGACreateSNES(iga,&snes);CHKERRQ(ierr);
+    if (!iga->collocation) {
+      Mat mat; KSP ksp;
+      ierr = SNESGetJacobian(snes,NULL,&mat,NULL,NULL);CHKERRQ(ierr);
+      ierr = MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+      ierr = MatSetOption(mat,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
+      ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
+      ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
+    }
     ierr = SNESSetTolerances(snes,PETSC_DEFAULT,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
-    ierr = SNESSolve(snes,0,x);CHKERRQ(ierr);
+    ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
     ierr = SNESDestroy(&snes);CHKERRQ(ierr);
   } else {
-    TS ts;
+    TS   ts;
     SNES snes;
     ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
     ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr);
     ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr);
     ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
     ierr = SNESSetTolerances(snes,PETSC_DEFAULT,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
+    if (!iga->collocation) {
+      Mat mat; KSP ksp;
+      ierr = TSGetIJacobian(ts,NULL,&mat,NULL,NULL);CHKERRQ(ierr);
+      ierr = MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+      ierr = MatSetOption(mat,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
+      ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
+      ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
+    }
     ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
-#if PETSC_VERSION_LE(3,3,0)
-    ierr = TSSolve(ts,x,NULL);CHKERRQ(ierr);
-#else
     ierr = TSSolve(ts,x);CHKERRQ(ierr);
-#endif
     ierr = TSDestroy(&ts);CHKERRQ(ierr);
   }
   PetscBool draw = PETSC_FALSE;
   ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
-  if (draw && dim < 3) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+  if (draw&&dim<3) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
   ierr = VecDestroy(&x);CHKERRQ(ierr);
 
   ierr = IGADestroy(&iga);CHKERRQ(ierr);

File demo/BratuFJ.F90

 use PetIGA
 implicit none
 
-type, bind(C) :: IGAUser
+type, bind(C) :: AppCtx
    real(kind=IGA_REAL_KIND) lambda
-end type IGAUser
+end type AppCtx
+
+integer, parameter :: rk = IGA_REAL_KIND
 
 contains
 
   use PetIGA
   implicit none
   type(IGAPoint), intent(in) :: p
-  type(IGAUser),  intent(in) :: ctx
+  type(AppCtx),   intent(in) :: ctx
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(out) :: FF(p%neq)
   real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:), hess_N(:,:,:)
 
   if (p%neq == 1) then
      ! collocation
-     FF(1) = - del2_u - ctx%lambda * exp(real(u))
+     FF(1) = - del2_u - ctx%lambda * exp(real(u,rk))
   else
      ! galerkin
      do a = 1, p%nen
-        FF(a) = dot_product(grad_N(:,a),real(grad_u)) - &
-             N(a) * ctx%lambda * exp(real(u))
+        FF(a) = + dot_product(grad_N(:,a),real(grad_u,rk)) &
+                - N(a) * ctx%lambda * exp(real(u,rk))
      end do
   end if
 
   use PetIGA
   implicit none
   type(IGAPoint), intent(in) :: p
-  type(IGAUser),  intent(in) :: ctx
+  type(AppCtx),   intent(in) :: ctx
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(out) :: JJ(p%nen,p%neq)
   real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:), hess_N(:,:,:)
         do i = 1, p%dim
            del2_N = del2_N + hess_N(i,i,a)
         end do
-        JJ(a,1) = - del2_N - N(a) * ctx%lambda * exp(real(u))
+        JJ(a,1) = - del2_N - N(a) * ctx%lambda * exp(real(u,rk))
      end do
   else
      ! galerkin
      do a = 1, p%nen
         do b = 1, p%nen
-           JJ(b,a) = dot_product(grad_N(:,a),grad_N(:,b)) - &
-                N(a) * N(b) * ctx%lambda * exp(real(u))
+           JJ(b,a) = + dot_product(grad_N(:,a),grad_N(:,b)) &
+                     - N(a) * N(b) * ctx%lambda * exp(real(u,rk))
         end do
      end do
   end if
   use PetIGA
   implicit none
   type(IGAPoint), intent(in) :: p
-  type(IGAUser),  intent(in) :: ctx
+  type(AppCtx),   intent(in) :: ctx
   real   (kind=IGA_REAL_KIND   ), intent(in), value :: dt, shift, t
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: VV(p%nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
   grad_u = IGA_Eval(grad_N,UU)
 
   do a = 1, p%nen
-     FF(a) = N(a) * real(v) + &
-             dot_product(grad_N(:,a),real(grad_u)) - &
-             N(a) * ctx%lambda * exp(real(u))
+     FF(a) = + N(a) * real(v,rk) &
+             + dot_product(grad_N(:,a),real(grad_u,rk)) &
+             - N(a) * ctx%lambda * exp(real(u,rk))
   end do
 
   ierr = 0
   use PetIGA
   implicit none
   type(IGAPoint), intent(in) :: p
-  type(IGAUser),  intent(in) :: ctx
+  type(AppCtx),   intent(in) :: ctx
   real   (kind=IGA_REAL_KIND   ), intent(in), value :: dt, shift, t
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: VV(p%nen)
   scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
 
   do a = 1, p%nen
      do b = 1, p%nen
-        JJ(b,a) = shift * N(a) * N(b) + &
-                  dot_product(grad_N(:,a),grad_N(:,b)) - &
-                  N(a) * N(b) * ctx%lambda  * exp(real(u))
+        JJ(b,a) = + shift * N(a) * N(b) &
+                  + dot_product(grad_N(:,a),grad_N(:,b)) &
+                  - N(a) * N(b) * ctx%lambda  * exp(real(u,rk))
      end do
   end do