00001 #include "../copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <stdio.h>
00037 #include <stdlib.h>
00038 #include <math.h>
00039 #include "../defs.h"
00040 #include "../athena.h"
00041 #include "../prototypes.h"
00042 #include "prototypes.h"
00043 #include "particle.h"
00044 #include "../globals.h"
00045
00046 #ifdef PARTICLES
00047
00048
00049
00050
00051
00052
00053 Real expr_dpar (const Grid *pG, const int i, const int j, const int k);
00054 Real expr_M1par(const Grid *pG, const int i, const int j, const int k);
00055 Real expr_M2par(const Grid *pG, const int i, const int j, const int k);
00056 Real expr_M3par(const Grid *pG, const int i, const int j, const int k);
00057 Real expr_V1par(const Grid *pG, const int i, const int j, const int k);
00058 Real expr_V2par(const Grid *pG, const int i, const int j, const int k);
00059 Real expr_V3par(const Grid *pG, const int i, const int j, const int k);
00060
00061
00062
00063
00064
00065
00066
00067
00068 void particle_to_grid(Grid *pG, Domain *pD, PropFun_t par_prop)
00069 {
00070 int i,j,k, is,js,ks, i0,j0,k0, i1,j1,k1, i2,j2,k2;
00071 int n0 = ncell-1;
00072 long p;
00073 Real drho;
00074 Real weight[3][3][3];
00075 Vector cell1;
00076 Grain *gr;
00077
00078
00079 if (pG->Nx1 > 1) cell1.x1 = 1.0/pG->dx1;
00080 else cell1.x1 = 0.0;
00081
00082 if (pG->Nx2 > 1) cell1.x2 = 1.0/pG->dx2;
00083 else cell1.x2 = 0.0;
00084
00085 if (pG->Nx3 > 1) cell1.x3 = 1.0/pG->dx3;
00086 else cell1.x3 = 0.0;
00087
00088
00089 for (k=klp; k<=kup; k++)
00090 for (j=jlp; j<=jup; j++)
00091 for (i=ilp; i<=iup; i++) {
00092 pG->Coup[k][j][i].grid_d = 0.0;
00093 pG->Coup[k][j][i].grid_v1 = 0.0;
00094 pG->Coup[k][j][i].grid_v2 = 0.0;
00095 pG->Coup[k][j][i].grid_v3 = 0.0;
00096 }
00097
00098
00099 for (p=0; p<pG->nparticle; p++) {
00100 gr = &(pG->particle[p]);
00101
00102 if ((*par_prop)(gr, &(pG->parsub[p]))) {
00103
00104 getweight(pG, gr->x1, gr->x2, gr->x3, cell1, weight, &is, &js, &ks);
00105
00106
00107 k1 = MAX(ks, klp); k2 = MIN(ks+n0, kup);
00108 j1 = MAX(js, jlp); j2 = MIN(js+n0, jup);
00109 i1 = MAX(is, ilp); i2 = MIN(is+n0, iup);
00110
00111 for (k=k1; k<=k2; k++) {
00112 k0 = k-k1;
00113 for (j=j1; j<=j2; j++) {
00114 j0 = j-j1;
00115 for (i=i1; i<=i2; i++) {
00116 i0 = i-i1;
00117
00118 #ifdef FEEDBACK
00119 drho = pG->grproperty[gr->property].m;
00120 #else
00121 drho = 1.0;
00122 #endif
00123 pG->Coup[k][j][i].grid_d += weight[k0][j0][i0]*drho;
00124 pG->Coup[k][j][i].grid_v1 += weight[k0][j0][i0]*drho*gr->v1;
00125 pG->Coup[k][j][i].grid_v2 += weight[k0][j0][i0]*drho*gr->v2;
00126 pG->Coup[k][j][i].grid_v3 += weight[k0][j0][i0]*drho*gr->v3;
00127
00128 }
00129 }
00130 }
00131 }
00132 }
00133
00134 return;
00135 }
00136
00137
00138
00139
00140
00141 void dump_particle_binary(Grid *pG, Domain *pD, Output *pOut)
00142 {
00143 int dnum = pOut->num;
00144 FILE *p_binfile;
00145 char *fname;
00146 long p, nout, my_id;
00147 int i,is,js,ks,h,init_id = 0;
00148 short pos;
00149 Vector cell1;
00150 Real weight[3][3][3];
00151 Real dpar,u1,u2,u3,cs;
00152 #ifdef FEEDBACK
00153 Real stiffness;
00154 #endif
00155 Grain *gr;
00156 float fdata[12];
00157
00158
00159 if((fname = ath_fname(NULL,pG->outfilename,num_digit,dnum,pOut->id,"lis"))
00160 == NULL){
00161 ath_error("[dump_particle_binary]: Error constructing filename\n");
00162 return;
00163 }
00164
00165 if((p_binfile = fopen(fname,"wb")) == NULL){
00166 ath_error("[dump_particle_binary]: Unable to open binary dump file\n");
00167 return;
00168 }
00169
00170
00171 particle_to_grid(pG, pD, property_all);
00172
00173
00174 if (pG->Nx1 > 1) cell1.x1 = 1.0/pG->dx1;
00175 else cell1.x1 = 0.0;
00176 if (pG->Nx2 > 1) cell1.x2 = 1.0/pG->dx2;
00177 else cell1.x2 = 0.0;
00178 if (pG->Nx3 > 1) cell1.x3 = 1.0/pG->dx3;
00179 else cell1.x3 = 0.0;
00180
00181
00182 for (p=0; p<pG->nparticle; p++)
00183 {
00184 gr = &(pG->particle[p]);
00185
00186
00187 getweight(pG, gr->x1, gr->x2, gr->x3, cell1, weight, &is, &js, &ks);
00188 #ifdef FEEDBACK
00189 h = getvalues(pG, weight, is, js, ks, &dpar,&u1,&u2,&u3,&cs,&stiffness);
00190 #else
00191 h = getvalues(pG, weight, is, js, ks, &dpar, &u1, &u2, &u3, &cs);
00192 #endif
00193
00194 pG->parsub[p].dpar = dpar;
00195 }
00196
00197
00198 nout = 0;
00199 for (p=0; p<pG->nparticle; p++)
00200 if ((*(pOut->par_prop))(&(pG->particle[p]), &(pG->parsub[p])))
00201 nout += 1;
00202
00203
00204
00205
00206 fdata[0] = (float)(pG->x1_0 + (pG->is + pG->idisp)*pG->dx1);
00207 fdata[1] = (float)(pG->x1_0 + (pG->ie+1 + pG->idisp)*pG->dx1);
00208 fdata[2] = (float)(pG->x2_0 + (pG->js + pG->jdisp)*pG->dx2);
00209 fdata[3] = (float)(pG->x2_0 + (pG->je+1 + pG->jdisp)*pG->dx2);
00210 fdata[4] = (float)(pG->x3_0 + (pG->ks + pG->kdisp)*pG->dx3);
00211 fdata[5] = (float)(pG->x3_0 + (pG->ke+1 + pG->kdisp)*pG->dx3);
00212 fdata[6] = (float)(par_getd("grid","x1min"));
00213 fdata[7] = (float)(par_getd("grid","x1max"));
00214 fdata[8] = (float)(par_getd("grid","x2min"));
00215 fdata[9] = (float)(par_getd("grid","x2max"));
00216 fdata[10] = (float)(par_getd("grid","x3min"));
00217 fdata[11] = (float)(par_getd("grid","x3max"));
00218
00219 fwrite(fdata,sizeof(float),12,p_binfile);
00220
00221
00222 fwrite(&(pG->partypes),sizeof(int),1,p_binfile);
00223
00224 for (i=0; i<pG->partypes; i++)
00225 {
00226 fdata[0] = (float)(pG->grproperty[i].rad);
00227 fwrite(fdata,sizeof(float),1,p_binfile);
00228 }
00229
00230
00231 fdata[0] = (float)pG->time;
00232 fdata[1] = (float)pG->dt;
00233 fwrite(fdata,sizeof(float),2,p_binfile);
00234
00235
00236
00237
00238 fwrite(&nout,sizeof(long),1,p_binfile);
00239
00240
00241 for (p=0; p<pG->nparticle; p++)
00242 {
00243 gr = &(pG->particle[p]);
00244 if ((*(pOut->par_prop))(gr,&(pG->parsub[p]))) {
00245
00246
00247 fdata[0] = (float)(gr->x1);
00248 fdata[1] = (float)(gr->x2);
00249 fdata[2] = (float)(gr->x3);
00250 fdata[3] = (float)(gr->v1);
00251 fdata[4] = (float)(gr->v2);
00252 fdata[5] = (float)(gr->v3);
00253
00254 fdata[6] = (float)(pG->parsub[p].dpar);
00255 my_id = gr->my_id;
00256 #ifdef MPI_PARALLEL
00257 init_id = gr->init_id;
00258 #endif
00259 pos = gr->pos;
00260
00261 fwrite(fdata,sizeof(float),7,p_binfile);
00262 fwrite(&(gr->property),sizeof(int),1,p_binfile);
00263 fwrite(&(my_id),sizeof(long),1,p_binfile);
00264 fwrite(&(init_id),sizeof(int),1,p_binfile);
00265
00266 }
00267 }
00268
00269 fclose(p_binfile);
00270
00271 return;
00272 }
00273
00274
00275
00276
00277
00278
00279 int property_all(const Grain *gr, const GrainAux *grsub)
00280 {
00281 return 1;
00282 }
00283
00284
00285
00286
00287
00288
00289
00290
00291 Real expr_dpar(const Grid *pG, const int i, const int j, const int k) {
00292 return pG->Coup[k][j][i].grid_d;
00293 }
00294
00295
00296 Real expr_M1par(const Grid *pG, const int i, const int j, const int k) {
00297 return pG->Coup[k][j][i].grid_v1;
00298 }
00299
00300
00301
00302 Real expr_M2par(const Grid *pG, const int i, const int j, const int k) {
00303 return pG->Coup[k][j][i].grid_v2;
00304 }
00305
00306
00307 Real expr_M3par(const Grid *pG, const int i, const int j, const int k) {
00308 return pG->Coup[k][j][i].grid_v3;
00309 }
00310
00311
00312 Real expr_V1par(const Grid *pG, const int i, const int j, const int k) {
00313 if (pG->Coup[k][j][i].grid_d>0.0)
00314 return pG->Coup[k][j][i].grid_v1/pG->Coup[k][j][i].grid_d;
00315 else return 0.0;
00316 }
00317
00318
00319 Real expr_V2par(const Grid *pG, const int i, const int j, const int k) {
00320 if (pG->Coup[k][j][i].grid_d>0.0)
00321 return pG->Coup[k][j][i].grid_v2/pG->Coup[k][j][i].grid_d;
00322 else return 0.0;
00323 }
00324
00325
00326 Real expr_V3par(const Grid *pG, const int i, const int j, const int k) {
00327 if (pG->Coup[k][j][i].grid_d>0.0)
00328 return pG->Coup[k][j][i].grid_v3/pG->Coup[k][j][i].grid_d;
00329 else return 0.0;
00330 }
00331
00332 #endif