[cig-commits] commit: Working version with new BC's
Mercurial
hg at geodynamics.org
Thu May 29 11:26:22 PDT 2008
changeset: 130:88579a7aa65a
user: Marc Spiegelman <mspieg at ldeo.columbia.edu>
date: Wed Mar 19 12:14:51 2008 -0400
files: src/FieldFunctions.c src/SolitaryWave2DSL.h
description:
Working version with new BC's
M src/FieldFunctions.c
Initialize: put in constant porosity buffer near top (might remove)
FormInitialGuess: modify for dirichlet BC
FormFunctionLocal: start by initializing for periodic or Neumann, the fix dirichlet
FormJacobianLocal: Get all BCs right (the first time round)
BicubicInterp: added ishift, jshift to move into domain for non-periodic BC's
InterpFields2DPeriodic: (use ishift=jshift=0)
InterpFields2DzInflow: check for inflow boundary and shift if necessary - requires recompilation of semi-lagrange
M src/SolitaryWave2DSL.h
changed prototype for BiCubicInterp for ishift,jshift
diff -r 22fca481af77 -r 88579a7aa65a src/FieldFunctions.c
--- a/src/FieldFunctions.c Mon Mar 17 10:57:19 2008 -0400
+++ b/src/FieldFunctions.c Wed Mar 19 12:14:51 2008 -0400
@@ -25,6 +25,7 @@ int Initialize(DMMG *dmmg)
DA da;
int ierr;
int i,j,xs,ys,xm,ym;
+ int jbuffer;
int ni,nj;
PetscReal dx, dz,x,z;
PetscReal amp, z0,eps;
@@ -48,6 +49,9 @@ int Initialize(DMMG *dmmg)
ierr = PetscViewerBinaryOpen(comm,param->inputFilename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
ierr = VecLoadIntoVector(viewer,user->Xold);CHKERRQ(ierr);
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
+ if (param->zBoundaryFlag) {
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: For inflow BC's make sure porosity is 1 near top boundary"); CHKERRQ(ierr);
+ }
} else { /* set up 1-D wave from internal code */
@@ -59,6 +63,9 @@ int Initialize(DMMG *dmmg)
z0= param->z0;
eps = (PetscReal) (param->noiseEps)/RAND_MAX;
param->waveSpeed = 2*amp+1; /* n=3 wave */
+
+ /* set up buffer layer of constant porosity 2 compaction length from top Boundary */
+ jbuffer= nj - 1 - (int ) 2./dz;
/* Compute initial guess */
for (j=ys; j<ys+ym; j++) {
@@ -66,7 +73,11 @@ int Initialize(DMMG *dmmg)
phi = (PetscReal) PhiSolitaryWave(fabs(z-z0), amp);
for (i=xs; i<xs+xm; i++) {
x = dx*i;
- u[j][i].phi = phi+eps*random();
+ if (param->zBoundaryFlag > 0 && j >= jbuffer) {
+ u[j][i].phi = 1.;
+ } else {
+ u[j][i].phi = phi+eps*random();
+ }
u[j][i].cmp=0. ;
}
}
@@ -108,21 +119,23 @@ int Initialize(DMMG *dmmg)
#define __FUNCT__ "FormInitialGuess"
PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
{
-/* DMMG dmmg = (DMMG) ptr; */
- /* AppCtx *user = (AppCtx*) dmmg[dmmg[0].nlevels-1].user; */ /* these should really be on fine grid */
+
AppCtx *user = (AppCtx*) dmmg->user; /* these should really be on fine grid */
TStepCtx *tsCtx = user->tsCtx;
+ PetscBag bag = user->bag;
+ Parameter *param;
DA da = (DA)dmmg->dm;
-/* Vec XLocal; */
Characteristic c = tsCtx->c;
PetscReal t,dt,halfDt;
int ierr;
- int i,j,xs,ys,xm,ym;
+ int i,j,xs,ys,xm,ym,mx,my;
Field **u, **ustar;
/* int nSmooths = 1; */
+ /* retrieve parameters from the bag */
+ ierr = PetscBagGetData(bag,(void**)¶m);CHKERRQ(ierr);
/* Time stepping parameters */
t=tsCtx->t;
@@ -145,6 +158,7 @@ PetscErrorCode FormInitialGuess(DMMG dmm
}
/* Get the fine grid info */
+ ierr = DAGetInfo(da,PETSC_NULL,&mx,&my,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
ierr = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(ierr);
ierr = DAVecGetArray(da,X,(void**)&u); CHKERRQ(ierr);
ierr = DAVecGetArray(da,user->Xstar,(void**)&ustar); CHKERRQ(ierr);
@@ -158,12 +172,17 @@ PetscErrorCode FormInitialGuess(DMMG dmm
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
u[j][i].phi = ustar[j][i].phi+ halfDt*(ustar[j][i].cmp+u[j][i].cmp);
-/* u[j][i].phi = ustar[j][i].phi+ halfDt*ustar[j][i].cmp; */
-/* u[j][i].cmp = 0.; */
- }
- }
-
-
+ }
+ }
+
+ /* fix dirichlet boundary */
+ if ( param->zBoundaryFlag == 1 && ys+ym == my ) {
+ j = ys + ym -1;
+ for (i=xs; i<xs+xm; i++) {
+ u[j][i].phi = ustar[j][i].phi;
+ }
+ }
+
ierr = DAVecRestoreArray(da,X,(void**)&u); CHKERRQ(ierr);
ierr = DAVecRestoreArray(da,user->Xstar,(void**)&ustar); CHKERRQ(ierr);
return 0;
@@ -182,7 +201,7 @@ PetscErrorCode FormFunctionLocal(DALocal
Field **xstar;
DA da=info->da;
- int i,j;
+ int i,j, im,ip,jm,jp;
int ierr, mx, my;
int xints,xinte,yints,yinte;
PetscReal nExp;
@@ -226,26 +245,61 @@ PetscErrorCode FormFunctionLocal(DALocal
/*
porosity and compaction rate equation for all points (periodic plus ghosts should require
- no boundary conditions)
+ no boundary conditions; initialize with all neumann than patch up dirichlet)
*/
for (j=yints; j<yinte; j++) {
+ jm = j-1;
+ jp = j+1;
+ if ( param->zBoundaryFlag == 1 ) {
+ if ( jp == my ) jp = j-1 ;
+ if ( jm < 0 ) jm = 1 ;
+ }
for (i=xints; i<xinte; i++) {
+
+ im = i-1;
+ ip = i+1;
+ if (param->xBoundaryFlag == 1 ) { /* use reflection in x */
+ if (im < 0 ) im = 1;
+ if (ip == mx ) ip = mx - 2;
+ }
+
/* porosity update equation */
f[j][i].phi = x[j][i].phi - xstar[j][i].phi - halfDt*(xstar[j][i].cmp+x[j][i].cmp);
/* Calculate permeabilities at half nodes */
- khlfE = pow((.5*(x[j][i].phi+x[j][i+1].phi)),nExp);
- khlfW = pow((.5*(x[j][i].phi+x[j][i-1].phi)),nExp);
- khlfN = pow((.5*(x[j][i].phi+x[j+1][i].phi)),nExp);
- khlfS = pow((.5*(x[j][i].phi+x[j-1][i].phi)),nExp);
+ khlfE = pow((.5*(x[j][i].phi+x[j][ip].phi)),nExp);
+ khlfW = pow((.5*(x[j][i].phi+x[j][im].phi)),nExp);
+ khlfN = pow((.5*(x[j][i].phi+x[jp][i].phi)),nExp);
+ khlfS = pow((.5*(x[j][i].phi+x[jm][i].phi)),nExp);
sum = khlfN + khlfS + dzdx2*(khlfE + khlfW) + dz2;
/* compaction rate equation */
- f[j][i].cmp = -(dzdx2*(khlfE*x[j][i+1].cmp+khlfW*x[j][i-1].cmp) + khlfN*x[j+1][i].cmp + khlfS*x[j-1][i].cmp) + sum*x[j][i].cmp + dz*(khlfN-khlfS);
-
- }
- }
+ f[j][i].cmp = -(dzdx2*(khlfE*x[j][ip].cmp+khlfW*x[j][im].cmp) + khlfN*x[jp][i].cmp + khlfS*x[jm][i].cmp) + sum*x[j][i].cmp + dz*(khlfN-khlfS);
+
+ }
+ }
+
+ /* patch up dirichlet condition if on top boundary */
+
+ if ( param->zBoundaryFlag == 1 ) {
+ if (yinte == my ) {
+ j= my - 1;
+ for (i=xints; i<xinte; i++) {
+ f[j][i].phi = x[j][i].phi - xstar[j][i].phi; /* this should set phi=1 */
+ f[j][i].cmp = x[j][i].cmp;
+ }
+ }
+ if ( yints == 0 ) { /*free-flux condition on the bottom */
+ j=0;
+ for (i=xints; i<xinte; i++) {
+ khlfN = pow((.5*(x[j][i].phi+x[j+1][i].phi)),nExp);
+ khlfS = pow(x[j][i].phi,nExp);
+ f[j][i].cmp += .5*dz*(khlfN-khlfS);
+ }
+ }
+ }
+
DAVecRestoreArray(da,user->Xstar,(void**)&xstar);
PetscFunctionReturn(0);
@@ -271,6 +325,7 @@ PetscErrorCode FormJacobianLocal(DALocal
int i,j,ierr;
+ int im,ip,jm,jp;
int xints,xinte,yints,yinte;
PetscReal nExp;
int mx,my;
@@ -282,6 +337,8 @@ PetscErrorCode FormJacobianLocal(DALocal
PetscFunctionBegin;
+ ierr = MatZeroEntries(jac); CHKERRQ(ierr);
+
/* retrieve parameters from the bag */
ierr = PetscBagGetData(bag,(void**)¶m);CHKERRQ(ierr);
nExp = param->nExp;
@@ -317,16 +374,31 @@ PetscErrorCode FormJacobianLocal(DALocal
/*
porosity and compaction rate stencils for all points (periodic plus ghosts should require
- no boundary conditions)
+ no boundary conditions; otherwise initialize with all neumann than patch up dirichlet))
*/
for (j=yints; j<yinte; j++) {
+ jm = j-1;
+ jp = j+1;
+ if ( param->zBoundaryFlag == 1 ) {
+ if ( jp == my ) jp = j-1 ;
+ if ( jm < 0 ) jm = 1 ;
+ }
for (i=xints; i<xinte; i++) {
- /* Calculate porosities and permeabilities at half nodes */
- phlfE = .5*(x[j][i].phi+x[j][i+1].phi);
- phlfW = .5*(x[j][i].phi+x[j][i-1].phi);
- phlfN = .5*(x[j][i].phi+x[j+1][i].phi);
- phlfS = .5*(x[j][i].phi+x[j-1][i].phi);
+
+ im = i-1;
+ ip = i+1;
+ if (param->xBoundaryFlag == 1 ) { /* use reflection in x */
+ if (im < 0 ) im = 1;
+ if (ip == mx ) ip = mx - 2;
+ }
+
+ /* Calculate porosities and permeabilities at half nodes */
+
+ phlfE = .5*(x[j][i].phi+x[j][ip].phi);
+ phlfW = .5*(x[j][i].phi+x[j][im].phi);
+ phlfN = .5*(x[j][i].phi+x[jp][i].phi);
+ phlfS = .5*(x[j][i].phi+x[jm][i].phi);
khlfE = pow(phlfE,nExp);
khlfW = pow(phlfW,nExp);
@@ -340,41 +412,49 @@ PetscErrorCode FormJacobianLocal(DALocal
khlfN*x[j+1][i].cmp + khlfS*x[j-1][i].cmp + sum*x[j][i].cmp - dz*(khlfN-khlfS);
*/
- nColumns = 5;
-
- row.i = i; row.j = j; row.c = 0; /* components of Jacobian for dF/dCmp */
- v[0] = -khlfS; col[0].i = i; col[0].j = j-1; col[0].c = 0;
- v[1] = -dzdx2*khlfW; col[1].i = i-1; col[1].j = j; col[1].c = 0;
- v[2] = sum; col[2].i = i; col[2].j = j; col[2].c = 0;
- v[3] = -dzdx2*khlfE; col[3].i = i+1; col[3].j = j; col[3].c = 0;
- v[4] = -khlfN; col[4].i = i; col[4].j = j+1; col[4].c = 0;
-
-
- if ( tsCtx->takeFullStep) {
- nColumns = 10;
- cmpC=x[j][i].cmp; /* components of Jacobian for dF/dPhi */
- v[5] = -halfNExp*khlfS/phlfS*(x[j-1][i].cmp-cmpC+dz); col[5].i = i; col[5].j = j-1; col[5].c = 1;
- v[6] = -halfNExp*khlfN/phlfN*(x[j+1][i].cmp-cmpC-dz); col[6].i = i; col[6].j = j+1; col[6].c = 1;
- v[7] = -dzdx2*halfNExp*khlfW/phlfW*(x[j][i-1].cmp-cmpC); col[7].i = i-1; col[7].j = j; col[7].c = 1;
- v[8] = -dzdx2*halfNExp*khlfE/phlfE*(x[j][i+1].cmp-cmpC); col[8].i = i+1; col[8].j = j; col[8].c = 1;
- v[9] = v[5]+v[6]+v[7]+v[8]; col[9].i = i; col[9].j = j; col[9].c = 1;
- }
-
- MatSetValuesStencil(jac,1,&row,nColumns,col,v,INSERT_VALUES);
+ row.i = i; row.j = j; row.c = 0; /* components of 11_Block of Jacobian for dF/dCmp */
+ if (param->zBoundaryFlag == 1 && j == my-1 ) { /* if dirichlet boundary */
+ nColumns = 1;
+ v[0]=1; col[0].i = i; col[0].j = j; col[0].c = 0;
+ } else {
+ nColumns = 5;
+
+ v[0] = -khlfS; col[0].i = i; col[0].j = jm; col[0].c = 0;
+ v[1] = -dzdx2*khlfW; col[1].i = im; col[1].j = j; col[1].c = 0;
+ v[2] = sum; col[2].i = i; col[2].j = j; col[2].c = 0;
+ v[3] = -dzdx2*khlfE; col[3].i = ip; col[3].j = j; col[3].c = 0;
+ v[4] = -khlfN; col[4].i = i; col[4].j = jp; col[4].c = 0;
+
+
+ if ( tsCtx->takeFullStep) {
+ nColumns = 10;
+ cmpC=x[j][i].cmp; /* components of 12_Block of Jacobian for dF/dPhi */
+ v[5] = -halfNExp*khlfS/phlfS*(x[jm][i].cmp-cmpC+dz); col[5].i = i; col[5].j = jm; col[5].c = 1;
+ v[6] = -halfNExp*khlfN/phlfN*(x[jp][i].cmp-cmpC-dz); col[6].i = i; col[6].j = jp; col[6].c = 1;
+ v[7] = -dzdx2*halfNExp*khlfW/phlfW*(x[j][im].cmp-cmpC); col[7].i = im; col[7].j = j; col[7].c = 1;
+ v[8] = -dzdx2*halfNExp*khlfE/phlfE*(x[j][ip].cmp-cmpC); col[8].i = ip; col[8].j = j; col[8].c = 1;
+ v[9] = v[5]+v[6]+v[7]+v[8]; col[9].i = i; col[9].j = j; col[9].c = 1;
+ }
+ }
+
+ MatSetValuesStencil(jac,1,&row,nColumns,col,v,ADD_VALUES);
/*
Jacobian entries for porosity update equation
f[j][i].phi = x[j][i].phi - xold[j][i].phi - halfDt*(xold[j][i].cmp+x[j][i].cmp);
*/
- nColumns = 2;
-
row.i = i; row.j = j; row.c = 1;
- v[0] = one; col[0].i = i; col[0].j = j; col[0].c = 1;
- v[1] = -halfDt; col[1].i = i; col[1].j = j; col[1].c = 0;
-
- MatSetValuesStencil(jac,1,&row,nColumns,col,v,INSERT_VALUES);
-
+ if (param->zBoundaryFlag == 1 && j == my-1 ) { /* if dirichlet boundary */
+ nColumns = 1;
+ v[0]=1; col[0].i = i; col[0].j = j; col[0].c = 1;
+ } else {
+ nColumns = 2;
+ v[0] = one; col[0].i = i; col[0].j = j; col[0].c = 1;
+ v[1] = -halfDt; col[1].i = i; col[1].j = j; col[1].c = 0;
+ }
+
+ MatSetValuesStencil(jac,1,&row,nColumns,col,v,ADD_VALUES);
}
}
@@ -413,6 +493,7 @@ PetscErrorCode InterpFields2DPeriodic(vo
AppCtx *user = (AppCtx *) ctx;
Field **x = (Field **) f;
int ni=user->grid->niFine, nj=user->grid->njFine;
+ PetscInt ishift=0, jshift=0;
PetscReal ir=ij_real[0], jr=ij_real[1];
int ierr;
@@ -422,7 +503,7 @@ PetscErrorCode InterpFields2DPeriodic(vo
if ( ir < 0 || ir > ni-1 || jr < 0 || jr> nj-1 ) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: InterpFields2D: ir=%g, jr=%g\n",ir,jr); CHKERRQ(ierr);
}
- BiCubicInterp((TwoVector **) x, ir, jr, (TwoVector *) field);
+ BiCubicInterp((TwoVector **) x, ir, jr, ishift, jshift, (TwoVector *) field);
return 0;
}
@@ -439,6 +520,7 @@ PetscErrorCode InterpFields2DzInflow(voi
AppCtx *user = (AppCtx *) ctx;
Field **x = (Field **) f;
int ni=user->grid->niFine, nj=user->grid->njFine;
+ PetscInt ishift=0, jshift=0;
PetscReal ir=ij_real[0], jr=ij_real[1];
int ierr;
@@ -454,20 +536,26 @@ PetscErrorCode InterpFields2DzInflow(voi
field[0] = 0.;
field[1] = 1.;
} else {
- BiCubicInterp((TwoVector **) x, ir, jr, (TwoVector *) field);
+ /* make sure to interpolate within bounds */
+ ishift =0; jshift = 0;
+ if ( ir < 1 ) ishift = 1;
+ if ( ir >= ni-2 ) ishift = -1;
+ if ( jr < 1) jshift = 1;
+ if ( jr >= nj-2) jshift = -1;
+ BiCubicInterp((TwoVector **) x, ir, jr, ishift, jshift, (TwoVector *) field);
}
return 0;
}
/*---------------------------------------------------------------------*/
#undef __FUNCT__
#define __FUNCT__ "BiCubicInterp"
-PetscErrorCode BiCubicInterp(TwoVector **x, PetscReal ir, PetscReal jr, TwoVector *xinterp)
+PetscErrorCode BiCubicInterp(TwoVector **x, PetscReal ir, PetscReal jr, PetscInt ishift, PetscInt jshift, TwoVector *xinterp)
/*---------------------------------------------------------------------*/
{
int im, jm, imm,jmm,ip,jp,ipp,jpp;
int k;
PetscReal il, jl, row1, row2, row3, row4;
- im = (int)floor(ir); jm = (int)floor(jr);
+ im = (int)floor(ir)+ishift; jm = (int)floor(jr)+jshift;
il = ir - im + 1.0; jl = jr - jm + 1.0;
imm = im-1; ip = im+1; ipp = im+2;
jmm = jm-1; jp = jm+1; jpp = jm+2;
diff -r 22fca481af77 -r 88579a7aa65a src/SolitaryWave2DSL.h
--- a/src/SolitaryWave2DSL.h Mon Mar 17 10:57:19 2008 -0400
+++ b/src/SolitaryWave2DSL.h Wed Mar 19 12:14:51 2008 -0400
@@ -120,7 +120,7 @@ PetscErrorCode InterpVelocity2D(void *,
PetscErrorCode InterpVelocity2D(void *, PetscReal [], PetscInt, PetscInt [], PetscScalar [], void *);
PetscErrorCode InterpFields2DPeriodic(void *, PetscReal [], PetscInt, PetscInt [], PetscReal [], void *);
PetscErrorCode InterpFields2DzInflow(void *, PetscReal [], PetscInt, PetscInt [], PetscReal [], void *);
-PetscErrorCode BiCubicInterp(TwoVector **, PetscScalar, PetscScalar, TwoVector*);
+PetscErrorCode BiCubicInterp(TwoVector **, PetscScalar, PetscScalar, PetscInt, PetscInt, TwoVector*);
PetscReal CubicInterp(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar);
More information about the cig-commits
mailing list