bmatrix.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "field/field.h"
#include "init/init.h"
#include "util/util.h"

int VertBSA(GridData *gdata,Vec vertbsa)
{
  MagparFunctionLogBegin;

  /* Calculate solid angles of boundary vertices */

  /* loop over all local elements */
  for (int i=0; i<gdata->ln_ele; i++) {
    /* loop over all vertices */
    for (int j=0; j<NV; j++) {
      PetscReal  t_omega;
      ierr = SolidAngle(
        gdata->vertxyz+gdata->elevert[NV*i+j]*ND,
        gdata->vertxyz+gdata->elevert[NV*i+(j+1+NV)%NV]*ND,
        gdata->vertxyz+gdata->elevert[NV*i+(j+2+NV)%NV]*ND,
        gdata->vertxyz+gdata->elevert[NV*i+(j+3+NV)%NV]*ND,
        &t_omega
      );CHKERRQ(ierr);

      int v;
      v=gdata->elevert[NV*i+j];
      ierr = VecSetValue(
        vertbsa,
        v,
        t_omega,
      ADD_VALUES);CHKERRQ(ierr);
    }
  }

  ierr = VecAssemblyBegin(vertbsa);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(vertbsa);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


int BMatrix(GridData *gdata,Mat *Bmat)
{
  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Using dense boundary matrix\n");

  PetscPrintf(PETSC_COMM_WORLD,
    "Size of dense B matrix: %i x %i = %g matrix elements = %g MB.\n",
    gdata->n_vert_bnd,
    gdata->n_vert_bnd,
    1.0*gdata->n_vert_bnd*gdata->n_vert_bnd,
    (gdata->n_vert_bnd/1024.0*gdata->n_vert_bnd/1024.0*sizeof(PetscReal))
  );

  /* create boundary matrix (boundary nodes only) */

  int ln_bmatrixrows,rowstart,rowend;
  ln_bmatrixrows=gdata->n_vert_bnd/size;
  rowstart=rank*(gdata->n_vert_bnd/size);
  rowend=rowstart+gdata->n_vert_bnd/size-1;
  if ((gdata->n_vert_bnd%size)/(rank+1)>=1) {
    ln_bmatrixrows++;
    rowstart+=rank;
    rowend+=rank+1;
  }
  else {
    rowstart+=gdata->n_vert_bnd%size;
    rowend+=gdata->n_vert_bnd%size;
  }
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
    "<%i>local rows of boundary matrix: %i (from %i to %i)\n",
    rank,ln_bmatrixrows,
    rowstart,rowend
  );
  PetscSynchronizedFlush(PETSC_COMM_WORLD);
  assert(ln_bmatrixrows==rowend-rowstart+1);

  Mat B;

  PetscPrintf(PETSC_COMM_WORLD,
    "Allocating %i MB of memory for boundary matrix on each processor (total %i MB)\n",
    (gdata->n_vert_bnd/1024*gdata->n_vert_bnd/1024*sizeof(PetscReal)+1)/size+1,
    (gdata->n_vert_bnd/1024*gdata->n_vert_bnd/1024*sizeof(PetscReal)+2)
  );

  ierr = MatCreateMPIDense(PETSC_COMM_WORLD,
    ln_bmatrixrows,ln_bmatrixrows,
    gdata->n_vert_bnd,gdata->n_vert_bnd,
    PETSC_NULL,
    &B);CHKERRQ(ierr);
  ierr = MatSetFromOptions(B);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Calculating boundary matrix...\n");

  int zeroes=0;

  /*
  hybrid FEM/BEM method from:
  T. R. Koehler, D. R. Fredkin,
  "Finite Element Methods for Micromagnetics"
  IEEE Trans. Magn. 28 (1992) 1239-1244

  matrix elements calculated according to:
  D. A. Lindholm,
  "Three-Dimensional Magnetostatic Fields from
  Point-Matched Integral Equations with Linearly Varying Scalar Sources",
  IEEE Trans. Magn. MAG-20 (1984) 2025-2032.
  */

  /* loop over all (!) vertices (R) */
  int i;
  for (i=0; i<gdata->n_vert; i++) {
    ierr = ProgressBar(i,gdata->n_vert,10);
    /* skip if R is not a boundary node */
    if (gdata->vertbndg2bnd[i]==C_INT) continue;

    /* loop over local surface triangles */

    for (int j=0; j<gdata->ln_bnd_fac; j++) {
      /* skip if R is a node of j */
      if (gdata->bndfacvert[j*NN+0]==i ||
          gdata->bndfacvert[j*NN+1]==i ||
          gdata->bndfacvert[j*NN+2]==i
      ) continue;

      PetscReal matele[ND];
      ierr = Bele(
                  gdata->vertxyz+ND*i,
                  gdata->vertxyz+ND*gdata->bndfacvert[j*NN+0],
                  gdata->vertxyz+ND*gdata->bndfacvert[j*NN+1],
                  gdata->vertxyz+ND*gdata->bndfacvert[j*NN+2],
                  matele
      );CHKERRQ(ierr);

      if (PetscAbsReal(matele[0])+PetscAbsReal(matele[1])+PetscAbsReal(matele[2])<D_EPS)
        continue;

      for (int k=0;k<NN;k++) {
        /* set (dense) full boundary matrix */
        ierr = MatSetValue(B,
          gdata->vertbndg2bnd[i],
          gdata->vertbndg2bnd[gdata->bndfacvert[j*NN+k]],
          matele[k],
        ADD_VALUES);CHKERRQ(ierr);
      }
    }

    /* flush caches to save memory (and avoid swapping) */
    if (i % 100 == 99) {
      ierr = MatAssemblyBegin(B,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
      ierr = MatAssemblyEnd(B,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
    }

  }
  ierr = ProgressBar(1,1,10);

  if (size>1) {
    ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"<%i> took %g s\n",rank,t_t2-t_t1);
    PetscPrintf(PETSC_COMM_WORLD,"waiting for other processors to complete...");
  }

  /* continue assembling to stay in sync with other processors */
  for (;i<gdata->n_vert;i++) {
    if (i % 100 == 99) {
      ierr = MatAssemblyBegin(B,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
      ierr = MatAssemblyEnd(B,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
    }
  }

  /* TODO not quite accurate */
  int tzeroes=0;
  ierr=MPI_Allreduce((int*)&zeroes,(int*)&tzeroes,1,MPI_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
  PetscInfo3(0,
    "total zero entries: %i (%i) of %i\n",
    (tzeroes+gdata->n_vert_bnd)/3,
    tzeroes,
    gdata->n_vert_bnd*gdata->n_vert_bnd
  );

  ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,
    "boundary matrix calculation took %g s = %g min\n",
    t_t2-t_t1,
    (t_t2-t_t1)/60.0
  );

  /* set diagonal matrix elements */

  Vec vertbsa; /* boundary solid angle (4*pi for interior vertices) */
  ierr = VecCreate(PETSC_COMM_WORLD,&vertbsa);CHKERRQ(ierr);
  ierr = VecSetSizes(vertbsa,gdata->ln_vert,gdata->n_vert);CHKERRQ(ierr);
  ierr = VecSetFromOptions(vertbsa);CHKERRQ(ierr);

  ierr = VertBSA(gdata,vertbsa);CHKERRQ(ierr);

  PetscReal *t_array;
  ierr = VecGetArray(vertbsa,&t_array);CHKERRQ(ierr);

  for (i=0; i<gdata->ln_vert; i++) {

    /* skip if R is not a boundary node */
    if (gdata->vertbndg2bnd[gdata->vertl2g[i]]==C_INT) continue;

    PetscReal matele;
    matele=t_array[i]/(4.0*PETSC_PI)-1.0;
    ierr = MatSetValue(B,
      gdata->vertbndg2bnd[gdata->vertl2g[i]],
      gdata->vertbndg2bnd[gdata->vertl2g[i]],
      matele,
    ADD_VALUES);CHKERRQ(ierr);
  }

  ierr = VecRestoreArray(vertbsa,&t_array);CHKERRQ(ierr);
  ierr = VecDestroy(vertbsa);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"B matrix assembly complete\n");
  ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = PrintMatInfoAll(B);CHKERRQ(ierr);

  /* scale B*=-1 since Au2*=-1 and ABu2*=1 is not */
  ierr = MatScale(B,-1.0);CHKERRQ(ierr);

  *Bmat=B;

/* DEBUG
  PetscPrintf(PETSC_COMM_WORLD,"B matrix::\n");
  ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
*/

/* DEBUG
  ierr = PetscViewerSetFormat(PETSC_VIEWER_DRAW_WORLD,PETSC_VIEWER_ASCII_DENSE);CHKERRQ(ierr);
  ierr = MatView(B,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);

  char fmesh[256];
  PetscViewer viewer;

  ierr = PetscSNPrintf(fmesh,255,"%s%s",gdata->simname,".ba");CHKERRQ(ierr);
  ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,fmesh,&viewer);CHKERRQ(ierr);
  ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_ASCII_DENSE);CHKERRQ(ierr);
  ierr = MatView(B,viewer);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
*/

  ierr = PetscGetTime(&t_t2);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}