calc_dMdt.c

Go to the documentation of this file.
00001 /*
00002     This file is part of magpar.
00003 
00004     Copyright (C) 2002-2010 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: calc_dMdt.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";
00026 
00027 #include "llg.h"
00028 #include "util/util.h"
00029 #include "field/field.h"
00030 #ifdef SPINTORQ
00031 #include "spintorq/spintorq.h"
00032 #endif
00033 
00034 int calc_dMdt(TS ts, PetscReal dt, Vec invec, Vec outvec, void *ptr)
00035 {
00036   GridData *gdata = (GridData*)ptr;
00037 
00038   MagparFunctionInfoBegin;
00039 
00040   /* DEBUG
00041   assert(gdata->M==invec);
00042   */
00043 
00044   ierr = Htot(gdata);CHKERRQ(ierr);
00045 
00046   PetscReal *ta_M, *ta_Htot, *ta_outvec;
00047   VecGetArray(gdata->M,&ta_M);
00048   VecGetArray(gdata->VHtot,&ta_Htot);
00049   VecGetArray(outvec,&ta_outvec);
00050 
00051 #ifdef SPINTORQ
00052   PetscReal  *ta_Hspintorque;
00053   PetscReal  *ta_Hspinforce;
00054   if (gdata->VHspintorque!=PETSC_NULL)  VecGetArray(gdata->VHspintorque,&ta_Hspintorque);
00055   else if (gdata->VHspinforce!=PETSC_NULL) VecGetArray(gdata->VHspinforce,&ta_Hspinforce);
00056 #endif
00057 
00058   for (int i=0;i<gdata->ln_vert;i++) {
00059     PetscReal *torq, damp[ND];
00060     PetscReal a1,alpha;
00061 
00062     torq=ta_outvec+ND*i;
00063 
00064     /* skip non-magnetic vertices and vertices with alpha >= RIGID_M_ALPHA */
00065     if (gdata->propdat[NP*gdata->vertprop[i]+4]<=0.0 ||
00066         gdata->propdat[NP*gdata->vertprop[i]+9]>=RIGID_M_ALPHA) {
00067       torq[0]=0.0;
00068       torq[1]=0.0;
00069       torq[2]=0.0;
00070       continue;
00071     }
00072 
00073     douter(ND,ta_M+ND*i,ta_Htot+ND*i,torq);
00074     douter(ND,ta_M+ND*i,torq,damp);
00075 
00076     alpha=gdata->propdat[NP*gdata->vertprop[i]+9];
00077     a1=1.0/(1.0+alpha*alpha);
00078 
00079     if (gdata->mode==2 && gdata->time<0) {
00080       /* relax with large damping and no precession */
00081       a1=2.0;
00082       my_dcopy(ND,damp,1,torq,1);
00083       my_dscal(ND,-a1,torq,1);
00084     }
00085     else {
00086       my_dscal(ND,-a1,torq,1);
00087       my_daxpy(ND,-a1*alpha,damp,1,torq,1);
00088     }
00089 
00090 #ifdef SPINTORQ
00091     /* FIXME: move to hspin*.c */
00092     PetscReal spintorq[ND], spindamp[ND];
00093     if (gdata->VHspintorque!=PETSC_NULL) {
00094       douter(ND,ta_M+ND*i,ta_Hspintorque+ND*i,spintorq);
00095       douter(ND,ta_M+ND*i,spintorq,spindamp);
00096       my_daxpy(ND,+a1*alpha,spintorq,1,torq,1);
00097       my_daxpy(ND,-a1,spindamp,1,torq,1);
00098     }
00099     else if (gdata->VHspinforce!=PETSC_NULL) {
00100       PetscReal  nonadiabatic;
00101 
00102       nonadiabatic=gdata->propdat_e[NPe*gdata->vertprop[i]+2];
00103 
00104       douter(ND,ta_M+ND*i,ta_Hspinforce+ND*i,spintorq);
00105       douter(ND,ta_M+ND*i,spintorq,spindamp);
00106       my_daxpy(ND,-a1*nonadiabatic+a1*alpha,spintorq,1,torq,1);
00107       my_daxpy(ND,-a1*alpha*nonadiabatic-a1,spindamp,1,torq,1);
00108     }
00109 #endif
00110   }
00111 
00112   VecRestoreArray(gdata->M,&ta_M);
00113   VecRestoreArray(gdata->VHtot,&ta_Htot);
00114   VecRestoreArray(outvec,&ta_outvec);
00115 
00116 #ifdef SPINTORQ
00117   if (gdata->VHspintorque!=PETSC_NULL)  VecRestoreArray(gdata->VHspintorque,&ta_Hspintorque);
00118   else if (gdata->VHspinforce!=PETSC_NULL) VecRestoreArray(gdata->VHspinforce,&ta_Hspinforce);
00119 #endif
00120 
00121   MagparFunctionProfReturn(0);
00122 }
00123 

magpar - Parallel Finite Element Micromagnetics Package
Copyright (C) 2002-2010 Werner Scholz