[petsc-dev] broken DMGetInterpolation_DA !!!!

[Barry Smith]

On Jul 21, 2011, at 5:00 PM, Jed Brown wrote:

> On Thu, Jul 21, 2011 at 14:48, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> 
>  Whoever wrote the DMGetInterpolation_DA_Q1() for when coordinates are provided in DA
> 
>   1) it is completely broken for the periodic case!
> 
> There is a serious semantic problem with periodic non-regular coordinates: we don't know how big the wrap element is. There is no way to determine this automatically, and there is no way to store it.

   So maybe generate an error message??? Producing the incorrect answer is not a good approach.

>  
> 
>   2) it doesn't actually use the coordinates provided by the user in the DA ; it uses uniform coordinates it cooks up on the fly
> 
> /* compute local coordinate arrays */
>    nxi  = ratioi + 1;
>    neta = ratioj + 1;
>    ierr = PetscMalloc(sizeof(PetscScalar)*nxi,&xi);CHKERRQ(ierr);
>    ierr = PetscMalloc(sizeof(PetscScalar)*neta,&eta);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));
>    }
> 
> These are reference coordinates. They get mapped using a Q1 local coordinate mapping. Dave wrote this code, but it was tested to give second order accuracy for non-uniform mapped grids.

   I completely do not understand the situation.  In the code

  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);
  }

but the variable coors and ccoors are NEVER used in the routine!  How can this routine supposedly support non-uniform mapped grids when the coordinates are never even taken into account? And why is 

 ierr = DMDAVecGetArray(ddf->da_coordinates,vcoors,&coors);CHKERRQ(ierr);
    ierr = DMDAVecGetArray(ddc->da_coordinates,cvcoors,&ccoors);CHKERRQ(ierr);

even called if the result is not used?   

What am I missing?




> 
> Are you using non-uniform periodic grids?

No, they are uniform. But I want to store the coordinates explicitly.


> If you have regular periodic grids, but you still want to store coordinates explicitly, then we could write some special case to detect that you set coordinates that happened to be uniform, therefore we can guess that the wrap cell should also be uniform, and the code for no provided coordinates would be correct. What do you want to do?