myllgjacobian.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2009 Werner Scholz
00005 
00006     www:   http://www.magpar.net/
00007     email: magpar(at)magpar.net
00008 
00009     magpar is free software; you can redistribute it and/or modify
00010     it under the terms of the GNU General Public License as published by
00011     the Free Software Foundation; either version 2 of the License, or
00012     (at your option) any later version.
00013 
00014     magpar is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017     GNU General Public License for more details.
00018 
00019     You should have received a copy of the GNU General Public License
00020     along with magpar; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 static char const Id[] = "$Id: myllgjacobian.c 2681 2009-07-31 04:30:53Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";
00026 
00027 #include "llg.h"
00028 #include "field/field.h"
00029 #include "util/util.h"
00030 
00031 extern Mat Ad2E_dMidMk;
00032 
00033 int myLLGJacobian(TS ts,PetscReal t,Vec global_in,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
00034 {
00035   GridData      *gdata = (GridData*)ctx;
00036   Mat           A_Jij = *AA;
00037 
00038   int           i;
00039 
00040   MagparFunctionInfoBegin;
00041 
00042 /* just init with unity matrix for testing */
00043 /*
00044   for (i=0; i<ND*gdata->ln_vert; i++) {
00045     ierr = MatSetValue(A_Jij,i,i,1.0,INSERT_VALUES);CHKERRQ(ierr);
00046   }
00047 
00048   PetscInfo(0,"A_Jij matrix assembly complete\n");
00049   ierr = MatAssemblyBegin(A_Jij,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00050   ierr = MatAssemblyEnd(A_Jij,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00051 
00052   *str = SAME_NONZERO_PATTERN;
00053 
00054   MagparFunctionProfReturn(0);
00055 */
00056 
00057   PetscReal *ta_M,*ta_Htot;
00058   ierr = VecGetArray(global_in,&ta_M);CHKERRQ(ierr);
00059   ierr = VecGetArray(gdata->VHtot,&ta_Htot);CHKERRQ(ierr);
00060 
00061 #if PETSC_VERSION < 300
00062   ierr = MatSetOption(A_Jij,MAT_ROWS_SORTED);CHKERRQ(ierr);
00063 #endif
00064 /*
00065   ierr = MatSetOption(A_Jij,MAT_NEW_NONZERO_LOCATION_ERR);CHKERRQ(ierr);
00066 */
00067 
00068   PetscInt       nonzeroes;
00069   const PetscInt *nonzcols;
00070 
00071   PetscTruth info;
00072   PetscOptionsHasName(PETSC_NULL,"-info",&info);
00073   if (info) {
00074     PetscPrintf(PETSC_COMM_WORLD,"Recalculating matrix elements of Jacobian...\n");
00075   }
00076 
00077   for (i=0; i<gdata->ln_vert; i++) {
00078     if (info) {
00079       ierr = ProgressBar(i,gdata->ln_vert,10);
00080     }
00081 
00082     ierr = MatGetRow(
00083       Ad2E_dMidMk,
00084       ND*gdata->vertl2g[i],
00085       &nonzeroes,
00086       (const PetscInt**)&nonzcols,
00087       PETSC_NULL
00088     );CHKERRQ(ierr);
00089     ierr = MatRestoreRow(
00090       Ad2E_dMidMk,
00091       ND*gdata->vertl2g[i],
00092       &nonzeroes,
00093       (const PetscInt**)&nonzcols,
00094       PETSC_NULL
00095     );CHKERRQ(ierr);
00096 
00097     /* row=nodes1=i col=nodes2=j */
00098 
00099     PetscReal *xi,*hi;
00100     xi=ta_M+ND*i;
00101     hi=ta_Htot+ND*i;
00102 
00103     int rows[ND];
00104     rows[0]=gdata->vertl2g[i]*ND;
00105     rows[1]=rows[0]+1;
00106     rows[2]=rows[0]+2;
00107 
00108     /* TODO: incorrect for mode 2: large damping, no precession for t<0 */
00109     PetscReal  a1,alpha;
00110     alpha=gdata->propdat[NP*gdata->vertprop[i]+9];
00111 
00112     /* skip calculation of Jacobian if Ms=0 or M is locked */
00113     if (gdata->propdat[NP*gdata->vertprop[i]+4]<=0.0 || alpha>=RIGID_M_ALPHA) continue;
00114 
00115     a1=1.0/(1.0+alpha*alpha);
00116 
00117     int j;
00118     j=-1;
00119     /* k increased by 1: matrix elements might be 0 by chance and therefore not stored! */
00120     for (int k=0; k<nonzeroes; k++) {
00121       /* skip this triplet, if we have already completed it */
00122       if (nonzcols[k]/ND==j) continue;
00123 
00124       j=nonzcols[k]/ND;
00125 
00126 /*    does not work: matrix elements might be 0 by chance and
00127       therefore not stored
00128 
00129       assert(nonzcols[k]%ND==0);
00130 */
00131 
00132       int cols[ND];
00133       cols[0]=ND*j;
00134       cols[1]=ND*j+1;
00135       cols[2]=ND*j+2;
00136 
00137       PetscReal dHdM[ND*ND];
00138       MatGetValues(Ad2E_dMidMk,ND,(PetscInt*)rows,ND,(PetscInt*)cols,dHdM);
00139 
00140       /* row/column index correct ?
00141          dHdM should be symmetric anyway !? no: due to MatDiagonalScale !
00142       */
00143 
00144       PetscReal  *dh0dxi;
00145       PetscReal  *dh1dxi;
00146       PetscReal  *dh2dxi;
00147 
00148       dh0dxi = dHdM+ND*0;
00149       dh1dxi = dHdM+ND*1;
00150       dh2dxi = dHdM+ND*2;
00151 
00152       /* TODO: explain this: why is this correct?
00153          without the else path: faster, but sphere_larmor does not preserve E, Mz!
00154 
00155          Local field hi[]!=0 only on local node (rows[0]==cols[0]).
00156          However, dHidMj!=0 also on neighboring node
00157       */
00158       /* TODO: incorrect for mode 2: large damping, no precession for t<0 */
00159       PetscReal Jij[ND*ND];
00160       if (rows[0]==cols[0]) {
00161         Jij[ND*0+0] = -(a1*(      -xi[2]*dh1dxi[0]+xi[1]*dh2dxi[0]))-a1*alpha*(   xi[1]*hi[1]+  xi[2]*hi[2]-xi[1]*xi[1]*dh0dxi[0]-xi[2]*xi[2]*dh0dxi[0]+xi[0]*xi[1]*dh1dxi[0]+xi[0]*xi[2]*dh2dxi[0]);
00162         Jij[ND*1+0] = -(a1*(-hi[2]+xi[2]*dh0dxi[0]-xi[0]*dh2dxi[0]))-a1*alpha*(   xi[1]*hi[0]-2*xi[0]*hi[1]+xi[0]*xi[1]*dh0dxi[0]-xi[0]*xi[0]*dh1dxi[0]-xi[2]*xi[2]*dh1dxi[0]+xi[1]*xi[2]*dh2dxi[0]);
00163         Jij[ND*2+0] = -(a1*( hi[1]-xi[1]*dh0dxi[0]+xi[0]*dh1dxi[0]))-a1*alpha*(   xi[2]*hi[0]-2*xi[0]*hi[2]+xi[0]*xi[2]*dh0dxi[0]+xi[1]*xi[2]*dh1dxi[0]-xi[0]*xi[0]*dh2dxi[0]-xi[1]*xi[1]*dh2dxi[0]);
00164 
00165         Jij[ND*0+1] = -(a1*( hi[2]-xi[2]*dh1dxi[1]+xi[1]*dh2dxi[1]))-a1*alpha*(-2*xi[1]*hi[0]+  xi[0]*hi[1]-xi[1]*xi[1]*dh0dxi[1]-xi[2]*xi[2]*dh0dxi[1]+xi[0]*xi[1]*dh1dxi[1]+xi[0]*xi[2]*dh2dxi[1]);
00166         Jij[ND*1+1] = -(a1*(       xi[2]*dh0dxi[1]-xi[0]*dh2dxi[1]))-a1*alpha*(   xi[0]*hi[0]+  xi[2]*hi[2]+xi[0]*xi[1]*dh0dxi[1]-xi[0]*xi[0]*dh1dxi[1]-xi[2]*xi[2]*dh1dxi[1]+xi[1]*xi[2]*dh2dxi[1]);
00167         Jij[ND*2+1] = -(a1*(-hi[0]-xi[1]*dh0dxi[1]+xi[0]*dh1dxi[1]))-a1*alpha*(   xi[2]*hi[1]-2*xi[1]*hi[2]+xi[0]*xi[2]*dh0dxi[1]+xi[1]*xi[2]*dh1dxi[1]-xi[0]*xi[0]*dh2dxi[1]-xi[1]*xi[1]*dh2dxi[1]);
00168 
00169         Jij[ND*0+2] = -(a1*(-hi[1]-xi[2]*dh1dxi[2]+xi[1]*dh2dxi[2]))-a1*alpha*(-2*xi[2]*hi[0]+  xi[0]*hi[2]-xi[1]*xi[1]*dh0dxi[2]-xi[2]*xi[2]*dh0dxi[2]+xi[0]*xi[1]*dh1dxi[2]+xi[0]*xi[2]*dh2dxi[2]);
00170         Jij[ND*1+2] = -(a1*( hi[0]+xi[2]*dh0dxi[2]-xi[0]*dh2dxi[2]))-a1*alpha*(-2*xi[2]*hi[1]+  xi[1]*hi[2]+xi[0]*xi[1]*dh0dxi[2]-xi[0]*xi[0]*dh1dxi[2]-xi[2]*xi[2]*dh1dxi[2]+xi[1]*xi[2]*dh2dxi[2]);
00171         Jij[ND*2+2] = -(a1*(      -xi[1]*dh0dxi[2]+xi[0]*dh1dxi[2]))-a1*alpha*(   xi[0]*hi[0]+  xi[1]*hi[1]+xi[0]*xi[2]*dh0dxi[2]+xi[1]*xi[2]*dh1dxi[2]-xi[0]*xi[0]*dh2dxi[2]-xi[1]*xi[1]*dh2dxi[2]);
00172       }
00173       else {
00174         Jij[ND*0+0] = -(a1*(      -xi[2]*dh1dxi[0]+xi[1]*dh2dxi[0]))-a1*alpha*(                            -xi[1]*xi[1]*dh0dxi[0]-xi[2]*xi[2]*dh0dxi[0]+xi[0]*xi[1]*dh1dxi[0]+xi[0]*xi[2]*dh2dxi[0]);
00175         Jij[ND*1+0] = -(a1*(       xi[2]*dh0dxi[0]-xi[0]*dh2dxi[0]))-a1*alpha*(                            +xi[0]*xi[1]*dh0dxi[0]-xi[0]*xi[0]*dh1dxi[0]-xi[2]*xi[2]*dh1dxi[0]+xi[1]*xi[2]*dh2dxi[0]);
00176         Jij[ND*2+0] = -(a1*(      -xi[1]*dh0dxi[0]+xi[0]*dh1dxi[0]))-a1*alpha*(                            +xi[0]*xi[2]*dh0dxi[0]+xi[1]*xi[2]*dh1dxi[0]-xi[0]*xi[0]*dh2dxi[0]-xi[1]*xi[1]*dh2dxi[0]);
00177 
00178         Jij[ND*0+1] = -(a1*(      -xi[2]*dh1dxi[1]+xi[1]*dh2dxi[1]))-a1*alpha*(                            -xi[1]*xi[1]*dh0dxi[1]-xi[2]*xi[2]*dh0dxi[1]+xi[0]*xi[1]*dh1dxi[1]+xi[0]*xi[2]*dh2dxi[1]);
00179         Jij[ND*1+1] = -(a1*(       xi[2]*dh0dxi[1]-xi[0]*dh2dxi[1]))-a1*alpha*(                            +xi[0]*xi[1]*dh0dxi[1]-xi[0]*xi[0]*dh1dxi[1]-xi[2]*xi[2]*dh1dxi[1]+xi[1]*xi[2]*dh2dxi[1]);
00180         Jij[ND*2+1] = -(a1*(      -xi[1]*dh0dxi[1]+xi[0]*dh1dxi[1]))-a1*alpha*(                            +xi[0]*xi[2]*dh0dxi[1]+xi[1]*xi[2]*dh1dxi[1]-xi[0]*xi[0]*dh2dxi[1]-xi[1]*xi[1]*dh2dxi[1]);
00181 
00182         Jij[ND*0+2] = -(a1*(      -xi[2]*dh1dxi[2]+xi[1]*dh2dxi[2]))-a1*alpha*(                            -xi[1]*xi[1]*dh0dxi[2]-xi[2]*xi[2]*dh0dxi[2]+xi[0]*xi[1]*dh1dxi[2]+xi[0]*xi[2]*dh2dxi[2]);
00183         Jij[ND*1+2] = -(a1*(       xi[2]*dh0dxi[2]-xi[0]*dh2dxi[2]))-a1*alpha*(                             xi[0]*xi[1]*dh0dxi[2]-xi[0]*xi[0]*dh1dxi[2]-xi[2]*xi[2]*dh1dxi[2]+xi[1]*xi[2]*dh2dxi[2]);
00184         Jij[ND*2+2] = -(a1*(      -xi[1]*dh0dxi[2]+xi[0]*dh1dxi[2]))-a1*alpha*(                             xi[0]*xi[2]*dh0dxi[2]+xi[1]*xi[2]*dh1dxi[2]-xi[0]*xi[0]*dh2dxi[2]-xi[1]*xi[1]*dh2dxi[2]);
00185       }
00186       ierr = MatSetValues(A_Jij,ND,(PetscInt*)rows,ND,(PetscInt*)cols,Jij,INSERT_VALUES);CHKERRQ(ierr);
00187     }
00188 
00189 #undef MATFLUSH
00190 #ifdef MATFLUSH
00191     /* flush caches to save memory (and avoid swapping) */
00192     if (i % 100 == 99) {
00193       ierr = MatAssemblyBegin(A_Jij,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00194       ierr = MatAssemblyEnd(A_Jij,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00195     }
00196 #endif
00197   }
00198 
00199 #ifdef MATFLUSH
00200   /* continue assembling to stay in sync with other processors
00201      (MatAssembly is collective!)
00202   */
00203   for (;i<gdata->n_vert;i++) {
00204     if (i % 100 == 99) {
00205       ierr = MatAssemblyBegin(A_Jij,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00206       ierr = MatAssemblyEnd(A_Jij,MAT_FLUSH_ASSEMBLY);CHKERRQ(ierr);
00207     }
00208   }
00209 #endif
00210 
00211   if (info) {
00212     ierr = ProgressBar(1,1,10);
00213   }
00214 
00215   PetscInfo(0,"A_Jij matrix assembly complete\n");
00216   ierr = MatAssemblyBegin(A_Jij,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00217   ierr = MatAssemblyEnd(A_Jij,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
00218 
00219 /*
00220   ierr = MatView(A_Jij,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
00221 */
00222 
00223   VecRestoreArray(global_in,&ta_M );
00224   VecRestoreArray(gdata->VHtot,&ta_Htot);
00225 
00226   /* TODO: optimize: */
00227 /*
00228   *str = SAME_NONZERO_PATTERN;
00229 */
00230   *str = DIFFERENT_NONZERO_PATTERN;
00231 
00232   MagparFunctionProfReturn(0);
00233 }
00234
magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2009 Werner Scholz