Source

PetIGA-FMM / fmm / fmm.c

Full commit
#include <assert.h>
#include "petiga.h"
#include "fmm.h"
#include "nbody.h"
#if PETSC_VERSION_LE(3,2,0)
#include <private/pcimpl.h>
#else
#include <petsc-private/pcimpl.h>
#endif

#define DENSE 0

#undef  __FUNCT__
#define __FUNCT__ "PCFMMSetIGA"
PetscErrorCode PCFMMSetIGA(PC pc, IGA iga)
{
  PC_FMM *fmm = (PC_FMM*)pc->data;
  PetscFunctionBegin;
  PetscValidHeaderSpecific(iga,IGA_CLASSID,2);
  if (iga) fmm->iga = iga;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "BEM_MatMult"
static PetscErrorCode BEM_MatMult(Mat G, Vec x, Vec y)
{
  PC_FMM *fmm;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = MatShellGetContext(G,(void**)&fmm);CHKERRQ(ierr);
  PetscInt i, nb = fmm->nb;
  PetscReal *xl, *yl;
  PetscReal *xm = fmm->xm;
  PetscReal *ym = fmm->ym;
  PetscReal *re = fmm->re;
  VecGetSize(x,&i);
  ierr = VecGetArray(x,&xl);CHKERRQ(ierr);
  ierr = VecGetArray(y,&yl);CHKERRQ(ierr);
  nbodyG(nb,xm,ym,yl,nb,xm,ym,xl,1e-2);
  for (i=0; i<nb; i++) {
    yl[i] += (1 - log(re[i] / 2)) * xl[i];
  }
  ierr = VecRestoreArray(x,&xl);CHKERRQ(ierr);
  ierr = VecRestoreArray(y,&yl);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCSetUp_FMM"
static PetscErrorCode PCSetUp_FMM(PC pc)
{
  PC_FMM *fmm = (PC_FMM*)pc->data;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  if (!fmm->iga) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCFMMSetIGA() must be called first");
  ierr = IGACreateNodeDM(fmm->iga,1,&(fmm->da));CHKERRQ(ierr);
  DMDALocalInfo info;
  ierr = DMDAGetLocalInfo(fmm->da,&info);CHKERRQ(ierr);
  fmm->nb = 0;
  if(info.xs == 0){
    fmm->nb += info.ym;
    if(info.ys == 0)               fmm->nb -= 1;
    if(info.my == info.ys+info.ym) fmm->nb -= 1;
  }
  if(info.mx == info.xs+info.xm){
    fmm->nb += info.ym;
    if(info.ys == 0)               fmm->nb -= 1;
    if(info.my == info.ys+info.ym) fmm->nb -= 1;
  }
  if(info.ys == 0              ) fmm->nb += info.xm;
  if(info.my == info.ys+info.ym) fmm->nb += info.xm;
  fmm->ni = info.xm*info.ym - fmm->nb;

  ierr = PetscMalloc(fmm->nb*sizeof(PetscInt),&(fmm->ip1));CHKERRQ(ierr);
  ierr = PetscMalloc3(fmm->ni,PetscReal,&(fmm->xi),
		      fmm->ni,PetscReal,&(fmm->yi),
		      fmm->ni,PetscReal,&(fmm->ri));CHKERRQ(ierr);
  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xb),
		      fmm->nb,PetscReal,&(fmm->yb),
		      fmm->nb,PetscReal,&(fmm->ub));CHKERRQ(ierr);
  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->xm),
		      fmm->nb,PetscReal,&(fmm->ym),
		      fmm->nb,PetscReal,&(fmm->um));CHKERRQ(ierr);
  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->dxe),
		      fmm->nb,PetscReal,&(fmm->dye),
		      fmm->nb,PetscReal,&(fmm->re));CHKERRQ(ierr);
  ierr = PetscMalloc3(fmm->nb,PetscReal,&(fmm->rhs),
		      fmm->nb,PetscReal,&(fmm->un),
		      fmm->ni,PetscReal,&(fmm->ui));CHKERRQ(ierr);
  ierr = VecCreateMPI(PETSC_COMM_WORLD,fmm->nb,PETSC_DETERMINE,&(fmm->u));CHKERRQ(ierr);
  ierr = VecDuplicate(fmm->u,&(fmm->f));CHKERRQ(ierr);
  ierr = MatCreate(PETSC_COMM_WORLD,&fmm->G);CHKERRQ(ierr);
  ierr = MatSetSizes(fmm->G,fmm->nb,fmm->nb,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
#if DENSE
  ierr = MatSetType(fmm->G,MATDENSE);CHKERRQ(ierr);
#else
  ierr = MatSetType(fmm->G,MATSHELL);CHKERRQ(ierr);
  ierr = MatShellSetOperation(fmm->G,MATOP_MULT,(void(*)(void))BEM_MatMult);CHKERRQ(ierr);
  ierr = MatShellSetContext(fmm->G,fmm);CHKERRQ(ierr);
#endif
  ierr = MatSetFromOptions(fmm->G);CHKERRQ(ierr);
  ierr = MatSetUp(fmm->G);CHKERRQ(ierr);

  PetscReal eps = .2 / (1<<((int)(log2(info.mx-1))));
  FMM_Init(eps*eps, 0);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCApply_FMM"
static PetscErrorCode PCApply_FMM(PC pc, Vec x, Vec y)
{
  PC_FMM *fmm = (PC_FMM*)pc->data;
  PetscErrorCode ierr;
  PetscFunctionBegin;

  DM da = fmm->da;

  PetscScalar **X,**Y;
  ierr = DMDAVecGetArray(da,x,&X);CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,y,&Y);CHKERRQ(ierr);
  DMDALocalInfo info;
  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);

  Vec u = fmm->u;
  Vec f = fmm->f;
  Mat G = fmm->G;
  IGA iga = fmm->iga;
  PetscInt ni = fmm->ni;
  PetscInt nb = fmm->nb;
  PetscReal *xi  = fmm->xi;
  PetscReal *yi  = fmm->yi;
  PetscReal *ri  = fmm->ri;
  PetscReal *ub  = fmm->ub;
  PetscReal *xm  = fmm->xm;
  PetscReal *ym  = fmm->ym;
  PetscReal *um  = fmm->um;
  PetscReal *dxe = fmm->dxe;
  PetscReal *dye = fmm->dye;
  PetscReal *re  = fmm->re;
  PetscReal *rhs = fmm->rhs;
  PetscReal *un  = fmm->un;
  PetscReal *ui  = fmm->ui;

  PetscInt i, j, ic = 0;
  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
    for (j=info.ys; j<info.ys+info.ym; j++) {
      if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1) {
	xi[ic] = (PetscReal)i / ( (PetscReal)(info.mx-1) );
	yi[ic] = (PetscReal)j / ( (PetscReal)(info.my-1) );
	ri[ic] = X[j][i] / 2 / M_PI;
	ic++;
      }
    }
  }
  assert( ni == ic );
  PetscReal eps = .2 / (1<<((int)(log2(info.mx-1))));

  ic = 0;
  j = 0;
  PetscInt q;
  PetscInt *start = iga->elem_start;
  PetscInt *width = iga->elem_width;
  PetscInt *sizes = iga->elem_sizes;
  IGABasis *BD    = iga->basis;
  PetscReal val_dirichlet = 0;
  PetscReal h     =  BD[0]->detJ[0]*2;
  for (j=0; j<2; j++){ // loop bottom and top
    if (start[1]+j*width[1] == j*sizes[1]) { // domain boundary?
      for (i=start[0]; i<start[0]+width[0]; i++) {
	for (q=0; q<1; q++) {
	  xm[ic]  =  (((PetscReal)i)+0.5)*h; //BD[0]->point[i*BD[0]->nqp+q];
	  ym[ic]  =  j*1;
	  um[ic]  = -val_dirichlet*0.5/M_PI;
	  dxe[ic] =  BD[0]->detJ[i]*2 *(j == 0 ? 1 : -1);
	  dye[ic] =  0;
	  re[ic]  =  fabs(dxe[ic]);
	  rhs[ic] =  um[ic]/re[ic]*M_PI;
	  //printf("%d%d %f %f %f %f %f %f %f\n",start[0],start[1],xm[ic],ym[ic],um[ic],dxe[ic],dye[ic],re[ic],rhs[ic]);
	  ic += 1;
	}
      }
    }
  }
  for (j=0; j<2; j++){ // loop left and right
    if (start[0]+j*width[0] == j*sizes[0]) { // domain boundary?
      for (i=start[1]; i<start[1]+width[1]; i++) {
	for (q=0; q<1; q++) {
	  xm[ic]  =  j*1;
	  ym[ic]  =  (((PetscReal)i)+0.5)*h; //BD[1]->point[i*BD[1]->nqp+q];
	  um[ic]  = -val_dirichlet*0.5/M_PI;
	  dxe[ic] =  0;
	  dye[ic] =  BD[1]->detJ[i]*2 *(j == 0 ? -1 : 1);
	  re[ic]  =  fabs(dye[ic]);
	  rhs[ic] =  um[ic]/re[ic]*M_PI;
	  //printf("%d%d %f %f %f %f %f %f %f\n",start[0],start[1],xm[ic],ym[ic],um[ic],dxe[ic],dye[ic],re[ic],rhs[ic]);
	  ic += 1;
	}
      }
    }
  }

  nb = ic; // this is ok, but we should improve memory allocation
  assert( nb == ic );

  //FMM(nb,xm,ym,rhs,ni,xi,yi,ri);
  nbodyG(nb,xm,ym,rhs,ni,xi,yi,ri,eps);
  //nbodyGn(nb,xm,ym,rhs,nb,xm,ym,um,dxe,dye,re,eps);

  for (i=0; i<nb; i++) {
    VecSetValue(f,i,rhs[i],INSERT_VALUES);
#if DENSE
    for (j=0; j<i; j++) {
      PetscReal dx = xm[i] - xm[j];
      PetscReal dy = ym[i] - ym[j];
      PetscReal r = sqrt(dx * dx + dy * dy + 1e-4);
      PetscReal value = -log(r);
      MatSetValue(G,i,j,value,INSERT_VALUES);
      MatSetValue(G,j,i,value,INSERT_VALUES);
    }
    PetscReal value = 1-log(re[i]/2);
    MatSetValue(G,i,i,value,INSERT_VALUES);
#endif
  }
  ierr = VecAssemblyBegin(f);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(f);CHKERRQ(ierr);
#if DENSE
  ierr = MatAssemblyBegin(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(G,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
#endif

  KSP ksp;
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,G,G,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSetOptionsPrefix(ksp,"pc_fmm_");CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  ierr = KSPSolve(ksp,f,u);CHKERRQ(ierr);
  ierr = VecGetArray(u,&un);CHKERRQ(ierr);

  for (i=0; i<nb; i++) {
    un[i] = -un[i];
  }

  for (i=0; i<ni; i++) {
    ui[i] = 0;
  }
  //FMM(ni,xi,yi,ui,ni,xi,yi,ri);
  //FMM(ni,xi,yi,ui,nb,xm,ym,un);
  nbodyG(ni,xi,yi,ui,ni,xi,yi,ri,eps);
  nbodyG(ni,xi,yi,ui,nb,xm,ym,un,eps);
  //nbodyGn(ni,xi,yi,ui,nb,xm,ym,um,dxe,dye,re,eps);

  ic = 0;
  for (i=info.xs; i<info.xs+info.xm; i++) { // Loop over interior points
    for (j=info.ys; j<info.ys+info.ym; j++) {
      if (i != 0 && i != info.mx-1 && j != 0 && j != info.my-1) {
        Y[j][i] = ui[ic];
        ic++;
      } else {
        Y[j][i] = 0;
      }
    }
  }

#if 0
  ic = 0;
  j = 0;
  for (i=info.xs; i<info.xs+info.xm-1; i++,ic++) { // Bottom points
    Y[j][i] = ub[ic];
  }
  i = info.xs+info.xm-1;
  for (j=info.ys; j<info.ys+info.ym-1; j++,ic++) { // Right points
    Y[j][i] = ub[ic];
  }
  j = info.ys+info.ym-1;
  for (i=info.xs+info.xm-1; i>info.xs; i--,ic++) { // Top points
    Y[j][i] = ub[ic];
  }
  i = 0;
  for (j=info.ys+info.ym-1; j>info.ys; j--,ic++) { // Left points
    Y[j][i] = ub[ic];
  }
#endif

  ierr = DMDAVecRestoreArray(da,x,&X);CHKERRQ(ierr);
  ierr = DMDAVecRestoreArray(da,y,&Y);CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCView_FMM"
static PetscErrorCode PCView_FMM(PC pc, PetscViewer viewer)
{
  PC_FMM *fmm = (PC_FMM*)pc->data;
  PetscBool      isascii;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
  if (!isascii) PetscFunctionReturn(0);
  if (!fmm->da) PetscFunctionReturn(0);
  ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
  ierr = PetscViewerASCIIPrintf(viewer,"DAView for dof grid:\n");CHKERRQ(ierr);
  ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
  ierr = DMView(fmm->da,viewer);CHKERRQ(ierr);
  ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
  ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCReset_FMM"
static PetscErrorCode PCReset_FMM(PC pc)
{
  PC_FMM *fmm = (PC_FMM*)pc->data;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = DMDestroy(&(fmm->da));CHKERRQ(ierr);
  ierr = VecDestroy(&(fmm->u));CHKERRQ(ierr);
  ierr = VecDestroy(&(fmm->f));CHKERRQ(ierr);
  ierr = MatDestroy(&(fmm->G));CHKERRQ(ierr);
  ierr = PetscFree3(fmm->xi,fmm->yi,fmm->ri);CHKERRQ(ierr);
  ierr = PetscFree3(fmm->xb,fmm->yb,fmm->ub);CHKERRQ(ierr);
  ierr = PetscFree3(fmm->xm,fmm->ym,fmm->um);CHKERRQ(ierr);
  ierr = PetscFree3(fmm->dxe,fmm->dye,fmm->re);CHKERRQ(ierr);
  ierr = PetscFree3(fmm->rhs,fmm->un,fmm->ui);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCDestroy_FMM"
static PetscErrorCode PCDestroy_FMM(PC pc)
{
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = PCReset_FMM(pc);CHKERRQ(ierr);
  ierr = PetscFree(pc->data);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

EXTERN_C_BEGIN
#undef  __FUNCT__
#define __FUNCT__ "PCCreate_FMM"
PetscErrorCode PCCreate_FMM(PC pc)
{
  PC_FMM *fmm = 0;
  PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = PetscNewLog(pc,PC_FMM,&fmm);CHKERRQ(ierr);
  pc->data = (void*)fmm;

  pc->ops->setup               = PCSetUp_FMM;
  pc->ops->reset               = PCReset_FMM;
  pc->ops->destroy             = PCDestroy_FMM;
  pc->ops->setfromoptions      = 0;/* PCSetFromOptions_FMM; */
  pc->ops->view                = PCView_FMM;
  pc->ops->apply               = PCApply_FMM;
  pc->ops->applytranspose      = 0;/* PCApplyTranspose_FMM; */

  PetscFunctionReturn(0);
}
EXTERN_C_END