| 1 | /* 2-dimensional Laplace equation solver with constant boundaries.
|
|---|
| 2 | * This operation 2nd order accurate. To verify with CIVL:
|
|---|
| 3 | * civl verify laplace.cvl
|
|---|
| 4 | */
|
|---|
| 5 |
|
|---|
| 6 | $input int rows; // the number of rows
|
|---|
| 7 | $assume(rows >= 1);
|
|---|
| 8 | $input int cols; // the number of columns
|
|---|
| 9 | $assume(cols >= 1);
|
|---|
| 10 | $input double h; // distance between two consecutive points in horizontal or vertical direction
|
|---|
| 11 | $assume(0<h && h<1);
|
|---|
| 12 | $input double u[rows][cols]; // initial state
|
|---|
| 13 | double result[rows][cols]; // new state after update
|
|---|
| 14 | // assume phi:R^2->R has 4 continuous derivatives in [-1,1]x[-1,1]:
|
|---|
| 15 | $abstract $differentiable(4, [-1,1][-1,1]) $real phi($real x, $real y);
|
|---|
| 16 |
|
|---|
| 17 | /* Performs one update step */
|
|---|
| 18 | void laplace() {
|
|---|
| 19 | /*@ loop invariant 1<=i && i<=rows-1;
|
|---|
| 20 | @ loop invariant $forall (int i0 : 1..i-1)
|
|---|
| 21 | @ $forall (int j:1..cols-1)
|
|---|
| 22 | @ result[i0][j] = (u[i0-1][j]+u[i0][j-1]-4*u[i0][j]+u[i0+1][j]+u[i0][j+1])/(h*h);
|
|---|
| 23 | @ loop assigns i, result[1..rows-2][1..cols-2];
|
|---|
| 24 | @*/
|
|---|
| 25 | for (int i = 1; i < rows-1; i++) {
|
|---|
| 26 | /*@ loop invariant 1<=j && j<=cols-1;
|
|---|
| 27 | @ loop invariant $forall (int j0 : 1..j-1)
|
|---|
| 28 | @ result[i][j0] = (u[i-1][j0]+u[i][j0-1]-4*u[i][j0]+u[i+1][j0]+u[i][j0+1])/(h*h);
|
|---|
| 29 | @ loop assigns j, result[i][1..cols-2];
|
|---|
| 30 | @*/
|
|---|
| 31 | for (int j = 1; j < cols-1; j++) {
|
|---|
| 32 | result[i][j] = (u[i-1][j]+u[i][j-1]-4*u[i][j]+u[i+1][j]+u[i][j+1])/(h*h);
|
|---|
| 33 | }
|
|---|
| 34 | }
|
|---|
| 35 | }
|
|---|
| 36 |
|
|---|
| 37 | void main() {
|
|---|
| 38 | $assume($forall (int i,j:($domain(2)){0..rows-1, 0..cols-1}) u[i][j] == phi(i*h, j*h));
|
|---|
| 39 | $assume(rows*h<=1 && cols*h<=1);
|
|---|
| 40 | laplace();
|
|---|
| 41 | $assert($uniform (int i,j:($domain(2)){1..rows-2, 1..cols-2})
|
|---|
| 42 | result[i][j] - $D[phi,{x,2}](i*h,j*h) - $D[phi,{y,2}](i*h,j*h) == $O(h*h));
|
|---|
| 43 | }
|
|---|