Commit c2e5abdf authored by alainm's avatar alainm
Browse files

changed a few loops but output data seems screwed


git-svn-id: https://svn.oca.eu/codes/DISK/fargOCA/branches/viscousTerms@2900 d5a9fe6a-16c7-4aab-97df-e63fc41ba8b9

Former-commit-id: 447dbd1d70a11dbe06fce58e681942428bb77a9e
parent 0dd8d276
......@@ -282,130 +282,121 @@ ComputeViscousTerms(PolarGrid *RadialVelocity,
real const *invdthetamed = InvDthetaMed->field();
real const onethird = 1.0/3.0;
#pragma omp parallel
{
#pragma omp for nowait
for (int i = 0; i < nr; i++) { /* Drr, Dtt, Dpp and divV computation */
for( int h = 0; h < ni; h++) {
int k = h+i*ni;
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int lip = l+ns*ni;
int ljp = l+1;
int lhp = l+ns;
if (j == ns-1) {
ljp = h*ns+i*ni*ns;
}
Drr[l] = (vr[lip]-vr[l])*InvDiffRsup[i];
Dtt[l] = (vt[ljp]-vt[l])*invdtheta[l]+0.5*(vr[lip]+vr[l])*InvRmed[i]+
(vp[lhp]+vp[l])*0.5*InvRmed[i]*CosPhiMed[h]/SinPhiMed[h];
Dpp[l] = (vp[lhp]-vp[l])*InvDphi[k]+0.5*(vr[lip]+vr[l])*InvRmed[i];
divergence[l] = (vr[lip]*Rsup[i]*Rsup[i]-vr[l]*Rinf[i]*Rinf[i])*InvDiffRsup[i]*InvRmed[i]*InvRmed[i];
divergence[l] += (vt[ljp]-vt[l])*invdtheta[l];
divergence[l] += (vp[lhp]*SinPhi[h+1]-vp[l]*SinPhi[h])*InvDphi[k]/SinPhiMed[h];
}
#pragma omp parallel for
for (int i = 1; i < nr-1; i++) { /* Drr, Dtt, Dpp and divV computation */
for( int h = 0; h < ni; h++) {
int k = h+i*ni;
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int lip = l+ns*ni;
int ljp = l+1;
int lhp = l+ns;
if (j == ns-1) {
ljp = h*ns+i*ni*ns;
}
Drr[l] = (vr[lip]-vr[l])*InvDiffRsup[i];
Dtt[l] = (vt[ljp]-vt[l])*invdtheta[l]+0.5*(vr[lip]+vr[l])*InvRmed[i]+
(vp[lhp]+vp[l])*0.5*InvRmed[i]*CosPhiMed[h]/SinPhiMed[h];
Dpp[l] = (vp[lhp]-vp[l])*InvDphi[k]+0.5*(vr[lip]+vr[l])*InvRmed[i];
divergence[l] = (vr[lip]*Rsup[i]*Rsup[i]-vr[l]*Rinf[i]*Rinf[i])*InvDiffRsup[i]*InvRmed[i]*InvRmed[i];
divergence[l] += (vt[ljp]-vt[l])*invdtheta[l];
divergence[l] += (vp[lhp]*SinPhi[h+1]-vp[l]*SinPhi[h])*InvDphi[k]/SinPhiMed[h];
}
}
#pragma omp for
for (int i = 1; i < nr; i++) { /* Drt computation */
for( int h = 0; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns + i*ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int lim = l-ns*ni;
Drt[l] = 0.5*(Rinf[i]*(vt[l]*InvRmed[i]-vt[lim]*InvRmed[i-1])*InvDiffRmed[i]+
(vr[l]-vr[ljm])*invdthetamed[l]);
}
}
#pragma omp parallel for
for (int i = 1; i < nr-1; i++) { /* Drt computation */
for( int h = 0; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns + i*ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int lim = l-ns*ni;
Drt[l] = 0.5*(Rinf[i]*(vt[l]*InvRmed[i]-vt[lim]*InvRmed[i-1])*InvDiffRmed[i]+
(vr[l]-vr[ljm])*invdthetamed[l]);
}
}
#pragma omp for
for (int i = 1; i < nr; i++) { /* Drp computation */
for(int h = 1; h < ni; h++) {
int k = h+i*ni;
for (int j = 0; j < ns; j++) {
int l = j+h*ns + i*ni*ns;
int lim = l-ns*ni;
int lhm = l-ns;
Drp[l] = 0.5*((vr[l]-vr[lhm])*InvDphiMed[k]+
Rinf[i]*(vp[l]/Rmed[i]-vp[lim]/Rmed[i-1])*InvDiffRmed[i]);
}
}
#pragma omp parallel for
for (int i = 1; i < nr-1; i++) { /* Drp computation */
for(int h = 1; h < ni; h++) {
int k = h+i*ni;
for (int j = 0; j < ns; j++) {
int l = j+h*ns + i*ni*ns;
int lim = l-ns*ni;
int lhm = l-ns;
Drp[l] = 0.5*((vr[l]-vr[lhm])*InvDphiMed[k]+
Rinf[i]*(vp[l]/Rmed[i]-vp[lim]/Rmed[i-1])*InvDiffRmed[i]);
}
}
#pragma omp for
for (int i = 0; i < nr; i++) { /* Dtp computation */
for(int h = 1; h < ni; h++) {
int k = h+i*ni;
for (int j = 0; j < ns; j++) {
int l = j+h*ns + i*ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int lhm = l-ns;
Dtp[l] = 0.5*((vp[l]-vp[ljm])*invdthetamed[l]+
SinPhi[h]*(vt[l]/SinPhiMed[h]-vt[lhm]/SinPhiMed[h-1])*InvDphiMed[k]);
}
}
#pragma omp parallel for
for (int i = 1; i < nr-1; i++) { /* Dtp computation */
for(int h = 1; h < ni; h++) {
int k = h+i*ni;
for (int j = 0; j < ns; j++) {
int l = j+h*ns + i*ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int lhm = l-ns;
Dtp[l] = 0.5*((vp[l]-vp[ljm])*invdthetamed[l]+
SinPhi[h]*(vt[l]/SinPhiMed[h]-vt[lhm]/SinPhiMed[h-1])*InvDphiMed[k]);
}
}
}
#pragma omp parallel
{
#pragma omp for nowait
for(int i = 0; i < nr; i++) { /* TAUrr TAUtt and TAUpp computation */
for(int h = 0; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
Trr[l] = 2.0*rho[l]*visc[l]*(Drr[l]-onethird*divergence[l]);
Ttt[l] = 2.0*rho[l]*visc[l]*(Dtt[l]-onethird*divergence[l]);
Tpp[l] = 2.0*rho[l]*visc[l]*(Dpp[l]-onethird*divergence[l]);
}
#pragma omp parallel for
for(int i = 1; i < nr-1; i++) { /* TAUrr TAUtt and TAUpp computation */
for(int h = 0; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
Trr[l] = 2.0*rho[l]*visc[l]*(Drr[l]-onethird*divergence[l]);
Ttt[l] = 2.0*rho[l]*visc[l]*(Dtt[l]-onethird*divergence[l]);
Tpp[l] = 2.0*rho[l]*visc[l]*(Dpp[l]-onethird*divergence[l]);
}
}
#pragma omp for
for (int i = 1; i < nr; i++) { /* TAUrt computation */
for(int h = 0; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int lim = l-ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int ljmim=ljm-ns*ni;
Trt[l] = 2.0*0.25*(rho[l]+rho[lim]+rho[ljm]+rho[ljmim])*visc[l]*Drt[l];
}
}
#pragma omp parallel for
for (int i = 1; i < nr-1; i++) { /* TAUrt computation */
for(int h = 0; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int lim = l-ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int ljmim=ljm-ns*ni;
Trt[l] = 2.0*0.25*(rho[l]+rho[lim]+rho[ljm]+rho[ljmim])*visc[l]*Drt[l];
}
}
#pragma omp for
for (int i = 1; i < nr; i++) { /* TAUrp computation */
for(int h = 1; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int lim = l-ni*ns;
int lhm = l-ns;
Trp[l] = 2.0*0.25*(rho[l]+rho[lhm]+rho[lim]+rho[lim-ns])*visc[l]*Drp[l];
}
}
#pragma omp parallel for
for (int i = 1; i < nr-1; i++) { /* TAUrp computation */
for(int h = 1; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int lim = l-ni*ns;
int lhm = l-ns;
Trp[l] = 2.0*0.25*(rho[l]+rho[lhm]+rho[lim]+rho[lim-ns])*visc[l]*Drp[l];
}
}
#pragma omp for
for (int i = 0; i < nr; i++) { /* TAUtp computation */
for(int h = 1; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int lhm = l-ns;
Ttp[l] = 2.0*0.25*(rho[l]+rho[lhm]+rho[ljm]+rho[ljm-ns])*visc[l]*Dtp[l];
}
}
#pragma omp parallel for
for (int i = 1; i < nr-1; i++) { /* TAUtp computation */
for(int h = 1; h < ni; h++) {
for (int j = 0; j < ns; j++) {
int l = j+h*ns+i*ni*ns;
int ljm = l-1;
if (j == 0) {
ljm = h*ns+i*ni*ns+ns-1;
}
int lhm = l-ns;
Ttp[l] = 2.0*0.25*(rho[l]+rho[lhm]+rho[ljm]+rho[ljm-ns])*visc[l]*Dtp[l];
}
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment