User_guide - apusok/FD-PDE GitHub Wiki

User guide

Manual and benchmarks

Documentation for the FD-PDE Framework and tests can be found here: FD-PDE Benchmarks.

Note: the directory tests/ may contain additional tests (we try our best to keep the manual updated).

FD-PDE Objects

The following FD-PDE objects are currently implemented

These be combined to solve for coupled systems of equations.

Introduction to FD-PDE

For any FD-PDE object (set with FDPDEType), we separate the unknowns (internal discretization) from the coefficients and boundary conditions (user defined). The simplest usage of the FD-PDE framework is in ./tests/test_fdpde.c:

FDPDE           fd;
  
PetscCall(FDPDECreate(PETSC_COMM_WORLD,nx,nz,0.0,1.0,0.0,1.0,FDPDE_STOKES,&fd));
PetscCall(FDPDESetUp(fd));

/* In between FDPDESetUp()/FDPDESolve(), the user may set: 
   - variable grid coordinates
   - form coefficient evaluation function 
   - BCs evaluation function 
   - SNES options
   - other model options */

// PetscCall(FDPDESolve(fd,NULL));
PetscCall(FDPDEView(fd)); 
PetscCall(FDPDEDestroy(&fd));

In the code above, we created an FD-PDE object fd of type FDPDE_STOKES on a 2-D rectangular grid of coordinates [0,1] x [0,1] that will be discretised in nx x nz cells. Before solving the system of equations contained inside fd, one must call FDPDESetUp(fd). After use, the object is destroyed amd memory is freed.

PetscCall() is the PETSc error handling functionality and good practice to use it. It will catch any warnings and errors in the funcions inside. Try to fix warnings as well, as they are future errors.

Structures inside FD-PDE objects

  • dm/x - a DMStag object to describe the discretization of unknowns, and its associated vector. Called with
FDPDEGetDM(fd,&dm);
FDPDEGetSolution(fd,&x);

Both must be destroyed after usage with

VecDestroy(&x);
DMDestroy(&dm);
  • dmcoeff/xcoeff - a DMStag object to describe the coefficients entering the discretization, and its associated vector. xcoeff must be set in FormCoefficient (see below). Both dmcoeff and xcoeff can be called outside FormCoefficient with
FDPDEGetCoefficient(fd,&dmcoeff,&xcoeff)

Note: This routine only gives pointers to dmcoeff/xcoeff so no need to destroy them after usage. They are freed during FDPDEDestroy(&fd).

  • bclist - a DMStagBCList object to describe the boundary conditions. bclist must be set in FormBCList (see below).

FormCoefficient

One can set the coefficient evaluation function with

// Set coefficients evaluation function
FDPDESetFunctionCoefficient(fd,FormCoefficient,coeff_description,usr);

The user will have to create a coefficient function of the form (mimics the FormFunction for residual evaluation required by SNES):

PetscErrorCode FormCoefficient(DM dm, Vec x, DM dmcoeff, Vec xcoeff, void *ctx)

where dm/x are the DM object and associated vector of the solution, and dmcoeff/coeff the coefficients. The user may pass a user context through the ctx argument, and a description. This is one of the most important functions that contains the physics of the problem!

FormBCList

One can set the boundary conditions evaluation function with

// Set BC evaluation function
FDPDESetFunctionBCList(fd,FormBCList,bc_description,usr);

The user will have to create a boundary conditions function to be evaluated (mimics the FormFunction for residual evaluation required by SNES):

PetscErrorCode FormBCList(DM dm, Vec x, DMStagBCList bclist, void *ctx)

where dm/x are the DM object and associated vector of the solution, and bclist is an object containing the boundary conditions dofs. The user may pass a user context through the ctx argument and a description. This is one of the most important functions that contains the boundary conditions of the problem!

How to set boundary conditions:

DMStagBCListGetValues(bclist,boundary_id,component_id,dof_id,&n,&idx,&x_bc,&value,&type);
for (k=0; k<n_bc; k++) {
    value_bc[k] = 1.0;
    type_bc[k] = BC_DIRICHLET;
}
DMStagBCListInsertValues(bclist,component_id,dof_id,&n,&idx,&x,&value,&type);

where:

  • boundary_id - the boundary label: left - 'w', right - 'e', down - 's', up - 'n',
  • component_id - the dof label: '.' vertex, '-' edge (horizontal), '|' edge (vertical), 'o' element
  • dof_id - component degree of freedom (DMStagStencil c)
  • n - count of boundary dof
  • idx - 1-D array containing the index
  • x_bc - 1-D array containing the coordinates
  • value - 1-D array containing the value
  • type - 1-D array containing BCtype: BC_NULL = 0, BC_DIRICHLET, BC_DIRICHLET_STAG, BC_NEUMANN, BC_NEUMANN_T, BC_PERIODIC

Variable grid spacing

Example in ./tests/test_stokes_solcx_vargrid.c

SOLVE

Once the FD-PDE object has been set-up (domain size and FDPDE type), and the FormCoefficient() and FormBCList() have been provided, the system of equations can be solved. The FDPDESolve() is intended for the user to have a simple way to solve the system. However, the user also has the freedom to modify solver options for more difficult problems, or even use SNES/PETSc routines instead. Here are a number of options:

Option 1 - use default FD-PDE solve, with initial guess=0, set command line snes options

FDPDESolve(fd,NULL);

Option 2 - use default FD-PDE solve, but set snes options with prefix

SNES snes;
FDPDEGetSNES(fd,&snes);
SNESSetOptionsPrefix(snes,"pv_");
FDPDESolve(fd,NULL);

Option 3 - use default FD-PDE solve, but set snes options with prefix (version 2)

SNESSetFromOptions(fd->snes);
SNESSetOptionsPrefix(fd->snes,"pv_");
FDPDESolve(fd,NULL);

Option 4 - change initial guess with

FDPDEGetSolutionGuess(fd,&xguess);
// change initial guess here <> before FDPDESolve()
FDPDESolve(fd,NULL);

Option 5 - example structure for time-dependent problems

-> initial guess/conditions
istep = 0;
// Time loop
while (istep <= nstep) {
  FDPDESolve(fd,NULL);
  FDPDEGetSolution(fd,&x);
  -> copy solution to old_solution
  -> may/not change initial guess for next time step; by default xguess = x after successful FDPDESolve()
}

Option 6 - let user have full control (skip FDPDESolve, use SNES routines). Advanced usage.

FDPDEGetSNES(fd,&snes);
-> set snes options
SNESSetFromOptions();
-> define xguess
SNESSolve(snes,NULL,x);
-> control outcome of SNES

Visualization

PetscBinary/Python

The default DMStag output routine is DMStagViewBinaryPython(). The majority of test files have various examples how to use them, and the python scripts in tests/python/ show how to load and visualize the data.

However, a very simple and general way to quickly inspect any DMStag solution written to file using PetscBinary/python, is using:

1. Imshow

For equidistant grids (fastest plotting version):

python
>>> import dmstagoutput as dmout
>>> dmout.general_output_imshow('out_solcx',None,'bilinear')

General usage:

general_output_imshow(fname,cmap_input,interp_input)

where:

More info in utils/dmstagoutput.py. These routines will create a .pdf for every dof field in a DMStag and associated vector. Example test for petsc/python binary output:

./tests/test_dmstagoutput_

2. Pcolormesh

For variable grids, fast version of pcolor() but cannot use interpolation straightforward (user needs to manipulate data):

python
>>> import dmstagoutput as dmout
>>> general_output_pcolormesh('out_solcx','RdBu')

General usage:

general_output_pcolormesh(fname,cmap_input)

where fname, cmap_input, interp_input are explained above.

VTK/Paraview

Tests test_stokes_solcx_vargrid.c and test_stokes_rt.c contain examples how to output data as VTK/XMF files that can be visualised using Paraview.

Back to HOME