Skip to content

Commit

Permalink
fix mem init
Browse files Browse the repository at this point in the history
  • Loading branch information
Pascal Frey committed Apr 29, 2021
1 parent 482421c commit 2e14ef5
Show file tree
Hide file tree
Showing 4 changed files with 57 additions and 32 deletions.
28 changes: 24 additions & 4 deletions sources/inout.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,8 +141,8 @@ int loadMesh(NSst *nsst) {

/* load initial solution */
int loadSol(NSst *nsst) {
double bufd[GmfMaxTyp],pmin,pmax,ddp;
float fp,buf[GmfMaxTyp];
double bufd[GmfMaxTyp];
int i,k,dim,ver,np,type,inm,typtab[GmfMaxTyp],offset;
char *ptr,data[128];

Expand Down Expand Up @@ -178,13 +178,17 @@ int loadSol(NSst *nsst) {
if ( nsst->info.verb != '0' ) fprintf(stdout," %s :",data);

/* read mesh solutions u_1,..,u_d,p */
pmin = DBL_MAX;
pmax = -DBL_MAX;
GmfGotoKwd(inm,GmfSolAtVertices);
if ( ver == GmfFloat ) {
for (k=0; k<nsst->info.np; k++) {
GmfGetLin(inm,GmfSolAtVertices,buf);
for (i=0; i<nsst->info.dim; i++)
nsst->sol.u[nsst->info.dim*k+i] = buf[i];
nsst->sol.p[k] = buf[i];
if ( nsst->sol.p[k] < pmin ) pmin = nsst->sol.p[k];
else if ( nsst->sol.p[k] > pmax ) pmax = nsst->sol.p[k];
}
}
else {
Expand All @@ -193,8 +197,15 @@ int loadSol(NSst *nsst) {
for (i=0; i<nsst->info.dim; i++)
nsst->sol.u[nsst->info.dim*k+i] = bufd[i];
nsst->sol.p[k] = bufd[i];
if ( nsst->sol.p[k] < pmin ) pmin = nsst->sol.p[k];
else if ( nsst->sol.p[k] > pmax ) pmax = nsst->sol.p[k];
}
}
/* rescale pressure field */
ddp = fabs(pmax - pmin);
for (k=0; k<nsst->info.np; k++)
nsst->sol.p[k] /= ddp;

if ( GmfStatKwd(inm,GmfIterations) ) {
GmfGotoKwd(inm,GmfIterations);
GmfGetLin(inm,GmfTime,&np);
Expand Down Expand Up @@ -222,7 +233,7 @@ int loadSol(NSst *nsst) {

/* save solution */
int saveSol(NSst *nsst,int it) {
double dbuf[GmfMaxTyp];
double dbuf[GmfMaxTyp],pmin,pmax,ddp;
float fbuf[GmfMaxTyp],tmpf;
int i,k,outm,type,typtab[GmfMaxTyp];
char *ptr,data[128],buf[64];
Expand Down Expand Up @@ -262,21 +273,30 @@ int saveSol(NSst *nsst,int it) {
typtab[0] = GmfVec;
typtab[1] = GmfSca;

/* rescale pressure field */
pmin = DBL_MAX;
pmax = -DBL_MAX;
for (k=0; k<nsst->info.np+nsst->info.np2; k++) {
if ( nsst->sol.p[k] < pmin ) pmin = nsst->sol.p[k];
else if ( nsst->sol.p[k] > pmax ) pmax = nsst->sol.p[k];
}
ddp = fabs(pmax - pmin);

/* write P1 sol */
GmfSetKwd(outm,GmfSolAtVertices,nsst->info.np+nsst->info.np2,type,typtab);
if ( nsst->info.ver == GmfFloat ) {
for (k=0; k<nsst->info.np+nsst->info.np2; k++) {
for (i=0; i<nsst->info.dim; i++)
fbuf[i] = nsst->sol.u[nsst->info.dim*k+i];
fbuf[i] = nsst->sol.p[k];
fbuf[i] = nsst->sol.p[k] / ddp;
GmfSetLin(outm,GmfSolAtVertices,fbuf);
}
}
else {
for (k=0; k<nsst->info.np+nsst->info.np2; k++) {
for (i=0; i<nsst->info.dim; i++)
dbuf[i] = nsst->sol.u[nsst->info.dim*k+i];
dbuf[i] = nsst->sol.p[k];
dbuf[i] = nsst->sol.p[k] / ddp;
GmfSetLin(outm,GmfSolAtVertices,dbuf);
}
}
Expand Down
9 changes: 5 additions & 4 deletions sources/nstokes.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -448,10 +448,6 @@ int main(int argc,char **argv) {
int ier,nel,nh;
char stim[32];

memset(&nsst,0,sizeof(NSst));
tminit(nsst.info.ctim,TIMEMAX);
chrono(ON,&nsst.info.ctim[0]);

/* trap exceptions */
signal(SIGABRT,excfun);
signal(SIGFPE,excfun);
Expand All @@ -462,6 +458,11 @@ int main(int argc,char **argv) {
signal(SIGBUS,excfun);

/* init structure */
memset(&nsst,0,sizeof(NSst));
memset(&nsst.info,0,sizeof(Info));
tminit(nsst.info.ctim,TIMEMAX);
chrono(ON,&nsst.info.ctim[0]);

memset(&nsst.mesh,0,sizeof(Mesh));
memset(&nsst.sol,0,sizeof(Sol));
nsst.sol.cl = (Cl*)calloc(NS_CL,sizeof(Cl));
Expand Down
41 changes: 24 additions & 17 deletions sources/nstokes1_2d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,11 +304,11 @@ static int matPNe_P1(double *a,double *b,double PNe[4][4]) {
static double w[2] ={1./3., 1./6.};

/* lenght(a,b) */
dx = b[0]-a[0];
dy = b[1]-a[1];
dx = b[0] - a[0];
dy = b[1] - a[1];
d = sqrt(dx*dx+dy*dy);
eps = 1.0 / NS_EPS;
/* normal to the edge a,b oriented according to the vector ab */
eps = 1.0 / NS_EPS2;
/* normal: orthogonal to the edge a,b oriented according to the vector ab */
n[0] = -dy / d;
n[1] = dx / d;

Expand Down Expand Up @@ -341,12 +341,13 @@ static int matPFe_P1(double *a,double *b,double PFe[4][4]) {
char i;

/* Distance between a and b */
dx = b[0]-a[0];
dy = b[1]-a[1];
dx = b[0] - a[0];
dy = b[1] - a[1];
d = sqrt(dx*dx+dy*dy);
/* Tangent to the edge a, b oriented according to the vector ab */
t[0] = dx/d;
t[1] = dy/d;
t[0] = dx / d;
t[1] = dy / d;

for (i=0; i<4; i++)
PFe[i][3-i]= d*w[1]*t[0]*t[1];
PFe[0][1]= d*w[0]*t[0]*t[1];
Expand All @@ -355,10 +356,10 @@ static int matPFe_P1(double *a,double *b,double PFe[4][4]) {
for (i=0; i<2; i++)
PFe[i][i]= d*w[0]*t[i]*t[i];
for (i=2; i<4; i++)
PFe[i][i]= d*w[0]*t[i-2]*t[i-2];
PFe[i][i]= d*w[0]*t[i-2]*t[i-2];
PFe[0][2]= PFe[2][0] = d*w[1]*t[0]*t[0];
PFe[1][3]= PFe[3][1] = d*w[1]*t[1]*t[1];

return(1);
}

Expand All @@ -373,7 +374,7 @@ static int slipon_2d(NSst *nsst,pCsr A) {
pEdge pa;
pCl pcl;
double *a,*b,PNe[4][4],PFe[4][4];
int j,k,l,dof,iga,igb,ier;
int j,k,l,dof,iga,igb;

dof = nsst->info.typ == P1 ? 2 : 3;

Expand All @@ -388,29 +389,35 @@ static int slipon_2d(NSst *nsst,pCsr A) {
a = &nsst->mesh.point[pa->v[0]].c[0];
b = &nsst->mesh.point[pa->v[1]].c[0];

ier = (nsst->info.typ == P1) ? matPNe_P1(a,b,PNe) : matPNe_P2(a,b,PNe);
if ( !ier ) continue;
if ( nsst->info.typ == P1 ) {
matPNe_P1(a,b,PNe);
matPFe_P1(a,b,PFe);
}
else {
matPNe_P2(a,b,PNe);
matPFe_P2(a,b,PFe);
}

/* P1b P1 */
for (j=0; j<2; j++) {
for (l=0; l<2; l++) {
csrPut(A,iga+j,iga+l,PNe[j][l]);
csrPut(A,igb+j,igb+l,PNe[dof+j][dof+l]);
csrPut(A,iga+j,igb+l,PNe[j][dof+l]);

csrPut(A,igb+j,igb+l,PNe[dof+j][dof+l]);
csrPut(A,igb+j,iga+l,PNe[dof+j][l]);
}
}

/* global friction term */
if ( fabs(pcl->u[0]) < NS_EPSD ) continue;

ier = (nsst->info.typ == P1) ? matPFe_P1(a,b,PFe) : matPFe_P2(a,b,PFe);
if ( !ier ) continue;
for (j=0; j<2; j++) {
for (l=0; l<2; l++) {
csrPut(A,iga+j,iga+l,pcl->u[0] * PFe[j][l]);
csrPut(A,igb+j,igb+l,pcl->u[0] * PFe[dof+j][dof+l]);
csrPut(A,iga+j,igb+l,pcl->u[0] * PFe[j][dof+l]);

csrPut(A,igb+j,igb+l,pcl->u[0] * PFe[dof+j][dof+l]);
csrPut(A,igb+j,iga+l,pcl->u[0] * PFe[dof+j][l]);
}
}
Expand Down
11 changes: 4 additions & 7 deletions sources/tools.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,22 +40,19 @@ int getMat(pSol sol,int ref,double *nu,double *rho) {
/* solve 3x3 definite symmetric positive */
static int sol3x3(double ma[6],double mb[3],double mx[3]) {
double aa,bb,cc,dd,ee,ff,det;
double vmin,vmax,maxx;
double vmin,vmax;
int k;

/* check ill-conditionned matrix */
vmin = vmax = fabs(ma[0]);
for (k=1; k<6; k++) {
maxx = fabs(ma[k]);
if ( maxx < vmin ) vmin = maxx;
if ( maxx > vmax ) vmax = maxx;
dd = fabs(ma[k]);
if ( dd < vmin ) vmin = dd;
else if ( dd > vmax ) vmax = dd;
}
if ( vmax < NS_EPSD ) return(0);
else if ( vmin > NS_EPS2 && vmin < vmax*NS_EPS2 ) return(0);

/* theory says: epsr = EPS2 * vmax; */
/* epsr = vmin == 0.0 ? EPS2 * vmax : vmax / vmin * EPS2; */

/* compute sub-dets */
aa = ma[3]*ma[5] - ma[4]*ma[4];
bb = ma[4]*ma[2] - ma[1]*ma[5];
Expand Down

0 comments on commit 2e14ef5

Please sign in to comment.