Skip to content

Quicker stateq #225

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
allegroLeiden merged 2 commits into from
Apr 21, 2017
Merged

Quicker stateq #225

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9cf2f7d
Slightly rearranged stateq.
imckstewart Apr 20, 2017
3a7817e
Moved unchanging calculation out of loop.
imckstewart Apr 20, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion src/lime.h
View file
Open in desktop
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
Open in desktop
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);
}