Skip to content

Commit

Permalink
Merge pull request #225 from allegroLeiden/quicker_stateq
Browse files Browse the repository at this point in the history 
Quicker stateq
  • Loading branch information
@allegroLeiden
allegroLeiden authored Apr 21, 2017
2 parents ea5327a + 3a7817e commit fdfaf98
Show file tree
Hide file tree
Showing 2 changed files with 104 additions and 86 deletions.
1 change: 0 additions & 1 deletion src/lime.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -337,7 +337,6 @@ void getclosest(double, double, double, long*, long*, double*, double*, double*)
double geterf(const double, const double);
void getjbar(int, molData*, struct grid*, const int, configInfo*, struct blendInfo, int, gridPointData*, double*);
void getMass(configInfo*, struct grid*, const gsl_rng*);
void getmatrix(int, gsl_matrix*, molData*, struct grid*, int, gridPointData*);
void getVelocities(configInfo *, struct grid *);
void getVelocities_pregrid(configInfo *, struct grid *);
void gridPopsInit(configInfo*, molData*, struct grid*);
Expand Down
189 changes: 104 additions & 85 deletions src/stateq.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,16 +12,110 @@
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_errno.h>

struct collPartMatrixType {
double *ctot;
gsl_matrix *colli;
};

/*....................................................................*/
void
getFixedMatrix(molData *md, int ispec, struct grid *gp, int id, struct collPartMatrixType **collPartMat){
int ipart,k,l,ti;

(*collPartMat) = malloc(sizeof(**collPartMat)*md[ispec].npart);

/* Initialize matrix with zeros */
for(ipart=0;ipart<md[ispec].npart;ipart++){
if(md[ispec].nlev<=0){
if(!silent) bail_out("Matrix initialization error in stateq");
exit(1);
}

(*collPartMat)[ipart].colli = gsl_matrix_alloc(md[ispec].nlev+1,md[ispec].nlev+1);
(*collPartMat)[ipart].ctot = malloc(sizeof(double)*md[ispec].nlev);
for(k=0;k<md[ispec].nlev+1;k++){
for(l=0;l<md[ispec].nlev+1;l++){
gsl_matrix_set((*collPartMat)[ipart].colli, k, l, 0.0);
}
}
}

/* Populate matrix with collisional transitions */
for(ipart=0;ipart<md[ispec].npart;ipart++){
struct cpData part = md[ispec].part[ipart];
double *downrates = part.down;
for(ti=0;ti<part.ntrans;ti++){
int coeff_index = ti*part.ntemp + gp[id].mol[ispec].partner[ipart].t_binlow;
double down = downrates[coeff_index]\
+ gp[id].mol[ispec].partner[ipart].interp_coeff*(downrates[coeff_index+1]\
- downrates[coeff_index]);
double up = down*md[ispec].gstat[part.lcu[ti]]/md[ispec].gstat[part.lcl[ti]]\
*exp(-HCKB*(md[ispec].eterm[part.lcu[ti]]-md[ispec].eterm[part.lcl[ti]])/gp[id].t[0]);
gsl_matrix_set((*collPartMat)[ipart].colli, part.lcu[ti], part.lcl[ti], down);
gsl_matrix_set((*collPartMat)[ipart].colli, part.lcl[ti], part.lcu[ti], up);
}

for(k=0;k<md[ispec].nlev;k++){
(*collPartMat)[ipart].ctot[k]=0.0;
for(l=0;l<md[ispec].nlev;l++)
(*collPartMat)[ipart].ctot[k] += gsl_matrix_get((*collPartMat)[ipart].colli,k,l);
}
}
}

/*....................................................................*/
void
getMatrix(int id, gsl_matrix *matrix, molData *md, struct grid *gp, int ispec, gridPointData *mp, struct collPartMatrixType *collPartMat){
int k,l,li,ipart,di;

/* Initialize matrix with zeros */
for(k=0;k<md[ispec].nlev+1;k++){
for(l=0;l<md[ispec].nlev+1;l++){
gsl_matrix_set(matrix, k, l, 0.);
}
}

/* Populate matrix with radiative transitions */
for(li=0;li<md[ispec].nline;li++){
k=md[ispec].lau[li];
l=md[ispec].lal[li];
gsl_matrix_set(matrix, k, k, gsl_matrix_get(matrix, k, k)+md[ispec].beinstu[li]*mp[ispec].jbar[li]+md[ispec].aeinst[li]);
gsl_matrix_set(matrix, l, l, gsl_matrix_get(matrix, l, l)+md[ispec].beinstl[li]*mp[ispec].jbar[li]);
gsl_matrix_set(matrix, k, l, gsl_matrix_get(matrix, k, l)-md[ispec].beinstl[li]*mp[ispec].jbar[li]);
gsl_matrix_set(matrix, l, k, gsl_matrix_get(matrix, l, k)-md[ispec].beinstu[li]*mp[ispec].jbar[li]-md[ispec].aeinst[li]);
}

for(k=0;k<md[ispec].nlev;k++){
for(ipart=0;ipart<md[ispec].npart;ipart++){
di = md[ispec].part[ipart].densityIndex;
if(di>=0)
gsl_matrix_set(matrix,k,k,gsl_matrix_get(matrix,k,k)+gp[id].dens[di]*collPartMat[ipart].ctot[k]);
}
for(l=0;l<md[ispec].nlev;l++){
if(k!=l){
for(ipart=0;ipart<md[ispec].npart;ipart++){
di = md[ispec].part[ipart].densityIndex;
if(di>=0)
gsl_matrix_set(matrix,k,l,gsl_matrix_get(matrix,k,l)-gp[id].dens[di]*gsl_matrix_get(collPartMat[ipart].colli,l,k));
}
}
}
gsl_matrix_set(matrix, md[ispec].nlev, k, 1.);
gsl_matrix_set(matrix, k, md[ispec].nlev, 0.);
}
}

/*....................................................................*/
void
stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
, struct blendInfo blends, int nextMolWithBlend, gridPointData *mp\
, double *halfFirstDs, _Bool *luWarningGiven){

int t,s,iter,status;
int t,s,iter,status,ipart;
double *opop,*oopop,*tempNewPop=NULL;
double diff;
char errStr[80];
struct collPartMatrixType *collPartMat=NULL;

gsl_matrix *matrix = gsl_matrix_alloc(md[ispec].nlev+1, md[ispec].nlev+1);
gsl_matrix *reduc = gsl_matrix_alloc(md[ispec].nlev, md[ispec].nlev);
Expand All @@ -42,10 +136,12 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
diff=1;
iter=0;

getFixedMatrix(md,ispec,gp,id,&collPartMat);

while((diff>TOL && iter<MAXITER) || iter<5){
getjbar(id,md,gp,ispec,par,blends,nextMolWithBlend,mp,halfFirstDs);

getmatrix(id,matrix,md,gp,ispec,mp);
getMatrix(id,matrix,md,gp,ispec,mp,collPartMat);
for(s=0;s<md[ispec].nlev;s++){
for(t=0;t<md[ispec].nlev-1;t++){
gsl_matrix_set(reduc,t,s,gsl_matrix_get(matrix,t,s));
Expand Down Expand Up @@ -94,6 +190,12 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
iter++;
}

for(ipart=0;ipart<md[ispec].npart;ipart++){
gsl_matrix_free(collPartMat[ipart].colli);
free(collPartMat[ipart].ctot);
}
free(collPartMat);

gsl_matrix_free(matrix);
gsl_matrix_free(reduc);
gsl_vector_free(rhVec);
Expand All @@ -104,88 +206,5 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
free(oopop);
}

void
getmatrix(int id, gsl_matrix *matrix, molData *md, struct grid *gp, int ispec, gridPointData *mp){
int ti,k,l,li,ipart,di;
struct getmatrix {
double *ctot;
gsl_matrix * colli;
} *partner;

partner = malloc(sizeof(struct getmatrix)*md[ispec].npart);

/* Initialize matrix with zeros */
for(ipart=0;ipart<md[ispec].npart;ipart++){
partner[ipart].colli = gsl_matrix_alloc(md[ispec].nlev+1,md[ispec].nlev+1);
if(md[ispec].nlev>0) partner[ipart].ctot = malloc(sizeof(double)*md[ispec].nlev);
else {
if(!silent)bail_out("Matrix initialization error in stateq");
exit(0);
}
for(k=0;k<md[ispec].nlev+1;k++){
for(l=0;l<md[ispec].nlev+1;l++){
gsl_matrix_set(matrix, k, l, 0.);
gsl_matrix_set(partner[ipart].colli, k, l, 0.);
}
}
}

/* Populate matrix with radiative transitions */
for(li=0;li<md[ispec].nline;li++){
k=md[ispec].lau[li];
l=md[ispec].lal[li];
gsl_matrix_set(matrix, k, k, gsl_matrix_get(matrix, k, k)+md[ispec].beinstu[li]*mp[ispec].jbar[li]+md[ispec].aeinst[li]);
gsl_matrix_set(matrix, l, l, gsl_matrix_get(matrix, l, l)+md[ispec].beinstl[li]*mp[ispec].jbar[li]);
gsl_matrix_set(matrix, k, l, gsl_matrix_get(matrix, k, l)-md[ispec].beinstl[li]*mp[ispec].jbar[li]);
gsl_matrix_set(matrix, l, k, gsl_matrix_get(matrix, l, k)-md[ispec].beinstu[li]*mp[ispec].jbar[li]-md[ispec].aeinst[li]);
}

/* Populate matrix with collisional transitions */
for(ipart=0;ipart<md[ispec].npart;ipart++){
struct cpData part = md[ispec].part[ipart];
double *downrates = part.down;
for(ti=0;ti<part.ntrans;ti++){
int coeff_index = ti*part.ntemp + gp[id].mol[ispec].partner[ipart].t_binlow;
double down = downrates[coeff_index]\
+ gp[id].mol[ispec].partner[ipart].interp_coeff*(downrates[coeff_index+1]\
- downrates[coeff_index]);
double up = down*md[ispec].gstat[part.lcu[ti]]/md[ispec].gstat[part.lcl[ti]]\
*exp(-HCKB*(md[ispec].eterm[part.lcu[ti]]-md[ispec].eterm[part.lcl[ti]])/gp[id].t[0]);
gsl_matrix_set(partner[ipart].colli, part.lcu[ti], part.lcl[ti], down);
gsl_matrix_set(partner[ipart].colli, part.lcl[ti], part.lcu[ti], up);
}

for(k=0;k<md[ispec].nlev;k++){
partner[ipart].ctot[k]=0.;
for(l=0;l<md[ispec].nlev;l++)
partner[ipart].ctot[k] += gsl_matrix_get(partner[ipart].colli,k,l);
}
}

for(k=0;k<md[ispec].nlev;k++){
for(ipart=0;ipart<md[ispec].npart;ipart++){
di = md[ispec].part[ipart].densityIndex;
if(di>=0)
gsl_matrix_set(matrix,k,k,gsl_matrix_get(matrix,k,k)+gp[id].dens[di]*partner[ipart].ctot[k]);
}
for(l=0;l<md[ispec].nlev;l++){
if(k!=l){
for(ipart=0;ipart<md[ispec].npart;ipart++){
di = md[ispec].part[ipart].densityIndex;
if(di>=0)
gsl_matrix_set(matrix,k,l,gsl_matrix_get(matrix,k,l)-gp[id].dens[di]*gsl_matrix_get(partner[ipart].colli,l,k));
}
}
}
gsl_matrix_set(matrix, md[ispec].nlev, k, 1.);
gsl_matrix_set(matrix, k, md[ispec].nlev, 0.);
}

for(ipart=0;ipart<md[ispec].npart;ipart++){
gsl_matrix_free(partner[ipart].colli);
free(partner[ipart].ctot);
}
free(partner);
}


0 comments on commit fdfaf98

Please sign in to comment.