hcubic.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 static char const Id[] = "$Id: hcubic.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00026 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00027
00028 #include "field.h"
00029
00030 static int doinit=1;
00031 static int fieldon=0;
00032 static Vec VHcubic;
00033 #ifdef EXCH
00034 static PetscReal Ecubic;
00035 #endif
00036
00037 int Hcubic_Init(GridData *gdata)
00038 {
00039 MagparFunctionLogBegin;
00040
00041 for (int i=0;i<gdata->n_prop;i++) {
00042 if (gdata->propdat[NP*i+16]==1 ||
00043 PetscAbsReal(gdata->propdat[NP*i+ 3]) > D_EPS) {
00044
00045
00046
00047 fieldon=1;
00048 break;
00049 }
00050 }
00051
00052 if (fieldon) {
00053 PetscPrintf(PETSC_COMM_WORLD,"Hcubic on\n");
00054 ierr = VecDuplicate(gdata->M,&VHcubic);CHKERRQ(ierr);
00055 }
00056 else {
00057 PetscPrintf(PETSC_COMM_WORLD,"Hcubic off\n");
00058 }
00059
00060 MagparFunctionLogReturn(0);
00061 }
00062
00063
00064 int Hcubic(GridData *gdata,Vec VHtotsum)
00065 {
00066 MagparFunctionInfoBegin;
00067
00068
00069 if (doinit) {
00070 ierr = Hcubic_Init(gdata);CHKERRQ(ierr);
00071 doinit=0;
00072 }
00073 if (!fieldon) {
00074 MagparFunctionInfoReturn(0);
00075 }
00076
00077 PetscReal *ta_Hcubic,*ta_M;
00078 ierr = VecGetArray(VHcubic,&ta_Hcubic);CHKERRQ(ierr);
00079 ierr = VecGetArray(gdata->M,&ta_M);CHKERRQ(ierr);
00080
00081 #ifdef EXCH
00082 PetscReal *ta_VMs3;
00083 ierr = VecGetArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00084 Ecubic=0.0;
00085 #endif
00086
00087
00088 for (int j=0; j<gdata->ln_vert; j++) {
00089 if (gdata->propdat[NP*gdata->vertprop[j]+16]) {
00090
00091
00092 PetscReal k1,k2,jsi;
00093 jsi=gdata->propdat[NP*gdata->vertprop[j]+4];
00094
00095 jsi /= (jsi*jsi+1.e-20);
00096 k1=gdata->propdat[NP*gdata->vertprop[j]+2]*jsi;
00097 k2=gdata->propdat[NP*gdata->vertprop[j]+3]*jsi;
00098
00099
00100 PetscReal ax,ay,az;
00101 ax=gdata->propdat[NP*gdata->vertprop[j]+6];
00102 ay=gdata->propdat[NP*gdata->vertprop[j]+7];
00103 az=gdata->propdat[NP*gdata->vertprop[j]+8];
00104
00105
00106 PetscReal bx,by,bz;
00107 bx=gdata->propdat[NP*gdata->vertprop[j]+10];
00108 by=gdata->propdat[NP*gdata->vertprop[j]+11];
00109 bz=gdata->propdat[NP*gdata->vertprop[j]+12];
00110
00111
00112 PetscReal cx,cy,cz;
00113 cx=gdata->propdat[NP*gdata->vertprop[j]+13];
00114 cy=gdata->propdat[NP*gdata->vertprop[j]+14];
00115 cz=gdata->propdat[NP*gdata->vertprop[j]+15];
00116
00117
00118 PetscReal mx,my,mz;
00119 mx=ta_M[ND*j+0];
00120 my=ta_M[ND*j+1];
00121 mz=ta_M[ND*j+2];
00122
00123
00124 PetscReal a1,a2,a3;
00125 a1=mx*ax+my*ay+mz*az;
00126 a2=mx*bx+my*by+mz*bz;
00127 a3=mx*cx+my*cy+mz*cz;
00128
00129 PetscReal a12,a22,a32;
00130 a12=a1*a1;
00131 a22=a2*a2;
00132 a32=a3*a3;
00133
00134
00135 PetscReal b1,b2,b3;
00136 b1=-2.0*a1*(k1*(a22+a32)+k2*a22*a32);
00137 b2=-2.0*a2*(k1*(a12+a32)+k2*a12*a32);
00138 b3=-2.0*a3*(k1*(a12+a22)+k2*a12*a22);
00139
00140
00141 ta_Hcubic[ND*j+0] = ax*b1+bx*b2+cx*b3;
00142 ta_Hcubic[ND*j+1] = ay*b1+by*b2+cy*b3;
00143 ta_Hcubic[ND*j+2] = az*b1+bz*b2+cz*b3;
00144
00145 #ifdef EXCH
00146 Ecubic +=
00147 (k1*(a12*a22+a12*a32+a22*a32)+k2*a12*a22*a32)*ta_VMs3[ND*j+0];
00148 }
00149 else if (gdata->propdat[NP*gdata->vertprop[j]+3] != 0.0) {
00150 PetscReal k2,jsi;
00151
00152
00153 jsi=gdata->propdat[NP*gdata->vertprop[j]+4];
00154 jsi /= (jsi*jsi+1.e-20);
00155
00156
00157 k2=gdata->propdat[NP*gdata->vertprop[j]+3]*jsi;
00158
00159
00160 PetscReal ax,ay,az;
00161 ax=gdata->propdat[NP*gdata->vertprop[j]+6];
00162 ay=gdata->propdat[NP*gdata->vertprop[j]+7];
00163 az=gdata->propdat[NP*gdata->vertprop[j]+8];
00164
00165 PetscReal a1,a2;
00166 a1=ax*ta_M[ND*j+0]+ay*ta_M[ND*j+1]+az*ta_M[ND*j+2];
00167
00168 a2=4.0*k2*a1*(1.0-a1*a1);
00169
00170
00171 Ecubic += k2*(1.0-a1*a1)*(1.0-a1*a1)*ta_VMs3[ND*j+0];
00172
00173
00174 ta_Hcubic[ND*j+0] = ax*a2;
00175 ta_Hcubic[ND*j+1] = ay*a2;
00176 ta_Hcubic[ND*j+2] = az*a2;
00177 #endif
00178 }
00179 else {
00180 ta_Hcubic[ND*j+0] = 0.0;
00181 ta_Hcubic[ND*j+1] = 0.0;
00182 ta_Hcubic[ND*j+2] = 0.0;
00183 }
00184 }
00185
00186 ierr = VecRestoreArray(gdata->M,&ta_M);CHKERRQ(ierr);
00187 ierr = VecRestoreArray(VHcubic,&ta_Hcubic);CHKERRQ(ierr);
00188 #ifdef EXCH
00189 ierr = VecRestoreArray(gdata->VMs3,&ta_VMs3);CHKERRQ(ierr);
00190 #endif
00191
00192 ierr = VecAXPY(VHtotsum,1.0,VHcubic);CHKERRQ(ierr);
00193
00194 MagparFunctionProfReturn(0);
00195 }
00196
00197
00198 int Hcubic_Energy(GridData *gdata,Vec VMom,PetscReal *energy)
00199 {
00200 PetscReal e;
00201
00202 MagparFunctionInfoBegin;
00203
00204 if (!fieldon) {
00205 *energy=0.0;
00206 MagparFunctionInfoReturn(0);
00207 }
00208
00209 #ifdef EXCH
00210 ierr = PetscGlobalSum(&Ecubic,&e,PETSC_COMM_WORLD);CHKERRQ(ierr);
00211 #else
00212
00213 ierr = VecDot(VMom,VHcubic,&e);CHKERRQ(ierr);
00214 e /= -4.0;
00215 #endif
00216
00217 *energy=e;
00218
00219 MagparFunctionProfReturn(0);
00220 }
00221