GEMS3K Code for Solving for GeoChemical Equilibria

How to Couple with GEMS3K at the TNodeArray Level

The input files to start calculations with GEMS3K can be obtained from running the GEM-Selektor package as described in Section 2.4 of the gemipm2k-man.html file (or in Technical Report), and can also be modified using any ASCII text editor. The file formats are described in detail in gemipm2k-inp.html file. The names of IPM, DCH and one or several DBR files should be collected in a list (txt file) with the file name/path given as a string in the first parameter of GEM_init() function call. Contents of this list usually look like this:

-t "Test-dch.dat" "Test-ipm.dat" "Test-dbr-0.dat"

These files should be located in the same directory where the list file is located.

Programming example fragment (see full example in /GEMS3K/main.cpp):

Code fragments below illustrate how the TNodeArray functionality can be used in implementing a new coupled RMT code (see also main.cpp file in nodearray-gem/ subdirectory). At the first step, two node arrays are allocated using the TNodeArray constructor. Then the initial distribution index list nodeIn[] is defined, and the GEMS3K input files are read in. The contents of nodeIn tell the program where to copy the contents of just read DBR file. In this way, the GEM_init() function initializes the C0 layer of the node array.

#include "nodearray.h"
const int nNodes = 10;   // set here how many nodes you need
const int nTimes = 20;    // set here how many time loops will be performed

int main( int argc, char* argv[] )
{
    . . . . . . . . . . . . . .
    TNodeArray* na;
    // The NodeArray must be allocated here
    TNodeArray::na = na = new TNodeArray( nNodes );

    // Prepare the index array for initial node system and boundary condition codes
    int nodeIn[nNodes] = { 0, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
    int Mode, RetCode;

    // Here we read the IPM, DCH and DBR files prepared from GEMS_PSI.
    // There must be two DBR files prepared - the first will be loaded into
    //   the node 0 and the second - into nodes 1 until 9.
    if( na->GEM_init( ipm_input_file_list_name, nodeIn ) )
          return 1;   // error reading files

The next step usually consists in getting direct access to node array layers:

    DATABRPTR* C0 = na->pNodT0(); // nodes at previous time point
    DATABRPTR* C1 = na->pNodT1(); // nodes at next time point
    bool* iaN = na->piaNode();     // indicators for AIA (true) or SIA (false) for the nodes

After that, the mass transport parameters can be set for each node. For instance, your own input file can be read here. The data for each node are inserted in the node array as follows:

    for( int in=0; in<nNodes; in++) // Loading mass transport parameters
    {
        C0[in]->NodeTypeHY = normal;
    C0[in]->Vt = 1.;
        C0[in]->vp = 1e-7; // cin >> C0[in]->vp;
        C0[in]->eps = 0;
        C0[in]->Km = 0;
        C0[in]->al = 0;
        C0[in]->Dif = 0;
        C0[in]->hDl = 0; // C0[j]->al*C0[j]->vp+C0[j]->Dif;
        C0[in]->nto = 0;
        // cin >> C0[in]->TC;
        iaN[in] = true;
     }
     C0[0]->NodeTypeHY = NBC3source;
     C0[nNodes-1]->NodeTypeHY = NBC3sink;

C0[in] is here a DATABR pointer to the node with index in. Hence, the identifiers like vp or Dif are just names of DATABR
structure fields (see databr.h file). In this loop, additional modifications to bulk composition or temperature can be added, if necessary. The nodes from 1 to 8 are set as normal nodes. The node with in=0 is set as a Cauchy source and the node with in=9 - as a Cauchy sink boundary condition.

To complete initialization of node array layers, a loop over nodes is performed, where the nodes are copied from level C0 to level C1, and chemical equilibria are initially calculated:

   for ( int in=0; in<nNodes; in++) // Copying node layer 0 to node layer 1
   {                                 // (for MT models that can reset time step dt)
        Mode = NEED_GEM_AIA;
        na->CopyNodeFromTo( in, nNodes, C0, C1 );
        RetCode = na->RunGEM( in, Mode ); // running GEM for the in-th node
        if( !(RetCode==OK_GEM_AIA || RetCode == OK_GEM_SIA ))
           return 5;
        na->CopyNodeFromTo( in, nNodes, C1, C0 );
   }

Now both layers are identical at time 0. The transport simulation time loop can begin.

   . . . . . . . . . . . . . . . . . .
   bool TimeStepAccept; // Flag set to true if the time step dt is accepted
   double ct = 0., dt = 1.;
   int xiSr = na->IC_name_to_xDB("Sr"); // Retrieving DATABR index for independent component Sr
   int xiCl = na->IC_name_to_xDB("Cl"); // same for Cl (assuming that both are present in DCH
                                         // and DBR files)
   if( xiSr < 0 || xiCl < 0 )
return 4;    // error - data for one of these elements is not present in DATABR

   FILE* logfile;    // Opening file for logging results
   logfile = fopen( logging_file_name, "w+" );
   if( !logfile)
       return -1;

   for( int it=0; it<nTimes; it++ ) // iterations over time
   {
     int in; // node index
     ct += dt; // current time increment
     iaN[0] = false;
     // Loop over nodes for calculating the mass transport step
     for( in=1; in<nNodes; in++ )
     {
       iaN[in] = false; // true; // SIA mode of GEMIPM2 is allowed for this node
       // In this place, add operators as function of tc and dt that transfer the matter
       // between the nodes according to transport model of choice.
       // Take care about the mass conservation.
       // in this example, we add SrCl2 to the old bIC vector of the previous node at
       // previous time to modify composition in this node at the next time step
     node1_bIC( in, xiSr ) = node0_bIC( in-1, xiSr ) += dt * C0[in]->vp * 1.;
     node1_bIC( in, xiCl ) = node0_bIC( in-1, xiCl ) += dt * C0[in]->vp * 2.;
        . . . . . . . . . . . . . . . . . . .
       // Uncomment this to test variable pressures and temperatures
       // C1[in]->TC = C0[in-1]->TC + 0.2;
       // C1[in]->P = C0[in-1]->P + 2.;
       // . . . . . . . . . . . . . . . . . .
     }
     // The above loop can be replaced by a function call that performs the MT integration step
     // Here the file output may be inserted
       . . . . . . . . . . . . . . . . . . . .
     // Loop over nodes for calculating the chemical equilibration step
     for( in=0; in<nNodes; in++ )
     {
        if( iaN[in] )
            C1[in]->NodeStatusCH = NEED_GEM_AIA; // Automatic (simplex) initial approximation AIA
        else
            C1[in]->NodeStatusCH = NEED_GEM_SIA; // Smart initial approximation SIA - faster
       // Here we can also check if the GEM calculation for this node is needed at all
       // e.g. by comparing this node in C1 layer with that in C0 layer
       if( in == 0 ) // First node (Cauchy source) is not re-equilibrated since it is not changed
          NeedGEM = false;
       else NeedGEM = true;
       node1_bIC( in, xiZz ) = 0.;   // zeroing charge off in node bulk composition
       if( !NeedGEM )
       {   // GEM calculation for this node not needed
           C1[in]->IterDone = 0; // number of GEMIPM iterations is set to 0
       }
       else {// Calling GEM IPM2 calculation
           Mode = C1[in]->NodeStatusCH; // mode of initial aproximation
         RetCode = na->RunGEM( in, Mode ); // running GEM for the node in (node layer C1)
           if( !(RetCode==OK_GEM_AIA || RetCode == OK_GEM_SIA ))
              return 5;
       }
     } // end of node equilibration loop
     . . . . . . . . . . . . . . . . . .
     // Here back-coupling of phase equilibrium data with transport parameters,
   // as well as various checks for the acceptance of the time step can be performed
     TimeStepAccept = true; // false;
   . . . . . . . . . . . . . . . . . .
     if( TimeStepAccept ) // Time step accepted
        for ( in=1; in<nNodes; in++) // Copying node layer C1 to node layer C0
          na->CopyNodeFromTo( in, nNodes, C1, C0 );
     else {
        if(dt > 1e-3)
        {
            ct-= dt;
            dt /= 2.; // Trying the time loop again with twice smaller time step
        }
        else
            return 6; // Error: Time step dt becomes too small
     }
     // Here the output for the current state at tc can be implemented
     . . . . . . . . . . . . . . . . . . .
     na->logProfileAqIC( logfile, it, ct, nNodes, 5 ); // logging elemental molarity profiles over
                                                       // for every fifth time step
     // na->logProfilePhVol( logfile, it, ct, nNodes, 5 ); // logging phase volumes at time step
     . . . . . . . . . . . . . . . . . . .
   dt = 1; // re-setting time step
   } // End of time loop

When the coupled modeling is finished, the whole cleanup of node array with its dynamic data is done by calling the destructor:

delete na;

As seen from the above code fragments, the TNodeArray class takes care of the most routine operations, thus letting the developer concentrate mainly on the transport and back-coupling model. The node.h and nodearray.h files contain many functions and macroses for accessing most of the chemical data stored for nodes at C0 and C1 levels.

End of File