Skip to content

Commit

Permalink
Cleanup of unused variables, restrict chebyshev.[ch] to C89.
Browse files Browse the repository at this point in the history
  • Loading branch information
@jedbrown
jedbrown committed Jul 17, 2008
1 parent c4838ef commit 7b5a751
Show file tree
Hide file tree
Showing 4 changed files with 27 additions and 34 deletions.
25 changes: 14 additions & 11 deletions chebyshev.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ PetscErrorCode MatCreateChebD1(MPI_Comm comm, Vec vx, Vec vy, unsigned flag, Mat
#define __FUNCT__ "ChebD1Mult"
PetscErrorCode ChebD1Mult(Mat A, Vec vx, Vec vy) {
PetscErrorCode ierr;
PetscInt i;
double *x, *y;
ChebD1Ctx *c;

Expand All @@ -47,7 +48,7 @@ PetscErrorCode ChebD1Mult(Mat A, Vec vx, Vec vy) {
int n = c->n - 1; // Gauss-Lobatto-Chebyshev points are numbered from [0..n]

fftw_execute_r2r(c->p_forward, x, c->work);
for (int i = 1; i < n; i++) c->work[i] *= (double)i;
for (i = 1; i < n; i++) c->work[i] *= (double)i;

fftw_execute_r2r(c->p_backward, c->work + 1, y + 1);

Expand All @@ -56,7 +57,7 @@ PetscErrorCode ChebD1Mult(Mat A, Vec vx, Vec vy) {
y[0] = 0.0;
y[n] = 0.0;
double s = 1.0;
for (int i = 1; i < n; i++) {
for (i = 1; i < n; i++) {
double I = (double)i;
y[i] /= 2.0 * n * sqrt(1.0 - PetscSqr(cos(I * pin)));
y[0] += I * c->work[i];
Expand Down Expand Up @@ -89,7 +90,7 @@ PetscErrorCode MatCreateCheb(MPI_Comm comm, int rank, int tr, int *dim, unsigned
Vec vx, Vec vy, Mat *A) {
PetscErrorCode ierr;
ChebCtx *c;
int n;
int n, r, ri, stride;
double *x, *y;

PetscFunctionBegin;
Expand All @@ -103,8 +104,8 @@ PetscErrorCode MatCreateCheb(MPI_Comm comm, int rank, int tr, int *dim, unsigned
c->tr = tr;

if (!(0 <= tr && tr < rank)) SETERRQ(PETSC_ERR_USER, "tdim out of range");
int stride = 1;
for (int r = rank-1, ri=rank-2; r >= 0; r--) {
stride = 1;
for (r = rank-1, ri=rank-2; r >= 0; r--) {
fftw_iodim *iod;
if (r == tr) {
iod = &(c->tdim);
Expand Down Expand Up @@ -141,6 +142,8 @@ PetscErrorCode MatCreateCheb(MPI_Comm comm, int rank, int tr, int *dim, unsigned
PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy) {
PetscErrorCode ierr;
double *x, *y;
bool done;
int i;
ChebCtx *c;

PetscFunctionBegin;
Expand All @@ -156,13 +159,13 @@ PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy) {
double N = (double)n;
ierr = PetscMemzero(ind, (c->rank - 1) * sizeof(int)); CHKERRQ(ierr);
int offset = 0;
for (bool done = false; !done; ) {
for (done = false; !done; ) {
int ix0 = offset;
int ixn = offset + n * c->tdim.is;
y[ix0] = 0.0;
y[ixn] = 0.0;
double s = 1.0;
for (int i = 1; i < n; i++) {
for (i = 1; i < n; i++) {
int ix = offset + i * c->tdim.is;
double I = (double)i;
c->work[ix] *= I;
Expand All @@ -180,8 +183,8 @@ PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy) {
double pin = PI / N;
ierr = PetscMemzero(ind, (c->rank - 1) * sizeof(int)); CHKERRQ(ierr);
offset = 0;
for (bool done = false; !done; ) {
for (int i = 1; i < n; i++) {
for (done = false; !done; ) {
for (i = 1; i < n; i++) {
int ix = offset + i * c->tdim.is;
double I = (double)i;
y[ix] /= 2 * n * sqrt(1.0 - PetscSqr(cos(I * pin)));
Expand All @@ -197,11 +200,11 @@ PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy) {


void perform_carry(ChebCtx *c, int *ind, int *offset, bool *done) {
int carry;
int carry, i;

*offset = 0;
carry = 1;
for (int i = c->rank-2; i >= 0; i--) {
for (i = c->rank-2; i >= 0; i--) {
ind[i] += carry;
if (ind[i] < c->dim[i].n) {
carry = 0;
Expand Down
4 changes: 4 additions & 0 deletions chebyshev.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
#include <fftw3.h>
#include <petscmat.h>

PETSC_EXTERN_CXX_BEGIN

#define PI 3.14159265358979323846

typedef struct {
Expand All @@ -31,4 +33,6 @@ PetscErrorCode MatCreateCheb(MPI_Comm comm, int rank, int tr, int *dims, unsigne
PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy);
PetscErrorCode ChebDestroy (Mat A);

PETSC_EXTERN_CXX_END

#endif // CHEBYSHEV_H
28 changes: 8 additions & 20 deletions stokes.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -546,7 +546,6 @@ PetscErrorCode StokesMatGetDiagonalSchur(Mat S, Vec y)
PetscErrorCode StokesMatMultPV(Mat A, Vec xG, Vec yG) // divergence of velocity
{
StokesCtx *ctx;
PetscReal **u, *v;
PetscErrorCode ierr;

PetscFunctionBegin;
Expand All @@ -559,7 +558,6 @@ PetscErrorCode StokesMatMultPV(Mat A, Vec xG, Vec yG) // divergence of velocity
#define __FUNCT__ "StokesDivergence"
PetscErrorCode StokesDivergence(StokesCtx *ctx, PetscTruth withDirichlet, Vec xG, Vec yG) // divergence of velocity
{
PetscReal **u, *v;
PetscErrorCode ierr;

ierr = VecZeroEntries(ctx->workV[0]);CHKERRQ(ierr);
Expand Down Expand Up @@ -752,7 +750,6 @@ PetscErrorCode StokesFunction(SNES snes, Vec xG, Vec yG, void *ctx)
#define __FUNCT__ "StokesJacobian"
PetscErrorCode StokesJacobian(SNES snes, Vec w, Mat *Ashell, Mat *Pshell, MatStructure *flag, void *void_ctx)
{
PetscErrorCode ierr;

PetscFunctionBegin;
// The nonlinear term has already been fixed up by StokesFunction() so we do nothing here.
Expand All @@ -767,7 +764,6 @@ PetscErrorCode StokesSetupDomain(StokesCtx *c, Vec *global)
MPI_Comm comm = c->comm;
PetscInt d = c->numDims, *dim = c->dim, m = productInt(d, dim), n = d*m, N=n+m;
StokesBoundaryMixed *mixed;
IS isG, isD;
PetscInt *ixLP, *ixPL, *ixLV, *ixVL, *ixLD, *ixDL, *ixPG, *ixGP, *ixVG, *ixGV;
PetscInt lp, lv, gp, gv, dv, g, im;
PetscReal *v, *w, *x;
Expand Down Expand Up @@ -1156,7 +1152,6 @@ PetscErrorCode StokesPCSetUp0(void *void_ctx)
PetscInt d = ctx->numDims, *dim = ctx->dim;
PetscReal *x, *eta, *deta, **strain;
PetscInt *ixL;
PetscTruth flg;
PetscErrorCode ierr;

PetscFunctionBegin;
Expand All @@ -1168,7 +1163,7 @@ PetscErrorCode StokesPCSetUp0(void *void_ctx)
{
PetscInt row, col[2*d+1], c, im=0;
PetscScalar v[2*d+1];
PetscScalar x0, xMM, xPP, xM, idxM, xP, idxP, idx, eM, deM, du0M, eP, deP, du0P;
PetscScalar x0, xMM, xPP, xM, idxM, xP, idxP, idx, eM, deM, eP, deP;
for (BlockIt it = BlockIt(d, dim); !it.done; it.next()) { // loop over local dof
const PetscInt i = it.i;
if (im < ctx->numMixed && i == ctx->mixed[im].localIndex) { // we are at a mixed boundary node
Expand Down Expand Up @@ -1237,13 +1232,13 @@ PetscErrorCode StokesPCSetUp1(void *void_ctx)
{
#define ORDER 3
#if (ORDER == 2)
const PetscReal qpoint[2] = { -0.57735026918962573, 0.57735026918962573 };
//const PetscReal qpoint[2] = { -0.57735026918962573, 0.57735026918962573 };
const PetscReal qweight[2] = { 1.0, 1.0 };
const PetscReal basis[2][2] = {{0.78867513459481287, 0.21132486540518708}, {0.21132486540518708, 0.78867513459481287}};
const PetscReal deriv[2][2] = {{-0.5, -0.5}, {0.5, 0.5}};
const PetscInt qdim[] = {2,2,2,2,2}; // 2 quadrature points in each direction
#elif (ORDER == 3)
const PetscReal qpoint[3] = {-0.7745966692414834, -2.4651903288156619e-31, 0.7745966692414834};
//const PetscReal qpoint[3] = {-0.7745966692414834, -2.4651903288156619e-31, 0.7745966692414834};
const PetscReal qweight[3] = {0.55555555555556, 0.88888888888889, 0.55555555555556};
const PetscReal basis[2][3] = {{0.887298334621,0.5,0.112701665379},{0.112701665379,0.5,0.887298334621}};
const PetscReal deriv[2][3] = {{-0.5, -0.5, -0.5},{0.5, 0.5, 0.5}};
Expand All @@ -1259,7 +1254,6 @@ PetscErrorCode StokesPCSetUp1(void *void_ctx)
PetscInt *ixL;
PetscInt row[N*d], col[N*d];
PetscReal A[N*d][N*d], M[N*d][N*d];
PetscTruth flg;
PetscErrorCode ierr;

PetscFunctionBegin;
Expand Down Expand Up @@ -1437,7 +1431,7 @@ PetscErrorCode StokesPCSetUp1(void *void_ctx)
PetscErrorCode StokesPCSetUp2(void *void_ctx)
{
StokesCtx *ctx = (StokesCtx *)void_ctx;
PetscInt d = ctx->numDims, *dim = ctx->dim, N = productInt(d,dim), m = d*(4*d+1);
PetscInt d = ctx->numDims, *dim = ctx->dim, m = d*(4*d+1);
PetscInt n, row, col[m];
PetscInt *ixL;
PetscScalar values[m];
Expand Down Expand Up @@ -1508,7 +1502,6 @@ PetscErrorCode StokesPCSetUp3(void *void_ctx)
PetscInt d = ctx->numDims, *dim = ctx->dim;
PetscReal *x, *Eta, *dEta, **Strain;
PetscInt *ixL;
PetscTruth flg;
PetscErrorCode ierr;

PetscFunctionBegin;
Expand All @@ -1526,10 +1519,10 @@ PetscErrorCode StokesPCSetUp3(void *void_ctx)
if (true) { // We are at an interior or Dirichlet node
if (ixL[node.i*d+0] < 0) continue;
if (ixL[node.i*d+1] < 0) continue;
const PetscReal J[2][2] = { x[node.shift(0,1)*d+0] - x[node.shift(0,-1)*d+0],
x[node.shift(0,1)*d+1] - x[node.shift(0,-1)*d+1],
x[node.shift(1,1)*d+0] - x[node.shift(1,-1)*d+0],
x[node.shift(1,1)*d+1] - x[node.shift(1,-1)*d+1] };
const PetscReal J[2][2] = { { x[node.shift(0,1)*d+0] - x[node.shift(0,-1)*d+0],
x[node.shift(0,1)*d+1] - x[node.shift(0,-1)*d+1] },
{ x[node.shift(1,1)*d+0] - x[node.shift(1,-1)*d+0],
x[node.shift(1,1)*d+1] - x[node.shift(1,-1)*d+1] } };
const PetscReal iJdet = 1.0 / (J[0][0]*J[1][1] - J[0][1]*J[1][0]);
const PetscReal Jinv[2*2] = { iJdet*J[1][1], -iJdet*J[0][1], -iJdet*J[1][0], iJdet*J[0][0] };
ierr = StokesComputeNodalJacobian(node, Jinv, x, Eta, dEta, Strain, nodeJac);CHKERRQ(ierr);
Expand Down Expand Up @@ -1581,7 +1574,6 @@ PetscErrorCode StokesComputeNodalJacobian(const BlockIt node, const PetscReal *J
using namespace CppAD;
const PetscInt d = node.numDims(), S = 2*d+1;
vector< AD<double> > vel(S*d), residual(d), flux(S*d*d);
PetscErrorCode ierr;

PetscFunctionBegin;
Independent(vel);
Expand Down Expand Up @@ -1963,7 +1955,6 @@ PetscErrorCode StokesExact2(PetscInt d, PetscReal *coord, PetscReal *value, Pets
#define __FUNCT__ "StokesExact3"
PetscErrorCode StokesExact3(PetscInt d, PetscReal *coord, PetscReal *value, PetscReal *rhs, void *ctx)
{
const PetscReal eta = 1.0;
PetscReal u, v, p;

PetscFunctionBegin;
Expand Down Expand Up @@ -2087,7 +2078,6 @@ PetscErrorCode StokesBoundary2(PetscInt d, PetscReal *coord, PetscReal *normal,
PetscErrorCode StokesBoundary3(PetscInt d, PetscReal *coord, PetscReal *normal, BdyType *type, PetscReal *value, void *void_ctx)
{
bool inside;
PetscErrorCode ierr;

PetscFunctionBegin;
inside = false;
Expand Down Expand Up @@ -2116,8 +2106,6 @@ PetscErrorCode StokesBoundary3(PetscInt d, PetscReal *coord, PetscReal *normal,
#define __FUNCT__ "StokesBoundary4"
PetscErrorCode StokesBoundary4(PetscInt d, PetscReal *coord, PetscReal *normal, BdyType *type, PetscReal *value, void *void_ctx)
{
bool inside;
PetscErrorCode ierr;

PetscFunctionBegin;
*type = DIRICHLET;
Expand Down
4 changes: 1 addition & 3 deletions util.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ class BlockIt {
i = 0;
done = s == 0;
}
BlockIt(const BlockIt &it) : d(it.d), i(it.i), done(it.done) {
BlockIt(const BlockIt &it) : done(it.done), i(it.i), d(it.d) {
dim = new int[d]; stride = new int[d]; ind = new int[d];
for (int i=0; i < d; i++) {
dim[i] = it.dim[i];
Expand Down Expand Up @@ -128,9 +128,7 @@ void zeroInt(int d, int v[]) {
#define __FUNCT__ "polyInterp"
PetscErrorCode polyInterp(const PetscInt n, const PetscReal *x, PetscScalar *w, const PetscReal x0, const PetscReal x1, PetscScalar *f0, PetscScalar *f1)
{
PetscScalar *tmp;
PetscInt o, e;
PetscReal y;

PetscFunctionBegin;
for (int di=1; di < n; di++) { // offset (column of table)
Expand Down

0 comments on commit 7b5a751

Please sign in to comment.