I know this question is not related to fenics, but I really have NO way to get any help 
…
I tried to convert a C code of PETsc to python version, but the python code cannot converge. I have checked many times but I don’t know where I am wrong. I will really appreciate if anyone can help me on this.
Here is the C code, from the book “PETSc for Partial Differential Equations”, and the python code, using the petsc4py.
static char help[] = "";
#include <petsc.h>
#include <stdio.h>
typedef struct {
PetscReal D0; // conductivity
} HeatCtx;
static PetscReal f_source(PetscReal x, PetscReal y) {
return 3.0 * PetscExpReal(-25.0 * (x-0.6) * (x-0.6))
* PetscSinReal(2.0*PETSC_PI*y);
}
static PetscReal gamma_neumann(PetscReal y) {
return PetscSinReal(6.0 * PETSC_PI * y);
}
extern PetscErrorCode Spacings(DMDALocalInfo*, PetscReal*, PetscReal*);
extern PetscErrorCode EnergyMonitor(TS, PetscInt, PetscReal, Vec, void*);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, PetscReal, PetscReal**,
PetscReal**, HeatCtx*);
extern PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo*, PetscReal, PetscReal**,
Mat, Mat, HeatCtx*);
int main(int argc,char **argv) {
HeatCtx user;
TS ts;
Vec u;
DM da;
DMDALocalInfo info;
PetscReal t0, tf;
PetscBool monitorenergy = PETSC_FALSE;
PetscCall(PetscInitialize(&argc,&argv,NULL,help));
user.D0 = 1.0;
PetscOptionsBegin(PETSC_COMM_WORLD, "ht_", "options for heat", "");
PetscCall(PetscOptionsReal("-D0","constant thermal diffusivity",
"heat.c",user.D0,&user.D0,NULL));
PetscCall(PetscOptionsBool("-monitor","also display total heat energy at each step",
"heat.c",monitorenergy,&monitorenergy,NULL));
PetscOptionsEnd();
//STARTDMDASETUP
PetscCall(DMDACreate2d(PETSC_COMM_WORLD,
DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR,
5, 4, PETSC_DECIDE,PETSC_DECIDE, // default to hx=hx=0.25 grid
1, 1, // degrees of freedom, stencil width
NULL,NULL,&da));
PetscCall(DMSetFromOptions(da));
PetscCall(DMSetUp(da));
PetscCall(DMCreateGlobalVector(da,&u));
//ENDDMDASETUP
//STARTTSSETUP
PetscCall(TSCreate(PETSC_COMM_WORLD,&ts));
PetscCall(TSSetProblemType(ts,TS_NONLINEAR));
PetscCall(TSSetDM(ts,da));
PetscCall(TSSetApplicationContext(ts,&user));
PetscCall(DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
(DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user));
PetscCall(DMDATSSetRHSJacobianLocal(da,
(DMDATSRHSJacobianLocal)FormRHSJacobianLocal,&user));
if (monitorenergy) {
PetscCall(TSMonitorSet(ts,EnergyMonitor,&user,NULL));
}
PetscCall(TSSetType(ts,TSBDF));
PetscCall(TSSetTime(ts,0.0));
PetscCall(TSSetMaxTime(ts,0.1));
PetscCall(TSSetTimeStep(ts,0.001));
PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP));
PetscCall(TSSetFromOptions(ts));
//ENDTSSETUP
// report on set up
PetscCall(TSGetTime(ts,&t0));
PetscCall(TSGetMaxTime(ts,&tf));
PetscCall(DMDAGetLocalInfo(da,&info));
PetscCall(PetscPrintf(PETSC_COMM_WORLD,
"solving on %d x %d grid for t0=%g to tf=%g ...\n",
info.mx,info.my,t0,tf));
// solve
PetscCall(VecSet(u,0.0)); // initial condition
PetscCall(TSSolve(ts,u));
PetscCall(VecDestroy(&u));
PetscCall(TSDestroy(&ts));
PetscCall(DMDestroy(&da));
PetscCall(PetscFinalize());
return 0;
}
PetscErrorCode Spacings(DMDALocalInfo *info, PetscReal *hx, PetscReal *hy) {
if (hx) *hx = 1.0 / (PetscReal)(info->mx-1);
if (hy) *hy = 1.0 / (PetscReal)(info->my); // periodic direction
return 0;
}
//STARTMONITOR
PetscErrorCode EnergyMonitor(TS ts, PetscInt step, PetscReal time, Vec u,
void *ctx) {
HeatCtx *user = (HeatCtx*)ctx;
PetscReal lenergy = 0.0, energy, dt, hx, hy, **au;
PetscInt i,j;
MPI_Comm com;
DM da;
DMDALocalInfo info;
PetscCall(TSGetDM(ts,&da));
PetscCall(DMDAGetLocalInfo(da,&info));
PetscCall(DMDAVecGetArrayRead(da,u,&au));
for (j = info.ys; j < info.ys + info.ym; j++) {
for (i = info.xs; i < info.xs + info.xm; i++) {
if ((i == 0) || (i == info.mx-1))
lenergy += 0.5 * au[j][i];
else
lenergy += au[j][i];
}
}
PetscCall(DMDAVecRestoreArrayRead(da,u,&au));
PetscCall(Spacings(&info,&hx,&hy));
lenergy *= hx * hy;
PetscCall(PetscObjectGetComm((PetscObject)(da),&com));
PetscCall(MPI_Allreduce(&lenergy,&energy,1,MPIU_REAL,MPIU_SUM,com));
PetscCall(TSGetTimeStep(ts,&dt));
//int rank;
//MPI_Comm_rank(MPI_COMM_WORLD, &rank);
//printf("rank = %d, dt = %8.4f, hx = %8.4f, dy = %8.4f\n", rank, dt, hx, hy);
PetscCall(PetscPrintf(PETSC_COMM_WORLD," energy = %9.2e nu = %8.4f\n",
energy,user->D0*dt/(hx*hy)));
return 0;
}
//ENDMONITOR
//STARTRHSFUNCTION
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info,
PetscReal t, PetscReal **au,
PetscReal **aG, HeatCtx *user) {
PetscInt i, j, mx = info->mx;
PetscReal hx, hy, x, y, ul, ur, uxx, uyy;
PetscCall(Spacings(info,&hx,&hy));
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
for (j = info->ys; j < info->ys + info->ym; j++) {
y = hy * j;
for (i = info->xs; i < info->xs + info->xm; i++) {
x = hx * i;
// apply Neumann b.c.s
ul = (i == 0) ? au[j][i+1] + 2.0 * hx * gamma_neumann(y)
: au[j][i-1];
ur = (i == mx-1) ? au[j][i-1] : au[j][i+1];
uxx = (ul - 2.0 * au[j][i]+ ur) / (hx*hx);
// DMDA is periodic in y
uyy = (au[j-1][i] - 2.0 * au[j][i]+ au[j+1][i]) / (hy*hy);
aG[j][i] = user->D0 * (uxx + uyy) + f_source(x,y);
}
}
// if (rank == 0){
// printf("%20.12f\n", aG[2][1]);
// }
return 0;
}
//ENDRHSFUNCTION
//STARTRHSJACOBIAN
PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo *info,
PetscReal t, PetscReal **au,
Mat J, Mat P, HeatCtx *user) {
PetscInt i, j, ncols;
const PetscReal D = user->D0;
PetscReal hx, hy, hx2, hy2, v[5];
MatStencil col[5],row;
PetscCall(Spacings(info,&hx,&hy));
hx2 = hx * hx; hy2 = hy * hy;
for (j = info->ys; j < info->ys+info->ym; j++) {
row.j = j; col[0].j = j;
for (i = info->xs; i < info->xs+info->xm; i++) {
// set up a standard 5-point stencil for the row
row.i = i;
col[0].i = i;
v[0] = - 2.0 * D * (1.0 / hx2 + 1.0 / hy2);
col[1].j = j-1; col[1].i = i; v[1] = D / hy2;
col[2].j = j+1; col[2].i = i; v[2] = D / hy2;
col[3].j = j; col[3].i = i-1; v[3] = D / hx2;
col[4].j = j; col[4].i = i+1; v[4] = D / hx2;
ncols = 5;
// if at the boundary, edit the row back to 4 nonzeros
if (i == 0) {
ncols = 4;
col[3].j = j; col[3].i = i+1; v[3] = 2.0 * D / hx2;
} else if (i == info->mx-1) {
ncols = 4;
col[3].j = j; col[3].i = i-1; v[3] = 2.0 * D / hx2;
}
PetscCall(MatSetValuesStencil(P,1,&row,ncols,col,v,INSERT_VALUES));
}
}
PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY));
if (J != P) {
PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY));
}
return 0;
}
//ENDRHSJACOBIAN
Python version
# -*- coding:utf-8 -*-
import math
import sys
import petsc4py
petsc4py.init(sys.argv)
from petsc4py import PETSc
from mpi4py import MPI
class Head2D:
def __init__(self, da, D0):
self.da = da
self.D0 = D0
self.rank = 0
self.local_vec_u = da.createLocalVec()
#Callable[[TS, float, Vec, Vec], None]
def FormRHSFunctionLocal(self, ts, t, vec_u, vec_g):
self.da.globalToLocal(vec_u, self.local_vec_u)
au = self.da.getVecArray(self.local_vec_u)
aG = self.da.getVecArray(vec_g)
mx, my = self.da.getSizes()
hx = 1.0 / (mx - 1)
hy = 1.0 / my
(xs, xe), (ys, ye) = self.da.getRanges()
for j in range(ys, ye):
y = hy * j
for i in range(xs, xe):
x = hx * i
ul = au[i+1, j] + 2.0*hx*gamma_neumann(y) if i==0 else au[i-1, j]
ur = au[i-1, j] if i == mx-1 else au[i+1, j]
uxx = (ul - 2.0*au[i, j] + ur)/(hx*hx)
uyy = (au[i, j-1] - 2.0*au[i, j] + au[i, j+1])/(hy*hy)
aG[i, j] = self.D0*(uxx + uyy) + f_source(x,y)
#aG[i, j] = -aG[i, j]
# Callable[[TS, float, Vec, Mat, Mat], None]
def FormRHSJacobianLocal(self, ts, t, vec_u, Mat_J, Mat_P):
row = PETSc.Mat.Stencil()
col = [PETSc.Mat.Stencil()] * 5
v = [0] * 5
mx, my = self.da.getSizes()
hx = 1.0 / (mx - 1)
hy = 1.0 / my
hx2 = hx * hx
hy2 = hy * hy
D = self.D0
(xs, xe), (ys, ye) = self.da.getRanges()
for j in range(ys, ye):
row.j = j
col[0].j = j
for i in range(xs, xe):
row.i = i
col[0].i = i
v[0] = - 2.0 * D * (1.0 / hx2 + 1.0 / hy2)
col[1].j = j-1; col[1].i = i; v[1] = D / hy2
col[2].j = j+1; col[2].i = i; v[2] = D / hy2
col[3].j = j; col[3].i = i-1; v[3] = D / hx2
col[4].j = j; col[4].i = i+1; v[4] = D / hx2
ncols = 5
if i == 0:
ncols = 4
col[3].j = j; col[3].i = i+1; v[3] = 2.0 * D / hx2
elif i == mx-1:
ncols = 4
col[3].j = j; col[3].i = i-1; v[3] = 2.0 * D / hx2
for ic in range(ncols):
Mat_P.setValueStencil(row, col[ic], v[ic])
Mat_P.assemble()
if Mat_J != Mat_P:
Mat_J.assemble()
return True
#alias of Callable[[TS, int, float, Vec], None]
def EnergyMonitor(self, ts, step, time, u):
self.da.globalToLocal(u, self.local_vec_u)
au = self.da.getVecArray(self.local_vec_u)
mx, my = self.da.getSizes()
lenergy = 0
(xs, xe), (ys, ye) = self.da.getRanges()
for j in range(ys, ye):
for i in range(xs, xe):
if (i == 0 or i == mx-1):
lenergy += 0.5*au[i, j]
else:
lenergy += au[i, j]
hx = 1.0 / (mx-1)
hy = 1.0 / my
lenergy *= (hx*hy)
energy = da.comm.tompi4py().allreduce(lenergy, MPI.SUM)
dt = ts.getTimeStep()
PETSc.Sys.Print(" energy = %9.2e nu = %8.4f"%(energy, self.D0*dt/(hx*hy)))
def f_source(x, y):
return 3.0 * math.exp(-25.0*(x-0.6)*(x-0.6)) * math.sin(2.0*math.pi*y)
def gamma_neumann(y):
return math.sin(6.0*math.pi*y)
OptDB = PETSc.Options("ht_")
D0 = OptDB.getReal("D0", 1.0)
monitor = OptDB.getBool("monitor", False)
mx = OptDB.getInt("mx", 5)
my = OptDB.getInt("my", 4)
# You can also :
da = PETSc.DMDA()
da.create(comm = PETSc.COMM_WORLD,
dim=2,
sizes=(mx, my),
proc_sizes=None, #PETSC_DECIDE..
boundary_type=(PETSc.DM.BoundaryType.NONE, PETSc.DM.BoundaryType.PERIODIC),
stencil_type=PETSc.DMDA.StencilType.STAR,
stencil_width=1,
dof=1)
da.setFromOptions()
da.setUp()
u = da.createGlobalVec()
heat = Head2D(da, D0)
heat.rank = PETSc.COMM_WORLD.getRank()
ts = PETSc.TS().create(comm=PETSc.COMM_WORLD)
ts.setProblemType(PETSc.TS.ProblemType.NONLINEAR)
ts.setDM(da)
ts.setRHSFunction(heat.FormRHSFunctionLocal)
ts.setRHSJacobian(heat.FormRHSJacobianLocal)
if monitor:
ts.setMonitor(heat.EnergyMonitor)
ts.setType(PETSc.TS.Type.BDF)
ts.setTime(0.0)
ts.setMaxTime(0.1)
ts.setTimeStep(0.001)
ts.setExactFinalTime(PETSc.TS.ExactFinalTime.MATCHSTEP)
ts.setFromOptions()
t0 = ts.getTime()
tf = ts.getMaxTime()
mx, my = ts.getDM().getSizes()
PETSc.Sys.Print("solving on %d x %d grid for t0=%g to tf=%g ..."%(mx, my, t0, tf))
u.set(0.0)
ts.solve(u)
u.destroy()
ts.destroy()
da.destroy()