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 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <math.h>
00025 #include "../defs.h"
00026 #include "../athena.h"
00027 #include "../prototypes.h"
00028 #include "prototypes.h"
00029 #include "particle.h"
00030 #include "../globals.h"
00031
00032 #ifdef PARTICLES
00033
00034
00035
00036
00037
00038 void grid_limit(Grid *pG, Domain *pD);
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 void init_particle(Grid *pG, Domain *pD)
00053 {
00054 int i, N1T, N2T, N3T, integratortype, interp, tsmode;
00055 Grain *GrArray;
00056 long size = 1000, size1 = 1, size2 = 1;
00057
00058
00059 grid_limit(pG, pD);
00060 N1T = iup-ilp+1;
00061 N2T = jup-jlp+1;
00062 N3T = kup-klp+1;
00063
00064
00065 pG->partypes = par_geti_def("particle","partypes",1);
00066
00067 if (pG->partypes < 0)
00068 ath_error("[init_particle]: Particle types must not be negative!\n");
00069
00070
00071 if(par_exist("particle","parnumcell"))
00072 {
00073
00074 size1 = N1T*N2T*N3T*(long)(pG->partypes*par_geti("particle","parnumcell"));
00075 if (size1 < 0)
00076 ath_error("[init_particle]: Particle number must not be negative!\n");
00077 }
00078
00079 if(par_exist("particle","parnumgrid"))
00080 {
00081
00082 size2 = (long)(pG->partypes*par_geti("particle","parnumgrid"));
00083 if (size2 < 0)
00084 ath_error("[init_particle]: Particle number must not be negative!\n");
00085
00086 size2 = (long)(size2/((double)(pG->Nx1*pG->Nx2*pG->Nx3))*N1T*N2T*N3T);
00087 }
00088
00089 size = MAX(size, MAX(size1, size2));
00090 pG->arrsize = (long)(1.2*size);
00091
00092 pG->particle = (Grain*)calloc_1d_array(pG->arrsize, sizeof(Grain));
00093 if (pG->particle == NULL) goto on_error;
00094
00095 pG->parsub = (GrainAux*)calloc_1d_array(pG->arrsize, sizeof(GrainAux));
00096 if (pG->parsub == NULL) goto on_error;
00097
00098
00099 pG->grproperty = (Grain_Property*)calloc_1d_array(pG->partypes,
00100 sizeof(Grain_Property));
00101 if (pG->grproperty == NULL) goto on_error;
00102
00103 grrhoa = (Real*)calloc_1d_array(pG->partypes, sizeof(Real));
00104 if (grrhoa == NULL) goto on_error;
00105
00106 tstop0 = (Real*)calloc_1d_array(pG->partypes, sizeof(Real));
00107 if (tstop0 == NULL) goto on_error;
00108
00109
00110 for (i=0; i<pG->partypes; i++) {
00111 tstop0[i] = 0.0;
00112 grrhoa[i] = 0.0;
00113 }
00114 alamcoeff = 0.0;
00115
00116
00117 for (i=0; i<pG->partypes; i++)
00118 pG->grproperty[i].integrator = par_geti_def("particle","integrator",2);
00119
00120
00121 interp = par_geti_def("particle","interp",2);
00122 if (interp == 1)
00123 {
00124 getweight = getwei_linear;
00125 ncell = 2;
00126 }
00127 else if (interp == 2)
00128 {
00129 getweight = getwei_TSC;
00130 ncell = 3;
00131 }
00132 else if (interp == 3)
00133 {
00134 getweight = getwei_QP;
00135 ncell = 3;
00136 }
00137 else
00138 ath_error("[init_particle]: Invalid interp value (should be 1, 2 or 3)!\n");
00139
00140
00141 tsmode = par_geti("particle","tsmode");
00142 if (tsmode == 1)
00143 get_ts = get_ts_general;
00144 else if (tsmode == 2)
00145 get_ts = get_ts_epstein;
00146 else if (tsmode == 3)
00147 get_ts = get_ts_fixed;
00148 else
00149 ath_error("[init_particle]: tsmode must be 1, 2 or 3!\n");
00150
00151
00152 pG->Coup = (GPCouple***)calloc_3d_array(N3T, N2T, N1T, sizeof(GPCouple));
00153 if (pG->Coup == NULL) goto on_error;
00154
00155 return;
00156
00157 on_error:
00158 ath_error("[init_particle]: Error allocating memory.\n");
00159 }
00160
00161
00162
00163
00164
00165 void particle_destruct(Grid *pG)
00166 {
00167 free_1d_array(pG->particle);
00168 free_1d_array(pG->parsub);
00169
00170 free_1d_array(pG->grproperty);
00171 free_1d_array(grrhoa);
00172
00173
00174 if (pG->Coup != NULL) free_3d_array(pG->Coup);
00175
00176 return;
00177 }
00178
00179
00180
00181
00182
00183 void particle_realloc(Grid *pG, long n)
00184 {
00185 pG->arrsize = MAX((long)(1.2*pG->arrsize), n);
00186
00187
00188 if ((pG->particle = (Grain*)realloc(pG->particle,
00189 pG->arrsize*sizeof(Grain))) == NULL)
00190 {
00191 ath_error("[init_particle]: Error re-allocating memory with array size\
00192 %ld.\n", n);
00193 }
00194
00195
00196 if ((pG->parsub = (GrainAux*)realloc(pG->parsub,
00197 pG->arrsize*sizeof(GrainAux))) == NULL)
00198 {
00199 ath_error("[init_particle]: Error re-allocating memory with array size\
00200 %ld.\n", n);
00201 }
00202
00203 return;
00204 }
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 void grid_limit(Grid *pG, Domain *pD)
00218 {
00219 int m1, m2, m3;
00220 int my_iproc, my_jproc, my_kproc;
00221
00222 if (pG->Nx1 > 1) m1 = 1;
00223 else m1 = 0;
00224
00225 if (pG->Nx2 > 1) m2 = 1;
00226 else m2 = 0;
00227
00228 if (pG->Nx3 > 1) m3 = 1;
00229 else m3 = 0;
00230
00231
00232 ilp = pG->is - m1*nghost;
00233 iup = pG->ie + m1*nghost;
00234
00235 jlp = pG->js - m2*nghost;
00236 jup = pG->je + m2*nghost;
00237
00238 klp = pG->ks - m3*nghost;
00239 kup = pG->ke + m3*nghost;
00240
00241
00242
00243
00244
00245 get_myGridIndex(pD, pG->my_id, &my_iproc, &my_jproc, &my_kproc);
00246
00247 if ((par_geti_def("grid","ibc_x1",4) == 2) && (my_iproc == 0))
00248 x1lpar = pG->x1_0 + (ilp+m1 + pG->idisp)*pG->dx1;
00249 else
00250 x1lpar = pG->x1_0 + (pG->is + pG->idisp)*pG->dx1;
00251
00252 if ((par_geti_def("grid","obc_x1",4) == 2) && (my_iproc == pD->NGrid_x1-1))
00253 x1upar = pG->x1_0 + (iup + pG->idisp)*pG->dx1;
00254 else
00255 x1upar = pG->x1_0 + (pG->ie + 1 + pG->idisp)*pG->dx1;
00256
00257 if ((par_geti_def("grid","ibc_x2",4) == 2) && (my_jproc == 0))
00258 x2lpar = pG->x2_0 + (jlp+m2 + pG->jdisp)*pG->dx2;
00259 else
00260 x2lpar = pG->x2_0 + (pG->js + pG->jdisp)*pG->dx2;
00261
00262 if ((par_geti_def("grid","obc_x2",4) == 2) && (my_jproc == pD->NGrid_x2-1))
00263 x2upar = pG->x2_0 + (jup + pG->jdisp)*pG->dx2;
00264 else
00265 x2upar = pG->x2_0 + (pG->je + 1 + pG->jdisp)*pG->dx2;
00266
00267 if ((par_geti_def("grid","ibc_x3",4) == 2) && (my_kproc == 0))
00268 x3lpar = pG->x3_0 + (klp+m3 + pG->kdisp)*pG->dx3;
00269 else
00270 x3lpar = pG->x3_0 + (pG->ks + pG->kdisp)*pG->dx3;
00271
00272 if ((par_geti_def("grid","obc_x3",4) == 2) && (my_kproc == pD->NGrid_x3-1))
00273 x3upar = pG->x3_0 + (kup + pG->kdisp)*pG->dx3;
00274 else
00275 x3upar = pG->x3_0 + (pG->ke + 1 + pG->kdisp)*pG->dx3;
00276
00277 if (pG->Nx1 == 1) x1upar += MAX(0.1*fabs(pG->x1_0), 1.);
00278 if (pG->Nx2 == 1) x2upar += MAX(0.1*fabs(pG->x2_0), 1.);
00279 if (pG->Nx3 == 1) x3upar += MAX(0.1*fabs(pG->x3_0), 1.);
00280
00281 return;
00282 }
00283
00284 #endif