[petsc-dev] Jed will not be happy
Barry Smith
bsmith at mcs.anl.gov
Sat Aug 13 19:02:42 CDT 2011
Jed,
You will need to explain to me all this code or I am going to have to rip it all out since Matt wants a release.
That is, all the DMDA interpolation stuff that was carefully put in under your supervision will be removed unless I understand why it is there. I asked you once before and you didn't really answer.
For example in the 3d interpolation it has a
if (!vcoors) {
then two blocks of code NEITHER of which actually use the coordinate information AT ALL. So why the two blocks of code? Is one wrong and one right and if so why not just get rid of the wrong one. BTW the second version does not work for periodic conditions which is bad! Is your intention to use the coordinate information on each level to generate an interpolation matrix that depends on the coordinates or not? Why is the vcoors flag used as a check here? Are you assuming that the vcoors existing has some special meaning? I understand that the code
Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
comes from the trilinear finite element basic functions to define an interpolation but since you are evaluating them on a uniform mesh
for (li=0; li<nxi; li++) {
xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
}
for (lj=0; lj<neta; lj++) {
eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
}
for (lk=0; lk<nzeta; lk++) {
zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
}
why bother? Isn't this always going to give you a uniform grid interpolation? Also src/dm/examples/tests/ex36.c seems to be suppose to test the interpolation for nonuniform grids but it doesn't seem to work
arry-smiths-macbook-pro:tests barrysmith$ ./ex36 -dim 2 -nl 1 -cmap 1
[2 x 2]=>[4 x 4], interp err = 1.4873e+00
and doesn't have any "nightly" output files so isn't a test at all. Is there a test that is actually run in the nightly for this?
I have been totally confused about this stuff for a long time and it needs to be resolved.
Barry
ierr = DMDAGetCoordinates(daf,&vcoors);CHKERRQ(ierr);
if (vcoors) {
ierr = DMDAGetGhostedCoordinates(dac,&cvcoors);CHKERRQ(ierr);
ierr = DMDAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr);
ierr = DMDAVecGetArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr);
}
/* loop over local fine grid nodes setting interpolation for those*/
if (!vcoors) {
for (l=l_start; l<l_start+p_f; l++) {
for (j=j_start; j<j_start+n_f; j++) {
for (i=i_start; i<i_start+m_f; i++) {
/* convert to local "natural" numbering and then to PETSc global numbering */
row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
i_c = (i/ratioi);
j_c = (j/ratioj);
l_c = (l/ratiok);
/*
Only include those interpolation points that are truly
nonzero. Note this is very important for final grid lines
in x and y directions; since they have no right/top neighbors
*/
x = ((double)(i - i_c*ratioi))/((double)ratioi);
y = ((double)(j - j_c*ratioj))/((double)ratioj);
z = ((double)(l - l_c*ratiok))/((double)ratiok);
nc = 0;
/* one left and below; or we are right on it */
col = dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
cols[nc] = idx_c[col]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
if (i_c*ratioi != i) {
cols[nc] = idx_c[col+dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. - (2.0*z-1.));
}
if (j_c*ratioj != j) {
cols[nc] = idx_c[col+m_ghost_c*dof]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
}
if (l_c*ratiok != l) {
cols[nc] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
}
if (j_c*ratioj != j && i_c*ratioi != i) {
cols[nc] = idx_c[col+(m_ghost_c+1)*dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. - (2.0*z-1.));
}
if (j_c*ratioj != j && l_c*ratiok != l) {
cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;
v[nc++] = .125*(1. - (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
}
if (i_c*ratioi != i && l_c*ratiok != l) {
cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. - (2.0*y-1.))*(1. + (2.0*z-1.));
}
if (i_c*ratioi != i && l_c*ratiok != l && j_c*ratioj != j) {
cols[nc] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof;
v[nc++] = .125*(1. + (2.0*x-1.))*(1. + (2.0*y-1.))*(1. + (2.0*z-1.));
}
ierr = MatSetValues(mat,1,&row,nc,cols,v,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
} else {
PetscScalar *xi,*eta,*zeta;
PetscInt li,nxi,lj,neta,lk,nzeta,n;
PetscScalar Ni[8];
/* compute local coordinate arrays */
nxi = ratioi + 1;
neta = ratioj + 1;
nzeta = ratiok + 1;
ierr = PetscMalloc(sizeof(PetscScalar)*nxi,&xi);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscScalar)*neta,&eta);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscScalar)*nzeta,&zeta);CHKERRQ(ierr);
for (li=0; li<nxi; li++) {
xi[li] = -1.0 + (PetscScalar)li*(2.0/(PetscScalar)(nxi-1));
}
for (lj=0; lj<neta; lj++) {
eta[lj] = -1.0 + (PetscScalar)lj*(2.0/(PetscScalar)(neta-1));
}
for (lk=0; lk<nzeta; lk++) {
zeta[lk] = -1.0 + (PetscScalar)lk*(2.0/(PetscScalar)(nzeta-1));
}
for (l=l_start; l<l_start+p_f; l++) {
for (j=j_start; j<j_start+n_f; j++) {
for (i=i_start; i<i_start+m_f; i++) {
/* convert to local "natural" numbering and then to PETSc global numbering */
row = idx_f[dof*(m_ghost*n_ghost*(l-l_start_ghost) + m_ghost*(j-j_start_ghost) + (i-i_start_ghost))]/dof;
i_c = (i/ratioi);
j_c = (j/ratioj);
l_c = (l/ratiok);
/* remainders */
li = i - ratioi * (i/ratioi);
if (i==mx-1){ li = nxi-1; }
lj = j - ratioj * (j/ratioj);
if (j==my-1){ lj = neta-1; }
lk = l - ratiok * (l/ratiok);
if (l==mz-1){ lk = nzeta-1; }
/* corners */
col = dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c)+m_ghost_c*(j_c-j_start_ghost_c)+(i_c-i_start_ghost_c));
cols[0] = idx_c[col]/dof;
Ni[0] = 1.0;
if ( (li==0) || (li==nxi-1) ) {
if ( (lj==0) || (lj==neta-1) ) {
if ( (lk==0) || (lk==nzeta-1) ) {
ierr = MatSetValue(mat,row,cols[0],Ni[0],INSERT_VALUES);CHKERRQ(ierr);
continue;
}
}
}
/* edges + interior */
/* remainders */
if (i==mx-1){ i_c--; }
if (j==my-1){ j_c--; }
if (l==mz-1){ l_c--; }
col = dof*(m_ghost_c*n_ghost_c*(l_c-l_start_ghost_c) + m_ghost_c*(j_c-j_start_ghost_c) + (i_c-i_start_ghost_c));
cols[0] = idx_c[col]/dof; /* one left and below; or we are right on it */
cols[1] = idx_c[col+dof]/dof; /* one right and below */
cols[2] = idx_c[col+m_ghost_c*dof]/dof; /* one left and above */
cols[3] = idx_c[col+(m_ghost_c+1)*dof]/dof; /* one right and above */
cols[4] = idx_c[col+m_ghost_c*n_ghost_c*dof]/dof; /* one left and below and front; or we are right on it */
cols[5] = idx_c[col+(m_ghost_c*n_ghost_c+1)*dof]/dof; /* one right and below, and front */
cols[6] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c)*dof]/dof;/* one left and above and front*/
cols[7] = idx_c[col+(m_ghost_c*n_ghost_c+m_ghost_c+1)*dof]/dof; /* one right and above and front */
Ni[0] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[1] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0-zeta[lk]);
Ni[2] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[3] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0-zeta[lk]);
Ni[4] = 0.125*(1.0-xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[5] = 0.125*(1.0+xi[li])*(1.0-eta[lj])*(1.0+zeta[lk]);
Ni[6] = 0.125*(1.0-xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
Ni[7] = 0.125*(1.0+xi[li])*(1.0+eta[lj])*(1.0+zeta[lk]);
for (n=0; n<8; n++) {
if( PetscAbsScalar(Ni[n])<1.0e-32) { cols[n]=-1; }
}
ierr = MatSetValues(mat,1,&row,8,cols,Ni,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
ierr = PetscFree(xi);CHKERRQ(ierr);
ierr = PetscFree(eta);CHKERRQ(ierr);
ierr = PetscFree(zeta);CHKERRQ(ierr);
}
More information about the petsc-dev
mailing list