Recent Forum Posts
From categories:
page »
Re: Optimizing C++ Code by SchrobiSchrobi, 10 Oct 2011 10:23

Following link is focused on embedded system, but it's worth looking in to it anyway, as some tricks are useful in our case.

"Embedded software often runs on processors with limited computation power, thus optimizing the code becomes a necessity. In this article we will explore the following optimization techniques for C and C++ code developed for Real-time and Embedded Systems."


Optimizing C++ Code by SchrobiSchrobi, 10 Oct 2011 10:18

Hey Thomas,

have you had a look in to Moeck et al (2009)?

It gives you a slip tendency Ts along the lines of what I am dealing with in CO2 sequestration. The slip tendency is the ratio of the resolved shear to resolved normal stress on the fault surface which has to exceed the frictional sliding resistance (given by the sliding friction coefficient \mus).

So what you want is to compute the traction vector (which is the force acting on a surface element whereas stress acts on the volume element) for given stress tensor field on each single element in CSMP which represents the fault. Given the traction vector, you can resolve the shear and normal component on the element and compare it against the failure criterion as outlined in Pollard and Fletcher (2005). Maybe it is not necessarly to use the traction vector to compute the normal and shear component and use the farfield stress tensor field directly.

Stephan used this approach for calculation on stress dependent permeability changes for a fracture network using Cruikshank et al (1991). However, the stress field is held constant and does not include fluid pressure coupling.

I am not sure whether a stress redistribution calculation is possible without including displacement. Haven't seen any mehtods on this so far, but will have a thorough look into this. A reference as a starting point would be helpful, you have one?



  1. Moeck, Inga, Grzegorz Kwiatek, and Günter Zimmermann. 2009. Slip tendency analysis, fault reactivation potential and induced seismicity in a deep geothermal reservoir. Journal of Structural Geology 31, no. 10 (October): 1174-1182. doi:10.1016/j.jsg.2009.06.012.
  2. Pollard, David D., and Raymond C. Fletcher. 2005. Fundamentals of structural geology. Cambridge University Press.
  3. Cruikshank, Kenneth M., Guozhu Zhao, and Arvid M. Johnson. 1991. Analysis of minor fractures associated with joints and faulted joints. Journal of Structural Geology 13, no. 8: 865-886. doi:10.1016/0191-8141(91)90083-U.

As for the positive divergence: I am plotting in log-scale, hence only the values > 0 are shown. There are values < 0 though. But as I said, for triangular FV I get divergence free fields.

Here is the volume flux for one of Mandefro's Bristol meshes (with 3 orders permeability difference between fractures and matrix)


and the essentially zero divergence field


So, all seems well for triangular meshes and the new FV scheme but not necessarily quadrilateral meshes and the old scheme.

Another thing that is strange for me in your pictures is, why the divergence is always positive and you don't have any negative value. It means that you are destroying mass in every control volume. This can't be correct!!!

I haven't worked with generic node center finite volume transport for a long time (more than 30 years, means practically not at all).
One thing that I want to address here is the divergence of velocity is equal to the total volume sink/source per volume per time, and it is equal to zero when you don't have sink or source in the domain, which I think it hasn't been considered in CSMP.
I have my own patch class which is my given name for the control volume in my new scheme and I have computed the flux balance over a patch and the divergence is in the range of numerical precision (epsilon) for triangular elements when you don't have a crazy permeability variation (like 1e-14 for matrix and 1e-8 for fracture). But for quadrilateral elements I have problem, and it comes from the pressure solution and it seems that the solver has not converged properly. I think if we tighten the convergence criteria for the solver the flux mismatch will vanish for the quadrilateral as well.
I never agree with adding or removing any mass to the CVs, and in my scheme I never do that and it is fine and I don't have any problem and it is conservative. I think there is something fishy with generic node center finite volume transport class. I propose to check it for a single phase case to find out if the problem is because of the mobilities that has been used in the transport algorithm or it comes from somewhere else.

CSMP seems to have problems creating divergence free velocity fields for quadrilateral FE meshes. Has anyone realised this? We solve, after all

\begin{align} \nabla\cdot\mathbf{v}=0 \end{align}


\begin{align} \nabla\cdot\left(\frac{k}{\mu}\nabla p\right)=0 \end{align}

and the resulting divergence, as computed with the FV scheme (Stephan's generic and my old one) must be numerically zero (this is the case for triangular meshes, see below). Now, I only get this for a homogeneous permeability field and quadrilateral finite elements, as shown in this image:


If I include a low-permeability region in the centre such that flow is around it, then the divergence becomes immediately non-zero.


Has somebody noticed this and/or a good explanation? We only realised this problem when Christine's Newton iterator got stuck trying to resolve these local source/sink terms that arise from the divergence.

Divergence free velocity fields by sebgeisebgei, 03 May 2011 20:08

Hi, has anyone investigated possibilities of including very simple fault mechanic "fakes" in CSMP (planar faults, linear elasticity, no explicit representation of displacements but rather applying some constitutive relations that fake single slip events once a failure criterion is matched and compute stress re-distribution after an event), i.e., something on a much, much simpler level than Adriana's work?

Reason is that in geothermal reservoir stimulation models such things have been used and we wondered if we could include them as well.

Thanks Sebastian, I fixed it now. I tried around, changed the placement of some nodes, tried different command line switches and then it suddenly worked. At some point I got errors from the function SplitSingleCornerElements, but they vanished also.
So, valid Triangle output is not automatically valid CSMP-input all the time, but as you said, small changes do the trick.

Gillian is really the expert in using Triangle and has recently done some changes to the Triangle interface and removed some bugs. One thing I noticed from time to time that the geometry is ok but somehow one certain mesh does not work but if I change the mesh size slightly (e.g., from FE area 1.0 to 0.9 or angle constrain from 30 degrees to 29.5 degrees) the next mesh suddenly works.

I have a problem with the TRIANGLE_Interface. It seems that I do something wrong when creating the triangle files.
The error is:

main: Enter name of 'Triangle' input file set: tetraMesh.1

MeshInterface::CheckTriangleOutput: '.ele' file contains no element attributes…

MeshInterface::CheckTriangleOutput 'tetraMesh.1' Input fileset apparently contains:
99013 nodes, 296148 segments=faces, 197136 triangles.

MeshInterface::ReadNodeDataFile: file 'tetraMesh.1.node' read successfully…

MeshInterface::ReadElementDataFile: file 'tetraMesh.1.ele' read successfully…

MeshInterface::ReadNeighborDataFile: file 'tetraMesh.1.neigh' read successfully…
csp5: ../CSP5.0/source_code/interfaces/CSP_TRIANGLE_Interface.cpp:936: void csp::TRIANGLE_Interface::FlagBoundaryNodesAccordingTo(csp::stl_index, csp::int32, const std::vector<unsigned int, std::allocator<unsigned int> >&, std::map<unsigned int, int, std::less<unsigned int>, std::allocator<std::pair<const unsigned int, int> > >&): Assertion ‘bit != bflags.end()’ failed.

Is there a guideline on how to create a proper triangle-mesh?

If you write darcy equation with gravity term for a single phase:


according to the picture above,


so, the velocity due to gravity force is -k/mu*ro*g*cos(alpha). This is the velocity along the line from gravitational force. To calculate the each component of velocity (in x and y directions), you should multiply the unit vector that connects points p1 and p2 by this value. If you add this velocity to the velocity due to pressure gradient then you have your total velocity for your lower dimensional element.

Re: Bug in 1D line elements by shahoshaho, 26 Nov 2010 14:49

Thanks, Shaho. I had another thought about the problem and here is what I come up with for 1D line elements in a 2D space using basic trigonometry. It is quite a bit different from your code and does not solve the issue that, possible, the calculation of the interpolation function derivatives for 1D line elements oriented in 2D space may be wrong. At least if you run my simple test above where you get different pressure gradients than the true solution.

Consider the following situation


Where the red vector is the velocity within a 1D fracture that has the same orientation. The fracture is represented as a 1D line element with starting point $P1$ (at location x1, y1 and hydrostatic pressure $\pi1$) and endpoint $P2$ (at location x2, y2 with hydrostatic pressure $\pi2$). For simplicity assume that permeability and viscosity are one and only a single phase fluid is present, which is at perfect hydrostatic equilibrium such that no fluid moves up or down and $\partial p/\partial z = \rho g$.

Now, because the fracture is oriented at some angle $\alpha$ to the vertical gravity vector $\mathbf{g}$, there will be a pressure gradient in the x and z direction, causing an "apparent velocity" $\mathbf{v}$ with an x and z component. Obviously, the z component $vz$ vanishes if we subtract $\rho g$ but the x component $vx$ remains, causing an apparent flow to the right, which clearly is non physical.

If we know $\alpha$ , we can compute vector $\mathbf{x}$ using $\tan\alpha=\mathbf{x}/\mathbf{g}$. Note that $\mathbf{x}$ is equivalent to the apparent $vx$ that we would like to subtract from $\mathbf{v}$ to correct for all hydrostatic effects, both in z and x direction. Any non-hydrostatic pressure effects will not be affected. We can obtain $\alpha$ from the scalar product of $\mathbf{g}$ and $\mathbf{v}$, where the orientation of $\mathbf{v}$ is given by $P1$ - $P1$. The scalar product allows us to compute the angle between two vectors as

\begin{align} \mathbf{v}\cdot\mathbf{g}=|v|\cdot|g|\cos\alpha \end{align}

Note that the scalar product only considers the length and orientation of the vectors intrinsically and can easily be solved for $\cos\alpha$ and we can compute $\alpha$ from the secant, i.e the inverse of $\cos\alpha$. Once we know $\alpha$ , we obtain $\mathbf{x}$ from $\tan\alpha=\mathbf{x}/\mathbf{g}$. Subtracting this value and $\mathbf{g}$ from the real velocity vector then corrects for gravity, not only in the vertical direction but also in the horizontal one. Some simple calculations on a piece and paper show that this corrects properly for gravity but leaves any other pressure gradients in order.

The same technique can be applied in 3D but we need another correction for the apparent offset in the y-direction.

Does this make sense?

Re: Bug in 1D line elements by sebgeisebgei, 26 Nov 2010 10:41

Shaho's answer per email frm 26 November

I already saw this bug in csmp and I have fixed it in the new csmp. The gravity projection is not correct and if you calculate the divergence of velocity field, it gives you a non-zero value. There is a bug also in the 3D case when you use surface fractures. I haven't fixed that part yet. I thing you should slowly move to the new code, because we are in the process of unit testing all classes and we found some bugs in other classes. This is the part of code that calculates the total velocity. You can change it to the old csmp code and use it.

void  ReservoirSimulator2D::Compute_vt_AtBaryCenterIncludingGravity()
    const uint32        VERTICAL_AXIS(1U);
    const double64      acc_gravity(9.80665);
    vector<double64>    un(2U);
    const Point<2U>     gvec( 0., -1. );
    ScalarVariable      flux;
    VectorVariable<2U>  velo;
    Region<2U,Element>&  gref(reservoir_model_.Region("Model"));
    for ( vector<Element<2U>*>::iterator 
          it=gref.ElementsBegin(); it!=gref.ElementsEnd(); it++ ) 
         // computing the critical relperm parameters
         relperm_model_->Initialize( *(*it) );
         relperm_model_->InitializeForBaryCenter( *(*it) );
         // computing the total velocity: vt = -k (lt grad p - (lw rho_w + lo rho_o) g)
         const double64 mob_t = relperm_model_->TotalMobility();
         velo = 0.;
         FiniteElementTraits<2U,Element>(*(*it)).dN_AtBaryCenter( DERIV_, 1U );
         for ( size_t i=0U; i<(*it)->Nodes(); i++ ) {
              double64 pf = (*it)->N(i)->Read( props_.pf1_key );
              pf         += (*it)->N(i)->Read( props_.hp_key );
              velo(0)  += pf  * -DERIV_(0,i) * mob_t;
              velo(1)  += pf  * -DERIV_(1,i) * mob_t;
         // taking into account the gravitational flow component for calculation of velocity of non-wetting phase
         //               k * (lambda_w * rho_w + lambda_o * rho_o) g
         const double64  gravity_term = (*it)->Read( props_.gt_key ) * -acc_gravity;
         if ( (*it)->FE()->IsSurfaceElement() ) velo(VERTICAL_AXIS) += gravity_term;
         // special case of lower dimensional elements embedded in 2D
         else {
              // LINE ELEMENT
              // ascertaining that the velocity vector remains aligned with the line element
              if ( (*it)->FE()->IsLineElement() ) {
                   assert( (*it)->FE()->Interpolation() == 1U ); // linear interpolation=two nodes
                   // 1. find unit vector in direction of line element
                   Point<2U> evec((*it)->N(1U)->Coordinate() - (*it)->N(0U)->Coordinate()); 
                   // 2. check direction of evec and correct it to upward
                   if ( evec[1U] < 0. ) evec *= -1.;
                   evec /= evec.Length();
                   // 3. scaling element aligned unit vector by flow vector gravitational component
                   evec *= evec[1U] * gravity_term;
                   // 4. add resulting velocity components to velocity vector
                   velo(0U) += evec[0U];
                   velo(1U) += evec[1U];
         // storing the computed velocity
        (*it)->Store( props_.vt_key, velo );
        // backing up the previous volume flux
        (*it)->Read( props_.vf1_key, flux );
        (*it)->Store( props_.vf0_key, flux );
        // volume flux
        flux = velo.Length();   
        (*it)->Store( props_.vf1_key, flux );
    extern bool global_verbose;
    if ( global_verbose ) printRangeOfVariable( reservoir_model_, "velocity" );
} // end Compute_vt_AtBaryCenterIncludingGravity


Re: Bug in 1D line elements by sebgeisebgei, 26 Nov 2010 10:32

There are some bugs in the CSP_Matrix class, most importantly the copy constructor does not resize the number of rows and matrices during the copy process. So, if you attempt to copy a large matrix to a smaller one (or indeed any two matrixes with different sizes) CSMP dies on you without warning. Further to that, I added some more functionality to the CSP_Matrix class such as LU decomposition and back-substitution, obtaining data from individual rows and columns etc. Please update CSP_Matrix.cpp and CSP_Matrix.h from the Edinburgh SVN, particularly the Leoben team such that these bugs do not get lost during the port.

When trying to run some simulations with gravity effects in 2D models with 1D fracture elements I stumbled across the fact that the velocities in the fractures where always extremely (read m/sec) high even if the system was completely hydrostatic, i.e. no fluid should move up or down. I included the correction for 1D fracture velocities under the influence of gravity as programmed in the CSMP reservoir simulator, which made things even worse although velocities were nicely aligned with the fracture (BTW, I don’t quite follow the logic what is done here). I then went on to reduce the problem to a homogeneous matrix with a single oblique fracture, which still gave ridiculously high velocities in the fracture at a perfectly homogeneous and uniform hydrostatic pressure where no fluid must move up or down. Further invistigations showed that there is something odd with the gradient calculations: The one I obtain using the interpolation function derivatives differ, sometimes a factor 10, from the ones I would calculate by hand.

I stripped down the code to the following stand-along test:

  stl_index DIM(2);
  IsoparametricLinearLineElement  elm(DIM,1); 
  Element<csp_float,2> e(&elm);
  Node<csp_float,2> n0, n1;
  e.ConnectTo( 0, n0 );
  e.ConnectTo( 1, n1 );  
  csp_float l = e.Volume();
  csp_float p[2], v[2];
  p[0] = 0.0; p[1] = 1.0;
  v[0] = v[1] = 0.0;
  DenseMatrix<csp_float,DM_MIN> DN, XY;
  for ( stl_index i=0; i<e.IntegrationPoints(); i++ ) {
      e.dN_AtIntegrationPoint( DN, i, 1 );
      for ( stl_index n=0; n<e.Nodes(); n++ )
          for ( stl_index j=0; j<DIM; j++ ) v[j] += p[n] * DN(j,n);  
  cout << "\n\nFE Gradient X: " << v[0] << ", should be " << (p[1]-p[0])/(n1.x()-n0.x());  
  cout << "\n\nFE Gradient Y: " << v[1] << ", should be " << (p[1]-p[0])/(n1.y()-n0.y());  

This is a 1D line element with two points at (0,0) and (1,1). There is a pressure of 0 at node 0 and a pressure of 1 at node 1. The length of the element is correctly calculated as square root of 2 and the pressure gradient in X and Y direction must be 1 (with the length of the vector between point 0 and 1 is square root of 2). Now, this is not what using the interpolation function derivatives yield. They yield gradients of 0.5 in X and Y directions. If I use coordinates of the nodes and the pressures from the actual mesh, I can recomputed the gradients that the “real” CSMP gives me, which differ from what I would think are the right ones. It only works if the line element is oriented horizontally or vertically.

Now, aside from the gravity effects, this could be very worrying because it implies that all fluid velocities in 2D DFM models where the fractures are 1D line elements are potentially not correct.

Am I doing something entirely stupid or is there an issue with the 1D line elements in 2D space?

Cheers, Sebastian

Bug in 1D line elements by sebgeisebgei, 25 Nov 2010 15:48

And some replies in this weeks issue of Nature (especially the last one is interesting as it refers to a new journal for publishing papers plus code)

Computer code: more credit needed
Computer code: incentives needed
Computer code: a model journal

Did you get anywhere with Stephan's/Julian's suggestion to evaluate the shape functions at point XYZ and test if they sum up to 1? Otherwise I have some code that I can dig out. Not sure though if it is fast enough for your frequent applications (but a good search structure in a stl map may help).

Hi Adriana

The idea would be to create a hybrid mesh 3D/2D of the phase diagram of H2O-NaCl. Phase boundaries would be represented by surfaces of, e.g., of triangles, the regions inbetween by volume elements like tetrahedra. I guess it'll be in the order of 10^5-10^6 elements to have sufficient resolution.

I'd need to pseudo-randomly access ANY xyz coordinate (BUT: for a given node in the simulation domain (not the lookup mesh) I'd probably have a good guess from the last time step in what region I should be).

As always, there is a trade-off to be considered: Currently, my 3D-lookup is a cartesian grid with variable resolution in certain regions. From given xyz (which are temperature-pressure-composition) I can jump to the correct position with one very quick operation and linearly interpolate in that rectangual cell.

This works fine, fast and cheap unless I am in a situation where a phase boundary (e.g. liquid/vapor) lies within that cell. Then there is a plethora of possible topologies how exactly that boundary lies in the cell (and whether there are one, two, or three boundaries). Figuring them out is possible but I'd like to avoid going through that pain. Up to now, I extrapolated across the boundary but this is probably not accurate enough for our upcoming strict transport CVFEM scheme (flow gets potentially unstable if fluid properties are even slightly inconsistent with each other).

So, creating an FE-based lookup table looked attractive as it (a) would give a very good geometric representation and (b) has intrisically consistent interpolation routines. The unattractive side would be the cost of finding the right element for interpolation.

One possibility I thought of is creating the mesh such that it'd essentially consist of a series of long prisms which have there long axis strictly parallel to one of the axes (in my case y) and are themselves then subdivided into smaller prism or tretrahedral elements. This would allow to quickly find the right "long" prism first and then loop within it to find the right spot along the y-axis. However, this is not well thought through yet.

page »
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License