hdemag.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 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
*/

static char const Id[] = "$Id: hdemag.c 2897 2009-12-04 03:56:51Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

/*
     Calculate magnetostatic field using a hybrid FE/BE method:

     D. R. Fredkin, T. R. Koehler,
     "Hybrid Method for Computing Demagnetizing Fields"
     IEEE Trans. Magn. 26 (1990) 415-417
*/

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

static int doinit=1;
static int fieldon=0;

static Mat Bmat;     
static Mat AdivM;    
static Mat Agradu;   
static Mat ABu1;     
static Mat ABu2;     
static Vec u2bnd;    
static Vec u1;       
static Vec u1bnd;    
static Vec u2;       
static Vec bu1;      
static Vec bu2;      
static Vec utot;     
KSP ksp_Au1;
KSP ksp_Au2;

#ifdef HDEMAGUEPOL
static Vec u1bak; 
static Vec u2bak; 
static Vec du1dt; 
static Vec du2dt; 
#endif

#include "bmatrix.c"

int Hdemag_Init(GridData *gdata)
{
  MagparFunctionLogBegin;

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

  gdata->VHdem=PETSC_NULL;

  PetscTruth flg;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-demag",(PetscInt*)&fieldon,&flg);CHKERRQ(ierr);
  if (flg!=PETSC_TRUE) {
    fieldon=1;
    PetscPrintf(PETSC_COMM_WORLD,
      "Option -demag not found, using default value: %i\n",
      fieldon
    );
  }

  if (!fieldon) {
    PetscPrintf(PETSC_COMM_WORLD,"Hdemag off\n");
    MagparFunctionLogReturn(0);
  }
  PetscPrintf(PETSC_COMM_WORLD,"Hdemag on\n");

/* activate this define to enforce Dirichlet boundary condition
   on one node (id given by U1DIRICHLET) for u1
*/
#if 0
#define U1DIRICHLET 0
  PetscPrintf(PETSC_COMM_WORLD,"Using 1 Dirichlet BC for u1 on node %i\n",U1DIRICHLET);
#else
  PetscPrintf(PETSC_COMM_WORLD,"No Dirichlet BC for u1 (floating)\n");
#endif

  /* Create stiffness matrix */

  /* Creates a sparse parallel matrix in AIJ format - speedup: 50 !!! */
  /* sensible values for preallocated memory !? */

  int ln_bmatrixrows;
  ln_bmatrixrows=gdata->n_vert_bnd/size;
  if ((gdata->n_vert_bnd%size)/(rank+1)>=1) ln_bmatrixrows++;

  int ln_vert;
  ln_vert=gdata->ln_vert;
  int n_vert;
  n_vert=gdata->n_vert;

  PetscPrintf(PETSC_COMM_WORLD,"Calculating sparsity pattern for matrix preallocation\n");

  /* get ownership range */
  int low,high;
  ierr = VecGetOwnershipRange(gdata->elevol,&low,&high);CHKERRQ(ierr);

  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
    "<%i>element ownership: %i - %i\n",
    rank,low,high-1
  );
  PetscSynchronizedFlush(PETSC_COMM_WORLD);

  int *elevertall;

#ifdef HAVE_ELEVERTALL
  elevertall=gdata->elevertall;
#else
  if (size>1) {
    ierr = PetscMalloc(NV*gdata->n_ele*sizeof(int),&elevertall);CHKERRQ(ierr);
    ierr = PetscMemcpy(elevertall+NV*low,gdata->elevert,NV*gdata->ln_ele*sizeof(int));CHKERRQ(ierr);

    int *p;
    p=elevertall;

    int recvcount[size];
    recvcount[rank]=NV*gdata->ln_ele;
    for (int i=0;i<size;i++) {
      ierr = MPI_Bcast(recvcount+i,1,MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      /* we could also use MPI_Send/Recv to get everything just to the first processor */
      ierr = MPI_Bcast(p,recvcount[i],MPI_INT,i,PETSC_COMM_WORLD);CHKERRQ(ierr);
      p+=recvcount[i];
    }
    assert(p==elevertall+NV*gdata->n_ele);
  }
  else {
    elevertall=gdata->elevert;
  }
#endif

  int *ia,*ja;
#ifdef HAVE_M2IJ
  ia=gdata->mesh2nodal_ia;
  ja=gdata->mesh2nodal_ja;
#else
  ierr = Mesh2Nodal(gdata->n_ele,gdata->n_vert,elevertall,&ia,&ja);CHKERRQ(ierr);
#endif

  /* get ownership range */
  ierr = VecGetOwnershipRange(gdata->M,&low,&high);CHKERRQ(ierr);
  low/=ND;
  high/=ND;

  /* create utot here, so we can use VecGetOwnershipRange below */
  ierr = VecCreate(PETSC_COMM_WORLD,&utot);CHKERRQ(ierr);
  ierr = VecSetSizes(utot,ln_vert,n_vert);CHKERRQ(ierr);
  ierr = VecSetFromOptions(utot);CHKERRQ(ierr);

#define XND 1
  int *d_nnz;
  int *o_nnz;
  ierr = PetscMalloc(ND*gdata->n_vert*sizeof(int),&d_nnz);CHKERRQ(ierr);
  ierr = PetscMalloc(ND*gdata->n_vert*sizeof(int),&o_nnz);CHKERRQ(ierr);
  /* calculate number of nonzeroes for each line */
  for (int i=0;i<gdata->n_vert;i++) {
    for (int j=0;j<XND;j++) {
      d_nnz[XND*i+j]=XND;
      o_nnz[XND*i+j]=0;
      for (int l=ia[i];l<ia[i+1];l++) {
        if (ja[l]>=low && ja[l]<high) {
          d_nnz[XND*i+j]+=1;
        }
        else {
          o_nnz[XND*i+j]+=1;
        }
      }
    }
  }

  int *d_nnzm;
  int *o_nnzm;
  ierr = PetscMalloc(ND*n_vert*sizeof(int),&d_nnzm);CHKERRQ(ierr);
  ierr = PetscMalloc(ND*n_vert*sizeof(int),&o_nnzm);CHKERRQ(ierr);

  /* connectivity matrix for merged mesh */
  int lowm,highm;

    lowm=low;
    highm=high;

    ierr = PetscMemcpy(d_nnzm,d_nnz,ND*n_vert*sizeof(int));CHKERRQ(ierr);
    ierr = PetscMemcpy(o_nnzm,o_nnz,ND*n_vert*sizeof(int));CHKERRQ(ierr);

  /* conventional MPIAIJ matrix */
  Mat Au1;
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ln_vert,ln_vert,
    n_vert,n_vert,
    0,d_nnzm+lowm,0,o_nnzm+lowm,
    &Au1
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(Au1);CHKERRQ(ierr);

  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ln_bmatrixrows,ln_vert,
    gdata->n_vert_bnd,n_vert,
    1,PETSC_NULL,1,PETSC_NULL,
    &ABu1
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(ABu1);CHKERRQ(ierr);

  Mat Au2;
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ln_vert,ln_vert,
    n_vert,n_vert,
    0,d_nnzm+lowm,0,o_nnzm+lowm,
    &Au2
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(Au2);CHKERRQ(ierr);

  /* preallocate with d_nnz for diag. and off-diag. because of boundary vertex columns */

  for (int i=lowm; i<highm; i++) {
    d_nnzm[i] = PetscMin(d_nnzm[i],ln_bmatrixrows);
    o_nnzm[i] = PetscMin(d_nnzm[i],gdata->n_vert_bnd-ln_bmatrixrows);
    assert(d_nnzm[i]>=0);
    assert(o_nnzm[i]>=0);
  }

  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ln_vert,ln_bmatrixrows,
    n_vert,gdata->n_vert_bnd,
    0,d_nnzm+lowm,0,o_nnzm+lowm,
    &ABu2
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(ABu2);CHKERRQ(ierr);


  for (int i=0; i<n_vert; i++) {
    d_nnzm[i] = 0;
    o_nnzm[i] = 0;
  }
  for (int i=0; i<gdata->n_vert; i++) {
    int i2;
    i2=i;
      if (i<lowm || i>=highm) continue;
    assert(d_nnz[i]>=0);
    assert(o_nnz[i]>=0);
    d_nnzm[i2] += ND*d_nnz[i]+ND*o_nnz[i];
    o_nnzm[i2] += ND*d_nnz[i]+ND*o_nnz[i];
    assert(i2<n_vert);
    assert(d_nnzm[i2]>=0);
    assert(o_nnzm[i2]>=0);
  }

  /* create matrix to calculate div M */
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ln_vert,ND*gdata->ln_vert,
    n_vert,ND*gdata->n_vert,
    0,d_nnzm+lowm,0,o_nnzm+lowm,
    &AdivM
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(AdivM);CHKERRQ(ierr);

  for (int i=ND*gdata->n_vert-1; i>=0; i--) {
    d_nnz[i] = d_nnz[i/ND]+o_nnz[i/ND];
    o_nnz[i] = d_nnz[i/ND]+o_nnz[i/ND];
    assert(d_nnz[i]>=0);
    assert(o_nnz[i]>=0);
  }

  /* create matrix to calculate Hdem = -grad u */
  ierr = MatCreateMPIAIJ(
    PETSC_COMM_WORLD,
    ND*gdata->ln_vert,ln_vert,
    ND*gdata->n_vert,n_vert,
    0,d_nnz+ND*low,0,o_nnz+ND*low,
    &Agradu
  );CHKERRQ(ierr);
  ierr = MatSetFromOptions(Agradu);CHKERRQ(ierr);

#ifndef HAVE_M2IJ
  ierr = PetscFree(ia);CHKERRQ(ierr);
  ierr = PetscFree(ja);CHKERRQ(ierr);
#endif
  ierr = PetscFree(d_nnz);CHKERRQ(ierr);
  ierr = PetscFree(o_nnz);CHKERRQ(ierr);
#ifndef HAVE_ELEVERTALL
  if (size>1) {
    ierr = PetscFree(elevertall);CHKERRQ(ierr);
  }
#endif
  ierr = PetscFree(d_nnzm);CHKERRQ(ierr);
  ierr = PetscFree(o_nnzm);CHKERRQ(ierr);

  PetscReal *ta_elevol;
  ierr = VecGetArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);

  PetscPrintf(PETSC_COMM_WORLD,"Calculating stiffness matrices...\n");
  for (int i=0; i<gdata->ln_ele; i++) {
    ierr = ProgressBar(i,gdata->ln_ele,10);

    /* get vertex ids */
    int t_v[NV];
    int t_w[NV];

    for (int j=0; j<NV; j++) {
      t_v[j]=gdata->elevert[i*NV+j];
      t_w[j]=t_v[j];
    }

    /* calculate gradients */
    PetscReal D_etaj[NV][ND];
    tetgrad(
      gdata->vertxyz+t_v[0]*ND,
      gdata->vertxyz+t_v[1]*ND,
      gdata->vertxyz+t_v[2]*ND,
      gdata->vertxyz+t_v[3]*ND,
      D_etaj
    );

    /* set matrix elements */

    for (int l=0;l<NV;l++) {
      for (int j=0;j<NV;j++) {
        /* calculate stiffness matrix element */
        PetscReal matele;
        matele=-ta_elevol[i]*my_ddot(ND,D_etaj[j],1,D_etaj[l],1);

        /* inserted many of these if clauses to avoid zero entries
           and save memory
        */
        if (PetscAbsReal(matele)>D_EPS) {
#ifdef U1DIRICHLET
          /* skip first row and first column, where we impose the
             Dirichlet BC: u1=0
          */
          if (t_w[j] != U1DIRICHLET &&
              t_w[l] != U1DIRICHLET) {
#endif
            ierr = MatSetValue(
              Au1,
              t_w[j],
              t_w[l],
              matele,
              ADD_VALUES
            );CHKERRQ(ierr);

            assert(t_w[j]<n_vert);
            assert(t_w[l]<n_vert);
#ifdef U1DIRICHLET
          }
#endif

          int tvj,tvl;
          tvj=t_v[j];
          tvl=t_v[l];
          if (gdata->vertbndg2bnd[tvj]==C_INT) {
            /* set non-Dirichlet rows (i.e. for interior nodes) only */
            if (gdata->vertbndg2bnd[tvl]==C_INT) {
              /* set non-Dirichlet columns only to keep Au2 symmetric
               * columns eliminated and worked into rhs through ABu2 (see below)
               */
              ierr = MatSetValue(
                Au2,
                t_w[j],
                t_w[l],
                matele,
                ADD_VALUES
              );CHKERRQ(ierr);
              assert(t_w[j]>=0);
              assert(t_w[l]>=0);
              assert(t_w[j]<n_vert);
              assert(t_w[l]<n_vert);
            }
            else {
              /* Dirichlet columns eliminated from Au2 and merged into RHS through ABu2 */
              ierr = MatSetValue(
                ABu2,
                t_w[j],
                gdata->vertbndg2bnd[tvl],
                -matele,
                ADD_VALUES
              );CHKERRQ(ierr);
              assert(t_w[j]<n_vert);
              assert(gdata->vertbndg2bnd[tvl]<gdata->n_vert_bnd);
            }
          }
        }
      } /* end for j loop */
    } /* end for l loop */

    for (int l=0;l<NV;l++) {
      for (int k=0;k<ND;k++) {
        PetscReal matele;
        matele=ta_elevol[i]/NV*D_etaj[l][k];
        if (PetscAbsReal(matele)>D_EPS) {
          for (int j=0;j<NV;j++) {
            ierr = MatSetValue(
              Agradu,
              ND*t_v[j]+k,
              t_w[l],
              matele,
              ADD_VALUES
            );CHKERRQ(ierr);
            assert(t_v[j]<gdata->n_vert);
            assert(t_w[l]<n_vert);

#ifdef U1DIRICHLET
            /* skip nodes with Dirichlet BC */
            if (t_w[l] != U1DIRICHLET) {
#endif
            /* skip non-magnetic elements */
            if (gdata->propdat[NP*gdata->eleprop[i]+4]>0.0) {
              ierr = MatSetValue(
                AdivM,
                t_w[l],
                ND*t_v[j]+k,
                matele*gdata->propdat[NP*gdata->eleprop[i]+4],
                ADD_VALUES
              );CHKERRQ(ierr);
              assert(t_w[l]<n_vert);
              assert(t_v[j]<gdata->n_vert);
            }
#ifdef U1DIRICHLET
            }
#endif
          }
        }
      }
    }
  }
  ierr = ProgressBar(1,1,10);

#undef MATFLUSH
#ifdef MATFLUSH
  ierr = MatAssemblyBegin(Au1,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Au1,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(Au2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Au2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(ABu2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(ABu2,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(Agradu,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Agradu,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyBegin(AdivM,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(AdivM,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
#endif

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

  ierr = VecRestoreArray(gdata->elevol,&ta_elevol);CHKERRQ(ierr);

    /* set diagonal matrix element for nodes with Dirichlet BC */
    if (!rank) {
#ifdef U1DIRICHLET
      ierr = MatSetValue(Au1,U1DIRICHLET,U1DIRICHLET,-1.0,ADD_VALUES);CHKERRQ(ierr);
#endif

      for (int i=0; i<gdata->n_vert; i++) {
        int j,k;
        j=i;
        k=j;
        if (gdata->vertbndg2bnd[j] >= 0) {
          assert(gdata->vertbndg2bnd[j]<gdata->n_vert_bnd);

          ierr = MatSetValue(Au2,k,k,-1.0,ADD_VALUES);CHKERRQ(ierr);
          ierr = MatSetValue(ABu2,k,gdata->vertbndg2bnd[j],-1.0,ADD_VALUES);CHKERRQ(ierr);
          ierr = MatSetValue(ABu1,gdata->vertbndg2bnd[j],k,1.0,INSERT_VALUES);CHKERRQ(ierr);
        }
      }
    }

#undef DOMATCOMPRESS
    PetscPrintf(PETSC_COMM_WORLD,"Au1 matrix assembly complete\n");
    ierr = MatAssemblyBegin(Au1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Au1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    /* TODO: MatCompress does not do anything */
#ifdef DOMATCOMPRESS
    ierr = MatCompress(Au1);CHKERRQ(ierr);
#endif
    ierr = PrintMatInfoAll(Au1);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"Au2 matrix assembly complete\n");
    ierr = MatAssemblyBegin(Au2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Au2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
#ifdef DOMATCOMPRESS
    ierr = MatCompress(Au2);CHKERRQ(ierr);
#endif
    ierr = PrintMatInfoAll(Au2);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"AdivM matrix assembly complete\n");
    ierr = MatAssemblyBegin(AdivM,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(AdivM,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
#ifdef DOMATCOMPRESS
    ierr = MatCompress(AdivM);CHKERRQ(ierr);
#endif
    ierr = PrintMatInfoAll(AdivM);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"Agradu matrix assembly complete\n");
    ierr = MatAssemblyBegin(Agradu,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Agradu,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
#ifdef DOMATCOMPRESS
    ierr = MatCompress(Agradu);CHKERRQ(ierr);
#endif
    ierr = PrintMatInfoAll(Agradu);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"ABu1 matrix assembly complete\n");
    ierr = MatAssemblyBegin(ABu1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(ABu1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
#ifdef DOMATCOMPRESS
    ierr = MatCompress(ABu1);CHKERRQ(ierr);
#endif
    ierr = PrintMatInfoAll(ABu1);CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD,"ABu2 matrix assembly complete\n");
    ierr = MatAssemblyBegin(ABu2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(ABu2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
#ifdef DOMATCOMPRESS
    ierr = MatCompress(ABu2);CHKERRQ(ierr);
#endif
    ierr = PrintMatInfoAll(ABu2);CHKERRQ(ierr);

  Vec vertvolrec3;
  ierr = VecDuplicate(gdata->M,&vertvolrec3);CHKERRQ(ierr);
  for (int i=0;i<ND;i++) {
    ierr = VecStrideScatter(gdata->vertvol,i,vertvolrec3,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecReciprocal(vertvolrec3);CHKERRQ(ierr);
  ierr = MatDiagonalScale(Agradu,vertvolrec3,PETSC_NULL);CHKERRQ(ierr);

  ierr = VecDestroy(vertvolrec3);CHKERRQ(ierr);

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

/* DEBUG
  PetscPrintf(PETSC_COMM_WORLD,"deb01: AdivM\n");
  ierr = MatView(AdivM,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb01: Au1\n");
  ierr = MatView(Au1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb01: ABu1\n");
  ierr = MatView(ABu1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb01: ABu2\n");
  ierr = MatView(ABu2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb01: Au2\n");
  ierr = MatView(Au2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb01: Agradu\n");
  ierr = MatView(Agradu,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
*/

    /* convert the stiffness matrix (if requested) */
#ifdef SUPERLU
    Mat Au1tmp;
    ierr = MatConvert(Au1,MATSUPERLU_DIST,MAT_INITIAL_MATRIX,&Au1tmp);CHKERRQ(ierr);
    ierr = MatDestroy(Au1);CHKERRQ(ierr);
    Au1=Au1tmp;
#elif defined(BLOCKSOLVE)
    Mat Au1tmp;
    ierr = MatConvert(Au1,MATMPIROWBS,MAT_INITIAL_MATRIX,&Au1tmp);CHKERRQ(ierr);
    ierr = MatDestroy(Au1);CHKERRQ(ierr);
    Au1=Au1tmp;
#endif

#ifdef SUPERLU
    Mat Au2tmp;
    ierr = MatConvert(Au2,MATSUPERLU_DIST,MAT_INITIAL_MATRIX,&Au2tmp);CHKERRQ(ierr);
    ierr = MatDestroy(Au2);CHKERRQ(ierr);
    Au2=Au2tmp;
#elif defined(BLOCKSOLVE)
    Mat Au2tmp;
    ierr = MatConvert(Au2,MATMPIROWBS,MAT_INITIAL_MATRIX,&Au2tmp);CHKERRQ(ierr);
    ierr = MatDestroy(Au2);CHKERRQ(ierr);
    Au2=Au2tmp;
#endif


#if PETSC_VERSION >= 300
    ierr = MatSetOption(Au1,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
    ierr = MatSetOption(Au2,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
    ierr = MatSetOption(Au1,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);
    ierr = MatSetOption(Au2,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);
#else
    ierr = MatSetOption(Au1,MAT_SYMMETRIC);CHKERRQ(ierr);
    ierr = MatSetOption(Au2,MAT_SYMMETRIC);CHKERRQ(ierr);
    ierr = MatSetOption(Au1,MAT_SYMMETRY_ETERNAL);CHKERRQ(ierr);
    ierr = MatSetOption(Au2,MAT_SYMMETRY_ETERNAL);CHKERRQ(ierr);
#endif

    ierr = VecDuplicate(utot,&u1);CHKERRQ(ierr);
    ierr = VecDuplicate(utot,&u2);CHKERRQ(ierr);

    ierr = VecDuplicate(utot,&bu1);CHKERRQ(ierr);
    ierr = VecDuplicate(utot,&bu2);CHKERRQ(ierr);

    ierr = VecCreate(PETSC_COMM_WORLD,&u1bnd);CHKERRQ(ierr);
    ierr = VecSetSizes(u1bnd,ln_bmatrixrows,gdata->n_vert_bnd);CHKERRQ(ierr);
    ierr = VecSetFromOptions(u1bnd);CHKERRQ(ierr);

    ierr = VecDuplicate(u1bnd,&u2bnd);CHKERRQ(ierr);


    /* set up VecScatter from u1 to u1bnd and u2bnd to u2 */

    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_Au1);CHKERRQ(ierr);
    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp_Au2);CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp_Au1,Au1,Au1,SAME_PRECONDITIONER);CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp_Au2,Au2,Au2,SAME_PRECONDITIONER);CHKERRQ(ierr);

    /* check if options are set old style
       FIXME: just for backward compatibility, remove in next release
    */
    ierr = PetscOptionsHasName(PETSC_NULL,"-ksp_type",&flg);CHKERRQ(ierr);
    if (flg==PETSC_TRUE) {
      PetscPrintf(PETSC_COMM_WORLD,"Warning: Options '-ksp_*' are deprecated! Please use -hdemag_u{1,2}_ksp_* instead!\n");
      ierr = KSPSetFromOptions(ksp_Au1);CHKERRQ(ierr);
      ierr = KSPSetFromOptions(ksp_Au2);CHKERRQ(ierr);
    }
    else { /* new style */

    /* evaluates
        PetscOptions "-hdemag_u1_ksp_type"
        PetscOptions "-hdemag_u1_ksp_atol"
        PetscOptions "-hdemag_u1_ksp_rtol"
        PetscOptions "-hdemag_u2_ksp_type"
        PetscOptions "-hdemag_u2_ksp_atol"
        PetscOptions "-hdemag_u2_ksp_rtol"
    */

    for (int i=1;i<=2;i++) {
      char str[256],ostr[256];

      /* by default use conjugate gradient solver */
      sprintf(ostr,"-hdemag_u%i_ksp_type",i);
      ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
      if (flg!=PETSC_TRUE) {
        sprintf(str,"cg");
        PetscOptionsSetValue(ostr,str);
        PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
      }

      if (size>1) {
        /* by default use BlockJacobi preconditioner */
        sprintf(ostr,"-hdemag_u%i_pc_type",i);
        ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
        if (flg!=PETSC_TRUE) {
          sprintf(str,"bjacobi");
          PetscOptionsSetValue(ostr,str);
          PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
        }

        /* by default use ICC preconditioner on subblocks */
        sprintf(ostr,"-hdemag_u%i_sub_pc_type",i);
        ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
        if (flg!=PETSC_TRUE) {
          sprintf(str,"icc");
          PetscOptionsSetValue(ostr,str);
          PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
        }

        /* by default use Manteuffel shift */
        sprintf(ostr,"-hdemag_u%i_sub_pc_factor_shift_positive_definite",i);
        ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
        if (flg!=PETSC_TRUE) {
          sprintf(str,"1");
          PetscOptionsSetValue(ostr,str);
          PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
        }
      }
      else {
        /* by default use ICC preconditioner */
        sprintf(ostr,"-hdemag_u%i_pc_type",i);
        ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
        if (flg!=PETSC_TRUE) {
          sprintf(str,"icc");
          PetscOptionsSetValue(ostr,str);
          PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
        }

        /* by default use Manteuffel shift */
        sprintf(ostr,"-hdemag_u%i_pc_factor_shift_positive_definite",i);
        ierr = PetscOptionsGetString(PETSC_NULL,ostr,str,256,&flg);CHKERRQ(ierr);
        if (flg!=PETSC_TRUE) {
          sprintf(str,"1");
          PetscOptionsSetValue(ostr,str);
          PetscPrintf(PETSC_COMM_WORLD,"Option %s not found, using default value %s!\n",ostr,str);
        }
      }
    }
    } /* new style */

    /* use prefix for options */
    ierr = KSPSetOptionsPrefix(ksp_Au1,"hdemag_u1_");CHKERRQ(ierr);
    /* possibly override defaults with settings from user defined options */
    ierr = KSPSetFromOptions(ksp_Au1);CHKERRQ(ierr);

    /* use prefix for options */
    ierr = KSPSetOptionsPrefix(ksp_Au2,"hdemag_u2_");CHKERRQ(ierr);
    /* possibly override defaults with settings from user defined options */
    ierr = KSPSetFromOptions(ksp_Au2);CHKERRQ(ierr);

    /* allow user to override tolerances of ksp_Au2 */
    ierr = PetscOptionsHasName(PETSC_NULL,"-hdemag_u2_ksp_rtol",&flg);CHKERRQ(ierr);
    /* set default values if undefined */
    if (flg!=PETSC_TRUE) {
      PetscReal rtol,atol,dtol;
      int maxits;
      ierr = KSPGetTolerances(ksp_Au1,&rtol,&atol,&dtol,&maxits);CHKERRQ(ierr);
      /* set stricter rtol for u2 */
      PetscPrintf(PETSC_COMM_WORLD,"Adjusting tolerances for ksp_Au2\n");
      ierr = KSPSetTolerances(ksp_Au2,rtol*1e-3,atol,dtol,maxits);CHKERRQ(ierr);
    }

    /* just to get correct info from KSPView below */
    ierr = KSPSetInitialGuessNonzero(ksp_Au1,PETSC_TRUE);CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(ksp_Au2,PETSC_TRUE);CHKERRQ(ierr);

    /* KSPSetUp only necessary to get more useful output from KSPView */
    ierr = KSPSetUp(ksp_Au1);CHKERRQ(ierr);
    ierr = KSPSetUp(ksp_Au2);CHKERRQ(ierr);

#if defined(SUPERLU) || defined(BLOCKSOLVE)
    ierr = KSPGetTolerances(ksp_Au1,&rtol,&atol,&dtol,&maxits);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"KSP for A*u1=divM:\n");
    PetscPrintf(PETSC_COMM_WORLD,"rtol: %g atol: %g dtol: %g maxits: %i\n",
      rtol,atol,dtol,maxits
    );
    ierr = KSPGetTolerances(ksp_Au2,&rtol,&atol,&dtol,&maxits);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"KSP for Au2*u2=u2bnd:\n");
    PetscPrintf(PETSC_COMM_WORLD,"rtol: %g atol: %g dtol: %g maxits: %i\n",
      rtol,atol,dtol,maxits
    );
#else
    PetscPrintf(PETSC_COMM_WORLD,"KSP for A*u1=divM:\n");
    ierr = KSPView(ksp_Au1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"KSP for Au2*u2=u2bnd:\n");
    ierr = KSPView(ksp_Au2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
#endif

    ierr = VecDuplicate(gdata->M,&gdata->VHdem);CHKERRQ(ierr);

#ifdef HDEMAGUEPOL
  ierr = VecDuplicate(utot,&u1bak);CHKERRQ(ierr);
  ierr = VecDuplicate(utot,&u2bak);CHKERRQ(ierr);
  ierr = VecDuplicate(utot,&du1dt);CHKERRQ(ierr);
  ierr = VecDuplicate(utot,&du2dt);CHKERRQ(ierr);
  ierr = VecZeroEntries(du1dt);CHKERRQ(ierr);
  ierr = VecZeroEntries(du2dt);CHKERRQ(ierr);
#endif

  ierr = BMatrix(gdata,&Bmat);CHKERRQ(ierr);

  MagparFunctionLogReturn(0);
}


int Hdemag(GridData *gdata,Vec VHtotsum)
{
#ifdef HDEMAGUEPOL
  static PetscReal tubak=0.0;
  PetscReal dt;
#endif

  int its1,its2;

  MagparFunctionInfoBegin;

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

  if (doinit) {
    ierr = Hdemag_Init(gdata);CHKERRQ(ierr);
    doinit=0;
    /* reset timer to exclude initialization  */
    ierr = PetscGetTime(&t_t1);CHKERRQ(ierr);
  }
  if (!fieldon) {
    MagparFunctionInfoReturn(0);
  }

  PetscLogDouble t_t3;
  t_t3=t_t1;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         Compute the matrix and right-hand-side vector that define
         the linear system, Au = b.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /* set bu1 = div(M) */
  ierr = MatMult(AdivM,gdata->M,bu1);CHKERRQ(ierr);

#ifdef HDEMAGUEPOL
  dt=gdata->time-tubak;
  if (gdata->mode==0 && PetscAbsReal(dt) > 1e-8){
    /* extrapolate u1 (weakened) */
    PetscReal dt2;
    dt2=dt*0.5; /* TODO: good factor? */
    ierr = VecWAXPY(u1,dt2,du1dt,u1bak);CHKERRQ(ierr);
  }
  else {
    dt=0.0;
  }
#endif

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

  PetscTruth oprofile;
  PetscOptionsHasName(PETSC_NULL,"-profile",&oprofile);
  if (oprofile) {
    PetscPrintf(PETSC_COMM_WORLD," timing: %s - div(M) - took %g s\n",__FUNCT__,t_t3-t_t2);
  } else {
    PetscInfo2(0," timing: %s - div(M) - took %g s\n",__FUNCT__,t_t3-t_t2);
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

#ifndef U1DIRICHLET
  /* shift potential u1 to set value on first node to zero
     this fixes problem with u1 "running away"
     TODO: better use ||u1||/n_vert ???
  */
  PetscReal u1a,*ta_u1;
  ierr = VecGetArray(u1,&ta_u1);CHKERRQ(ierr);
  u1a = -ta_u1[0];
  ierr = VecRestoreArray(u1,&ta_u1);CHKERRQ(ierr);
  if (fabs(u1a)>1e6) {
    ierr = VecShift(u1,u1a);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"u1 shifted by %g\n",u1a);
  }
#endif

  /* do we need to call KSPSetInitialGuess* before every KSPSolve?
     KSP tutorial ex 9 indicates that
  */
  ierr = KSPSetInitialGuessNonzero(ksp_Au1,PETSC_TRUE);CHKERRQ(ierr);

/* TODO: better performance with this!?
  ierr = KSPSetInitialGuessKnoll(ksp_Au1,PETSC_TRUE);CHKERRQ(ierr);
*/

  ierr = KSPSolve(ksp_Au1,bu1,u1);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp_Au1,(PetscInt*)&its1);CHKERRQ(ierr);

  KSPConvergedReason reason;
  ierr = KSPGetConvergedReason(ksp_Au1,&reason);CHKERRQ(ierr);

  PetscInfo2(0,"Au1*u1=divM: %i its, reason: %i\n",its1, reason);
  if (reason<0) {
    PetscReal rnorm;

/*  does not work:

    PCGetFactoredMatrix
    Gets the factored matrix from the preconditioner context.
    This routine is valid only for the LU, incomplete LU, Cholesky, and incomplete Cholesky methods.

    PC pc;
    Mat M;
    Vec D;
    ierr = KSPGetPC(ksp_Au1,&pc);CHKERRQ(ierr);
    ierr = PCGetFactoredMatrix(pc,&M);CHKERRQ(ierr);
    ierr = VecDuplicate(bu1,&D);CHKERRQ(ierr);
    ierr = MatGetDiagonal(M,D);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,"Pivots:\n");
    ierr = VecView(D,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
*/

    ierr = KSPGetResidualNorm(ksp_Au1,&rnorm);CHKERRQ(ierr);
    SETERRQ3(PETSC_ERR_CONV_FAILED,
      "Warning: Au1*u1=divM: KSP did not converge at t=%g ns - reason: %i, rnorm: %g\n",
      gdata->time*gdata->tscale*1e9,reason,rnorm
    );
  }

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

  if (oprofile) {
    PetscPrintf(PETSC_COMM_WORLD," timing: %s - u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
  } else {
    PetscInfo2(0," timing: %s - u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
  }

/* skip the calculation of u2 and the rest, if u1 has not changed!
   with CG and ICC its1 is always >=1
   but LU solver checks residual first and may give its==0

   PETSc 2.3.0:
   great speed up for cg-icc
   does not work with lu
*/
if (its1>0) {
  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                          calculate u2=B*u1
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /* copy boundary values of u1 into u1bnd */
  ierr = MatMult(ABu1,u1,u1bnd);CHKERRQ(ierr);

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

  if (oprofile) {
    PetscPrintf(PETSC_COMM_WORLD," timing: %s - u1 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
  } else {
    PetscInfo2(0," timing: %s - u1 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
  }

  /* calculate u2=B*u1 (boundary vertices) */
  ierr = MatMult(Bmat,u1bnd,u2bnd);CHKERRQ(ierr);

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

  if (oprofile) {
    PetscPrintf(PETSC_COMM_WORLD," timing: %s - u2=B*u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
  } else {
    PetscInfo2(0," timing: %s - u2=B*u1 - took %g s\n",__FUNCT__,t_t3-t_t2);
  }

  /* copy boundary values of u2 (in u2bnd) into global load vector bu2 */
  ierr = MatMult(ABu2,u2bnd,bu2);CHKERRQ(ierr);

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

  if (oprofile) {
    PetscPrintf(PETSC_COMM_WORLD," timing: %s - u2 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
  } else {
    PetscInfo2(0," timing: %s - u2 scatter - took %g s\n",__FUNCT__,t_t3-t_t2);
  }

#ifdef HDEMAGUEPOL
    /* extrapolate u2 (weakened) */
/*
  if (gdata->mode==0 && PetscAbsReal(dt) > 1e-8){
    PetscReal dt2;
    dt2=dt*0.5;
    ierr = VecWAXPY(u2,dt2,du2dt,u2bak);CHKERRQ(ierr);
  }
*/
#endif

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = KSPSetInitialGuessNonzero(ksp_Au2,PETSC_TRUE);CHKERRQ(ierr);

/* TODO: better performance with this!?
  ierr = KSPSetInitialGuessKnoll(ksp_Au2,PETSC_TRUE);CHKERRQ(ierr);
*/

  ierr = KSPSolve(ksp_Au2,bu2,u2);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp_Au2,(PetscInt*)&its2);CHKERRQ(ierr);
  ierr = KSPGetConvergedReason(ksp_Au2,&reason);CHKERRQ(ierr);

  PetscInfo2(0,"Au2*u2=bu2: %i its, reason: %i\n",its2, reason);
  if (reason<0) {
    PetscReal rnorm;

    ierr = KSPGetResidualNorm(ksp_Au2,&rnorm);CHKERRQ(ierr);
    PetscPrintf(PETSC_COMM_WORLD,
      "Warning: Au2*u2=bu2: KSP did not converge at t=%g ns - reason: %i, rnorm: %g\n",
      gdata->time*gdata->tscale*1e9,reason,rnorm
    );
  }

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

  if (oprofile) {
    PetscPrintf(PETSC_COMM_WORLD," timing: %s - u2 - took %g s\n",__FUNCT__,t_t3-t_t2);
  } else {
    PetscInfo2(0," timing: %s - u2 - took %g s\n",__FUNCT__,t_t3-t_t2);
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      calculate H = -grad u
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = VecWAXPY(utot,1.0,u1,u2);CHKERRQ(ierr);
  ierr = MatMult(Agradu,utot,gdata->VHdem);CHKERRQ(ierr);

#ifdef HDEMAGUEPOL
  if (gdata->mode==0 && PetscAbsReal(dt) > 1e-8){
    dt=1.0/dt;

    /* calculate rate of change of u1 */
    ierr = VecCopy(u1,du1dt);CHKERRQ(ierr);
    ierr = VecAXPY(du1dt,c_mone,u1bak);CHKERRQ(ierr);
    ierr = VecScale(du1dt,dt);CHKERRQ(ierr);
    ierr = VecCopy(u1,u1bak);CHKERRQ(ierr);

    /* calculate rate of change of u2 */
/*
    ierr = VecCopy(u2,du2dt);CHKERRQ(ierr);
    ierr = VecAXPY(du2dt,c_mone,u2bak);CHKERRQ(ierr);
    ierr = VecScale(du2dt,dt);CHKERRQ(ierr);
    ierr = VecCopy(u2,u2bak);CHKERRQ(ierr);
*/

    tubak=gdata->time;
  }
#endif

/* endif (its1>0) */
}
else {
  PetscInfo2(0,
    "Info (t=%g ns): Skipped calculation of u2: %i\n",
    gdata->time*gdata->tscale*1e9,
    its1
  );
}

  ierr = VecAXPY(VHtotsum,1.0,gdata->VHdem);CHKERRQ(ierr);

/* DEBUG
  PetscPrintf(PETSC_COMM_WORLD,"deb02: bu1\n");
  ierr = VecView(bu1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb02: u1\n");
  ierr = VecView(u1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb02: u1bnd\n");
  ierr = VecView(u1bnd,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb02: u2bnd\n");
  ierr = VecView(u2bnd,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb02: bu2\n");
  ierr = VecView(bu2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb02: u2\n");
  ierr = VecView(u2,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb02: utot\n");
  ierr = VecView(utot,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"deb02: VHdem\n");
  ierr = VecView(gdata->VHdem,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
*/

  MagparFunctionProfReturn(0);
}


int Hdemag_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
{
  PetscReal e;

  MagparFunctionInfoBegin;

  if (!fieldon) {
    *energy=0.0;
    MagparFunctionInfoReturn(0);
  }

  ierr = VecDot(VMom,gdata->VHdem,&e);CHKERRQ(ierr);
  e /= -2.0;

  *energy=e;

  MagparFunctionProfReturn(0);
}