View Issue Details [ Jump to Notes ] [ Print ]
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 
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.
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.
ProjectSandia
Topic Name12602-cth-derived-density
Typeincorrect functionality
Attached Filespdf file icon spyvars.pdf [^] (91,423 bytes) 2011-03-30 18:32

 Relationships

  Notes
(0025978)
Greg Weirs (developer)
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 (developer)
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 (manager)
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 (manager)
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 (manager)
2011-11-03 13:53

 Setting as todo.
(0027790)
Utkarsh Ayachit (administrator)
2011-12-09 14:20

 merged in master (wherever applicable)
(0028261)
Greg Weirs (developer)
2012-02-13 16:28

 Customer review: works as expected.
(0028301)
Alan Scott (manager)
2012-02-18 00:02

 Closing as per Greg.

 Issue History
Date Modified Username Field Change
2011-03-10 16:03 Greg Weirs New Issue
2011-03-17 11:22 Robert Maynard Assigned To => Robert Maynard
2011-03-17 11:22 Robert Maynard Status backlog => tabled
2011-03-29 15:00 Utkarsh Ayachit Target Version => 3.10.1
2011-03-30 18:32 Greg Weirs File Added: spyvars.pdf
2011-03-30 18:34 Greg Weirs Note Added: 0025978
2011-03-31 11:24 Utkarsh Ayachit Target Version 3.10.1 => 3.12
2011-06-13 09:17 Robert Maynard Note Added: 0026846
2011-06-13 09:17 Robert Maynard Status tabled => @80@
2011-06-13 09:17 Robert Maynard Fixed in Version => 3.12
2011-06-13 09:17 Robert Maynard Resolution open => fixed
2011-06-16 13:10 Zack Galbreath Category Feature Request => Feature
2011-07-22 14:16 Utkarsh Ayachit Project => Sandia
2011-11-01 22:13 Alan Scott Note Added: 0027646
2011-11-01 22:13 Alan Scott Status customer review => closed
2011-11-03 13:51 Alan Scott Note Added: 0027656
2011-11-03 13:51 Alan Scott Status closed => backlog
2011-11-03 13:51 Alan Scott Resolution fixed => reopened
2011-11-03 13:53 Alan Scott Note Added: 0027657
2011-11-03 13:53 Alan Scott Status backlog => todo
2011-12-06 15:02 Robert Maynard Topic Name => 12602-cth-derived-density
2011-12-06 15:02 Robert Maynard Type => incorrect functionality
2011-12-06 15:04 Robert Maynard Status todo => gatekeeper review
2011-12-09 14:20 Utkarsh Ayachit Status gatekeeper review => customer review
2011-12-09 14:20 Utkarsh Ayachit Note Added: 0027790
2011-12-09 14:21 Utkarsh Ayachit Fixed in Version 3.12 => git-master
2012-02-08 17:21 Utkarsh Ayachit Fixed in Version git-master => 3.14
2012-02-13 16:28 Greg Weirs Note Added: 0028261
2012-02-18 00:02 Alan Scott Note Added: 0028301
2012-02-18 00:02 Alan Scott Status customer review => closed
2012-02-18 00:02 Alan Scott Resolution reopened => fixed

Copyright © 2000 - 2018 MantisBT Team