MantisBT - ParaView
View Issue Details
IDProjectCategoryView StatusDate SubmittedLast Update
0011960ParaViewFeaturepublic2011-03-10 16:032012-02-18 00:02
ReporterGreg Weirs 
Assigned ToRobert Maynard 
PriorityhighSeverityminorReproducibilityN/A
StatusclosedResolutionfixed 
PlatformOSOS Version
Product VersionDevelopment 
Target VersionFixed in Version3.14 
ProjectSandia
Topic Name12602-cth-derived-density
Typeincorrect functionality
Summary0011960: CTH reader: Optionally calculate "derived variables" density and material densities.
DescriptionIn spcth files, derived variables are not saved in the file because they can (usually) be calculated from the other variables that are written to the file. However, the spcth reader does not do this. The feature request is to optionally compute several derived variables if the variables they depend on are in the file and if the user requests them (through the variable checkboxes on the preperties panel.) Formulas given in Additional Information.

The motivation for this request is that the Prism plugin is much more effective when the material densities are available, so this request is focused on just these quantities. A more complete solution would be to optionally compute all the derived variables, which would be straightforward and not too much more work after computing the material densities.
Steps To Reproduce
Additional InformationHow to compute the density DENS and the material densities DENSM:
All the variables under discussion are cell variables. In spcth files, material variables end in M and the specific material is appended as +n, where n is the material number. So, e.g. there is a material pressure for each material and these pressures are named PM+1, PM+2, ..., PM+N, for n=1, ..., N.

DENS depends on the material masses M (M+1, M+2, ..., M+N), and the cell volume VOL.
First calculate the cell volume VOL. This computation depends on the coordinate system (distinguished by "igm" in the header section of the file) as well as the coordinates:

if (igm==11) Cvol=PI*(x[i+1]*x[i+1]-x[i]*x[i]); # 1D cyl
else if (igm==12) Cvol = 4*PI/3*(x[i+1]*x[i+1]*x[i+1]-x[i]*x[i]*x[i]); # 1D sph
else if (igm==20) Cvol = (y[j+1]-y[j])*(x[i+1]-x[i]); # 2D cart
else if (igm==21) Cvol = PI*(y[j+1]-y[j])*(x[i+1]*x[i+1]-x[i]*x[i]); # 2D cyl
else Cvol = (z[k+1]-z[k])*(y[j+1]-y[j])*(x[i+1]-x[i]); # 3D cart

Then, density:
DENS = sum(M+1, M+2, ...)/VOL

and material densities:
for n = 1,N:
  DENSM+n = M+n/( VOL*VOLM+n) # "+" is not an operator
end
TagsNo tags attached.
Relationships
Attached Filespdf spyvars.pdf (91,423) 2011-03-30 18:32
https://www.vtk.org/Bug/file/8790/spyvars.pdf
Issue History
Date ModifiedUsernameFieldChange
2011-03-10 16:03Greg WeirsNew Issue
2011-03-17 11:22Robert MaynardAssigned To => Robert Maynard
2011-03-17 11:22Robert MaynardStatusbacklog => tabled
2011-03-29 15:00Utkarsh AyachitTarget Version => 3.10.1
2011-03-30 18:32Greg WeirsFile Added: spyvars.pdf
2011-03-30 18:34Greg WeirsNote Added: 0025978
2011-03-31 11:24Utkarsh AyachitTarget Version3.10.1 => 3.12
2011-06-13 09:17Robert MaynardNote Added: 0026846
2011-06-13 09:17Robert MaynardStatustabled => @80@
2011-06-13 09:17Robert MaynardFixed in Version => 3.12
2011-06-13 09:17Robert MaynardResolutionopen => fixed
2011-06-16 13:10Zack GalbreathCategoryFeature Request => Feature
2011-07-22 14:16Utkarsh AyachitProject => Sandia
2011-11-01 22:13Alan ScottNote Added: 0027646
2011-11-01 22:13Alan ScottStatuscustomer review => closed
2011-11-03 13:51Alan ScottNote Added: 0027656
2011-11-03 13:51Alan ScottStatusclosed => backlog
2011-11-03 13:51Alan ScottResolutionfixed => reopened
2011-11-03 13:53Alan ScottNote Added: 0027657
2011-11-03 13:53Alan ScottStatusbacklog => todo
2011-12-06 15:02Robert MaynardTopic Name => 12602-cth-derived-density
2011-12-06 15:02Robert MaynardType => incorrect functionality
2011-12-06 15:04Robert MaynardStatustodo => gatekeeper review
2011-12-09 14:20Utkarsh AyachitStatusgatekeeper review => customer review
2011-12-09 14:20Utkarsh AyachitNote Added: 0027790
2011-12-09 14:21Utkarsh AyachitFixed in Version3.12 => git-master
2012-02-08 17:21Utkarsh AyachitFixed in Versiongit-master => 3.14
2012-02-13 16:28Greg WeirsNote Added: 0028261
2012-02-18 00:02Alan ScottNote Added: 0028301
2012-02-18 00:02Alan ScottStatuscustomer review => closed
2012-02-18 00:02Alan ScottResolutionreopened => fixed

Notes
(0025978)
Greg Weirs   
2011-03-30 18:34   
This is how spy calculates derived variables. spy is the viz system for CTH. I also attached the pages from the spy manual with the variable names and meanings - handy for labels, etc.


    /* velocity magnitude */
    case STM_VMAG:
       vx = blk->CField[STM_VX][k][j][i];
       if (ndim==1) {
         scratch[i] = fabs(vx);
       } else if (ndim==2) {
         vy = blk->CField[STM_VY][k][j][i];
         scratch[i]= sqrt(vx*vx + vy*vy);
       } else {
         vy = blk->CField[STM_VY][k][j][i];
         vz = blk->CField[STM_VZ][k][j][i];
         scratch[i]= sqrt(vx*vx + vy*vy + vz*vz);
       }
       break;

    /* specific kinetic energy, per unit mass - skip the mass multiplier */
    case STM_KE:
       vx = blk->CField[STM_VX][k][j][i];
       if (ndim==1) {
         scratch[i] = 0.5*(vx*vx);
       } else if (ndim==2) {
         vy = blk->CField[STM_VY][k][j][i];
         scratch[i]= 0.5*(vx*vx + vy*vy);
       } else {
         vy = blk->CField[STM_VY][k][j][i];
         vz = blk->CField[STM_VZ][k][j][i];
         scratch[i]= 0.5*(vx*vx + vy*vy + vz*vz);
       }
       break;
 
    /* density */
    case STM_DENS:

       /* total cell mass */
       Mass=0;
       for (m=0; m<nmat; m++) Mass+= blk->MField[1][m][k][j][i];

       /* cell volume */

       if (igm==11) Cvol=PI*(x[i+1]*x[i+1]-x[i]*x[i]);
       else if (igm==12) Cvol = 4*PI/3*(x[i+1]*x[i+1]*x[i+1]-x[i]*x[i]*x[i]);
       else if (igm==20) Cvol = (y[j+1]-y[j])*(x[i+1]-x[i]);
       else if (igm==21) Cvol = PI*(y[j+1]-y[j])*(x[i+1]*x[i+1]-x[i]*x[i]);
       else Cvol = (z[k+1]-z[k])*(y[j+1]-y[j])*(x[i+1]-x[i]);

       scratch[i] = Mass/Cvol;
/*
       if (blk->CField[0][k][j][i]<1) {
        scratch[i] = Mass/Cvol/(1-blk->CField[0][k][j][i]);
       } else {
        scratch[i]=0;
       }
*/

       break;

     /* Material volume fraction */
     case STM_MVOL:
       if (blk->CField[0][0][0]!=NULL) {
        scratch[i] = 1 - blk->CField[0][k][j][i];
       } else {
        for (m=0; m<nmat; m++) scratch[i] += blk->MField[0][m][k][j][i];
       }
       break;

     /* IE - total material specific energy (mat energy per unit mass * mat_mass)/total_mass */
     case STM_IE:

       /* total cell mass */
       Mass=0;
       for (m=0; m<nmat; m++) {
         Mass+= blk->MField[1][m][k][j][i];
         scratch[i]+=blk->MField[1][m][k][j][i]*blk->MField[4][m][k][j][i];
       }
       if (Mass>0) scratch[i]/=Mass;
       break;

     /* XX stress */
     case STM_XXSTRESS:
       scratch[i] = blk->CField[STM_P][k][j][i] - blk->CField[STM_XXDEV][k][j][i];
       break;

     /* YY stress */
     case STM_YYSTRESS:
       scratch[i] = blk->CField[STM_P][k][j][i] - blk->CField[STM_YYDEV][k][j][i];
       break;

     /* ZZ stress */
     case STM_ZZSTRESS:
       scratch[i] = blk->CField[STM_P][k][j][i] + blk->CField[STM_XXDEV][k][j][i] + blk->CField[STM_YYDEV][k][j][i];
       break;

     /* ZZ deviator */
     case STM_ZZDEV:
       scratch[i] = -blk->CField[STM_XXDEV][k][j][i] - blk->CField[STM_YYDEV][k][j][i];
       break;

     /* J2P (now, for 2-D only) */
     case STM_J2P:
       if(igm < 20) {
         if(igm == 11) {
           scratch[i] = sqrt(3*(blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_XXDEV][k][j][i]
                                + blk->CField[STM_YYDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]
                                + blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]));
         } else {
           scratch[i] = 3*fabs(blk->CField[STM_XXDEV][k][j][i])/2;
         }
       } else if (igm < 30) {
           scratch[i] = sqrt(3*(blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_XXDEV][k][j][i]
                                + blk->CField[STM_YYDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]
                                + blk->CField[STM_XYDEV][k][j][i]*blk->CField[STM_XYDEV][k][j][i]
                                + blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]));
       } else {
           scratch[i] = sqrt(3*(blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_XXDEV][k][j][i]
                                + blk->CField[STM_YYDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]
                                + blk->CField[STM_XYDEV][k][j][i]*blk->CField[STM_XYDEV][k][j][i]
                                + blk->CField[STM_YZDEV][k][j][i]*blk->CField[STM_YZDEV][k][j][i]
                                + blk->CField[STM_XZDEV][k][j][i]*blk->CField[STM_XZDEV][k][j][i]
                                + blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]));
       }
       break;

     /* J2PP (now, for 2-D only) */
     case STM_J2PP:
       if (blk->CField[STM_P][k][j][i]>0) {
        if(igm < 20) {
         if(igm == 11) {
           scratch[i] = sqrt(3*(blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_XXDEV][k][j][i]
                                + blk->CField[STM_YYDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]
                                + blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]))/
                                  blk->CField[STM_P][k][j][i];
         } else {
           scratch[i] = 3*fabs(blk->CField[STM_XXDEV][k][j][i])/(2*blk->CField[STM_P][k][j][i]);
         }
        } else if (igm < 30) {
           scratch[i] = sqrt(3*(blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_XXDEV][k][j][i]
                                + blk->CField[STM_YYDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]
                                + blk->CField[STM_XYDEV][k][j][i]*blk->CField[STM_XYDEV][k][j][i]
                                + blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]))/
                                  blk->CField[STM_P][k][j][i];
        } else {
           scratch[i] = sqrt(3*(blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_XXDEV][k][j][i]
                                + blk->CField[STM_YYDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]
                                + blk->CField[STM_XYDEV][k][j][i]*blk->CField[STM_XYDEV][k][j][i]
                                + blk->CField[STM_YZDEV][k][j][i]*blk->CField[STM_YZDEV][k][j][i]
                                + blk->CField[STM_XZDEV][k][j][i]*blk->CField[STM_XZDEV][k][j][i]
                                + blk->CField[STM_XXDEV][k][j][i]*blk->CField[STM_YYDEV][k][j][i]))/
                                  blk->CField[STM_P][k][j][i];
        }
       } else {
       scratch[i]=-1;
       }
       break;

     /* CVMAG */
     case STM_CVMAG:
       ip1=min(i+1,blk->Nx-1);
       vx = 0.5*(blk->CField[STM_VX][k][j][i]+blk->CField[STM_VX][k][j][ip1]);
       if (ndim==1) {
         scratch[i] = fabs(vx);
       } else if (ndim==2) {
         vy = 0.5*(blk->CField[STM_VY][k][j][i]+blk->CField[STM_VY][k][jp1][i]);
         scratch[i]= sqrt(vx*vx + vy*vy);
       } else {
         vy = 0.5*(blk->CField[STM_VY][k][j][i]+blk->CField[STM_VY][k][jp1][i]);
         vz = 0.5*(blk->CField[STM_VZ][k][j][i]+blk->CField[STM_VZ][kp1][j][i]);
         scratch[i]= sqrt(vx*vx + vy*vy + vz*vz);
       }
       break;

     /* CVX, cell centered velocity */
     case STM_CVX:
       ip1=min(i+1,blk->Nx-1);
       vx = 0.5*(blk->CField[STM_VX][k][j][i]+blk->CField[STM_VX][k][j][ip1]);
       scratch[i] = vx;
       break;

     /* CVY, cell centered velocity*/
     case STM_CVY:
       if (ndim>=2) {
         vy = 0.5*(blk->CField[STM_VY][k][j][i]+blk->CField[STM_VY][k][jp1][i]);
         scratch[i]= vy;
       } else {
         scratch[i]= 0.;
       }
       break;

     /* CVZ, cell centered velocity */
     case STM_CVZ:
       if (ndim==2) {
         vz = 0.5*(blk->CField[STM_VZ][k][j][i]+blk->CField[STM_VZ][kp1][j][i]);
         scratch[i]= vz;
       } else {
         scratch[i]= 0.;
       }
       break;

     /* DIVV */
     case STM_DIVV:
       ip1=min(i+1,blk->Nx-1);
       vx = blk->CField[STM_VX][k][j][i];
       vxp = blk->CField[STM_VX][k][j][ip1];
       if (ndim==1) {
         if(igm == 10) {
           scratch[i] = (vxp-vx)/(x[i+1]-x[i]);
         } else if (igm == 11) {
           scratch[i] = 2.0*(vxp*x[i+1]-vx*x[i])/((x[i+1]+x[i])*(x[i+1]-x[i]));
         } else if (igm == 12) {
           scratch[i] = 4.0*(vxp*x[i+1]*x[i+1]-vx*x[i]*x[i])/((x[i+1]+x[i])*(x[i+1]+x[i])*(x[i+1]-x[i]));
         }
       } else if (ndim==2) {
         vy = blk->CField[STM_VY][k][j][i];
         vyp = blk->CField[STM_VY][k][jp1][i];
         if(igm == 20) {
           scratch[i] = (vxp-vx)/(x[i+1]-x[i]) + (vyp-vy)/(y[j+1]-y[j]);
         } else {
           scratch[i] = 2.0*(vxp*x[i+1]-vx*x[i])/((x[i+1]+x[i])*(x[i+1]-x[i])) +
             (vyp-vy)/(y[j+1]-y[j]);
         }
       } else {
         vy = blk->CField[STM_VY][k][j][i];
         vyp = blk->CField[STM_VY][k][jp1][i];
         vz = blk->CField[STM_VZ][k][j][i];
         vzp = blk->CField[STM_VZ][kp1][j][i];
         scratch[i] = (vxp-vx)/(x[i+1]-x[i]) + (vyp-vy)/(y[j+1]-y[j]) + (vzp-vz)/(z[i+1]-z[i]);
       }
       break;


     /* Cell Temperature in K*/
     case STM_TK:
       scratch[i] = blk->CField[STM_T][k][j][i]*EV2K;
       break;

     /* Material Temperature in K*/
     case STM_TKM:
       scratch[i] = blk->MField[3][mat_id][k][j][i]*EV2K;
       break;
 
     /* material density */
     case STM_DENSM:

       /* material mass */
       Mass = blk->MField[1][mat_id][k][j][i];

       /* cell volume */

       if (igm==11) Cvol=PI*(x[i+1]*x[i+1]-x[i]*x[i]);
       else if (igm==12) Cvol = 4*PI/3*(x[i+1]*x[i+1]*x[i+1]-x[i]*x[i]*x[i]);
       else if (igm==20) Cvol = (y[j+1]-y[j])*(x[i+1]-x[i]);
       else if (igm==21) Cvol = PI*(y[j+1]-y[j])*(x[i+1]*x[i+1]-x[i]*x[i]);
       else Cvol = (z[k+1]-z[k])*(y[j+1]-y[j])*(x[i+1]-x[i]);

       if(Mass>0.) {
         scratch[i] = Mass/(Cvol*blk->MField[0][mat_id][k][j][i]);
       } else {
         scratch[i] = 0.;
       }
       break;

     /* angular momenta, 3D only! */
     case STM_LX:
       yc=0.5*(y[j]+y[j+1])-center_of_mass[1];
       zc=0.5*(z[k]+z[k+1])-center_of_mass[2];
       vy = blk->CField[STM_VY][k][j][i];
       vyp = blk->CField[STM_VY][k][jp1][i];
       vz = blk->CField[STM_VZ][k][j][i];
       vzp = blk->CField[STM_VZ][kp1][j][i];
       scratch[i] = 0.5*(yc*(vz+vzp)-zc*(vy+vyp));
       break;

     case STM_LY:
       xc=0.5*(x[i]+x[i+1])-center_of_mass[0];
       zc=0.5*(z[k]+z[k+1])-center_of_mass[2];
       vx = blk->CField[STM_VX][k][j][i];
       vxp = blk->CField[STM_VX][k][j][ip1];
       vz = blk->CField[STM_VZ][k][j][i];
       vzp = blk->CField[STM_VZ][kp1][j][i];
       scratch[i] = 0.5*(zc*(vx+vxp)-xc*(vz+vzp));
       break;

     case STM_LZ:
       xc=0.5*(x[i]+x[i+1])-center_of_mass[0];
       yc=0.5*(y[j]+y[j+1])-center_of_mass[1];
       vx = blk->CField[STM_VX][k][j][i];
       vxp = blk->CField[STM_VX][k][j][ip1];
       vy = blk->CField[STM_VY][k][j][i];
       vyp = blk->CField[STM_VY][k][jp1][i];
       scratch[i] = 0.5*(xc*(vy+vyp)-yc*(vx+vxp));
       break;

     /* Cell locations */
     case STM_X:
       scratch[i] = x[i];
       break;

     case STM_Y:
       if (ndim>=2) scratch[i] = y[j]; else scratch[i]=0;
       break;

     case STM_Z:
       if (ndim==3) scratch[i] = z[k]; else scratch[i]=0;
       break;

     case STM_XC:
       scratch[i] = 0.5*(x[i]+x[i+1]);
       break;

     case STM_YC:
       if (ndim>=2) scratch[i] = 0.5*(y[j]+y[j+1]); else scratch[i]=0;
       break;

     case STM_ZC:
       if (ndim==3) scratch[i] = 0.5*(z[k]+z[k+1]); else scratch[i]=0;
       break;
   }
  }
(0026846)
Robert Maynard   
2011-06-13 09:17   
Committed changes to master to support derived density to SpyPlot Reader.

author Robert Maynard <robert.maynard@kitware.com>
Mon, 13 Jun 2011 13:15:43 +0000 (09:15 -0400)
committer ParaView Topic Stage <kwrobot@kitware.com>
Mon, 13 Jun 2011 13:15:43 +0000 (09:15 -0400)
commit 03583a21194a7c45f6d602016bb52a5813130e4f
tree 4cf9a55a5d866927249c192270f1fb043564b1ed tree | snapshot
parent fae1aede7e313388ee4a011521d3a74972026dda commit | diff
parent 25864a5c172351b0a3aa8c3941dff1005e7ed41b commit | diff
Merge topic '11960-derived-var-support'

25864a5 ENH: SpyPlot Reader now has user option to compute derived vars.
8abbec7 ENH: The derived variable algorithm doesnt require all materials to be loaded.
0e2d9ee BUG: Properly handle down converted material volume fractions.
905aa2c ENH: Adds in support for derived density to be per material.
23401a1 ENH: Reworked the Derived Density code to be faster.
e0178d9 ENH: Added derived density variable for blocks that have bad ghost cells.
1c50fe0 ENH: First draft at adding derived density variable to spy plot.
520125c ENH: Giving vtkSpyPlotUniReader more methods to make it easier to read.
08c3f8a LEAK: Fix a memory leak in the spy plot stream reader.
75df86d BUG: Change the istream buffer which was creating debug issues in MSVC.
(0027646)
Alan Scott   
2011-11-01 22:13   
Tested Linux client, remote server, 3.12.0. Assuming calculations are correct - I don't know how to test them.
(0027656)
Alan Scott   
2011-11-03 13:51   
As per discussion at weekly meeting, we are possibly reading material volume fraction incorrectly. This must be fixed.
(0027657)
Alan Scott   
2011-11-03 13:53   
Setting as todo.
(0027790)
Utkarsh Ayachit   
2011-12-09 14:20   
merged in master (wherever applicable)
(0028261)
Greg Weirs   
2012-02-13 16:28   
Customer review: works as expected.
(0028301)
Alan Scott   
2012-02-18 00:02   
Closing as per Greg.