Dissecting the Poisson Codes

Objective: At the end of this tutorial you will understand how a simple linear solver works in PetIGA.

Assumptions: I assume that you are familiar with C and in particular understand pointers.

Detailed Look at Poisson1D.c

In this tutorial we will look at the Poisson codes and explain each portion in more detail. It is my hope that in explaining how the demo codes work, this will help you see how you can use PetIGA to address your own problems. In this demo I will be discussing this demo code in sections. It might be useful to first scan through this code to familiarize yourself with its organization.

In most PetIGA applications, you only need to include petiga.h. Including this file will also automatically include petsc.h.

#include "petiga.h"

The first function we encounter is the System function. This function is the key to understanding how to use PetIGA to solve your own problems. For a linear, steady problem, this function represents the evaluation of the bilinear and linear form at a quadrature point.

At a glance, a few things may strike you as strange particularly if you are not familiar with PETSc. First, you will note the #undef and #define lines before the function declaration. This is a mechanism that PETSc uses such that when its functions fail, the entire call stack is returned to the user. This is typically something that you only get when running the debugger. Having it always return on error reduces the time it takes you to locate errors.

You will also see variable types PetscReal and PetscScalar instead of the expected double. These are PETSc types which for most applications are actually of the type double. However, you can configure PETSc to use, for example, quadruple precision or complex numbers. In this case, using PETSc types means that your code automatically ports to these situations with no additional work for you.

The function has as its arguments a IGAPoint, two PetscScalars and a void pointer. For now, ignore the existence of the void pointer. The IGAPoint p is the quadrature/collocation point. We have built in all the information that you will need to evaluate the bilinear and linear form at this point. The two PetscScalar pointers are for the discretized bilnear and linear form respectively. The variable K can be thought of as a flattened point contribution to the element stiffness matrix. The variable F is the point contribution to the element load vector.

The basis functions are part of the IGAPoint p and so we need to simply assign some pointers to their location. The pointers N0 and N1 represent the 0th and 1st order derivative of the basis functions respectively. The pointers used and their casting is a bit technical and ugly, a consequence of using the C language. However, we emphasize that it is sufficient to copy these lines into your own code and understand what it is the pointers address. If a geometry is used, then the basis is mapped and/or made rational if using NURBS. Otherwise, these are just the B-spline basis functions in the parametric space. Your code does not need to know the difference, it is all handled internally.

Then we simply code a double loop over the local number of basis functions, nen. As a convention, we use a to loop over the test functions and b to loop over the trial functions. A dot product of the basis function derivatives is assembled into the stiffness matrix and the unit body force is assembled into the load vector. The quadrature weights and potential Jacobian of the mapping are applied internally.

#undef  __FUNCT__
#define __FUNCT__ "System"
PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
{
  const PetscReal *N0,(*N1)[1];
  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
  PetscInt a,b,nen=p->nen;
  for (a=0; a<nen; a++) {
    PetscReal Na   = N0[a];
    PetscReal Na_x = N1[a][0];
    for (b=0; b<nen; b++) {
      PetscReal Nb_x = N1[b][0];
      K[a*nen+b] = Na_x * Nb_x;
    }
    F[a] = Na * 1.0;
  }
  return 0;
}

Below this we code the program’s main function. As in all PETSc codes, we intially call an initialization routine. Consider this a rite-of-passage which is required for PETSc’s internals and parallelism. However, another oddity may jump out at you. You will see that we have defined an error code ierr and then after the PETSc call, we postpend a function call CHKERRQ(ierr). If you glance at the rest of the code, you will see this pattern appears at most every line of the code (yes they are ugly, be patient, eventually you won’t notice them). Every function in PETSc/PetIGA returns such a code which is really just an integer whose value corresponds to a error type. The function call at the end is really a macro which will cause the application to fail and report the proper error messages/stack.

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char *argv[]) {

  PetscErrorCode  ierr;
  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);

Now we create an IGA, your one-stop-shop for the IGA-discretization. We set the dimensions of the space as well as the number of degrees of freedom per basis. Basic usuage of an IGA in PetIGA assume that the same function space is used for all test/trial spaces in the problem. This is a constraint that we will later show how to relax. The XXXSetFromOptions calls in PETSc are very important. This call is what populates the list of options that you see when you run the code with the -help option. The final call is an internal setup of many data structures which we need to compute. This is considered a final step of creating the IGA. There are many other ways to initialize an IGA, this way is just the simplest and uses the default discretizations.

IGA iga;
ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
ierr = IGASetDim(iga,1);CHKERRQ(ierr);
ierr = IGASetDof(iga,1);CHKERRQ(ierr);
ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
ierr = IGASetUp(iga);CHKERRQ(ierr);

In this block we set the boundary conditions. Again, there are many ways you can set boundary conditions in PetIGA. Here we will set a single value on all basis functions on the boundaries. The IGABoundary represents a side in 1D, a boundary edge in 2D, and a boundary face in 3D. You get the boundary object from the IGA itself by specifying a direction (which parametric dimension) and then a side (0 or 1 reflecting the min and max of the axis). With the IGABoundary in hand, then you can set a constant value on all associated basis functions of a particular component of the solution. Since this is a scalar problem, we set the 0th component.

IGABoundary bnd;
PetscInt dir=0,side;
PetscScalar value = 1.0;
for (side=0; side<2; side++) {
  ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
  ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
}

Now we setup and compute the linear system of our problem, Ax=b. Both Mat and Vec are PETSc types which encapsulate sparsity and parallel communication. You do not need to understand how these work internally–they are given to you by the IGA. The call to IGASetUserSystem associates the System function we wrote earlier to the discretization. The IGA uses this function in the next call to assemble the linear system.

Mat A;
Vec x,b;
ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);

Then we need a KSP to solve the linear system. This we also get from our IGA. This merely ensures that the KSP and the IGA use the same communicator (a detail of parallelism which you can ignore). Then we need to set the problem operators by KSPSetOperators. The first assignment of A is for the problem operator. The second assignment of A specifies the matrix from which PETSc will create a preconditioner. In this case, and in many cases, we want the preconditioner to be created from original system matrix. We also specify a flag for the KSP that our preconditioner has the same nonzero pattern as the matrix A. As before with the IGA, we add the XXXSetFromOptions call which will allows us to control everything about the KSP from the commandline. Finally, we solve.

KSP ksp;
ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);

This call to file:VecView is optional. It is what will draw a rough plot of the degrees of freedom, as described in Post-processing Results with igakit

ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);

Finally, we destroy the objects that we created and issue a finalize call to PETSc.

 ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
 ierr = MatDestroy(&A);CHKERRQ(ierr);
 ierr = VecDestroy(&x);CHKERRQ(ierr);
 ierr = VecDestroy(&b);CHKERRQ(ierr);
 ierr = IGADestroy(&iga);CHKERRQ(ierr);
 ierr = PetscFinalize();CHKERRQ(ierr);
 return 0;
}

That is all that you need to do, everything else in handled internally in the PetIGA library. We feel that our approach is beneficial because it leads to shorter codes which focus on the physics of the problem. We also emphasize that this code will work in serial or in parallel without the author needing to know anything about parallelism aside from having to initialize objects with communicators.

Higher Dimensions and Dimension Independent

How different is the code when changing dimensions of the problem? There is very little which depends on the problem dimension in the application code. If we compare to PETIGA_DIR/demo/Poisson2d.c, we have a few lines in our System function which are different. The pointers to the shape function gradients are now 2D:

const PetscReal *N0,(*N1)[1];

const PetscReal *N0,(*N1)[2];

and the dot product is now 2D:

K[a*nen+b] = Na_x * Nb_x;

K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;

In main we also have a few lines:

ierr = IGASetDim(iga,1);CHKERRQ(ierr);

ierr = IGASetDim(iga,2);CHKERRQ(ierr);

and we need to set more boundary conditions. In fact, this reflects that a dimension independent code is simple to obtain.