GEMS3K Solver of GeoChemical Equilibria
GEMS3K Solver of GeoChemical Equilibria
2. Data exchange using disk files only
3. Data exchange in computer memory
The main.cpp file in "standalone/gemcalc/" folder shows how to perform a batch calculation of chemical equilibria using the input data given in DCH, IPM and one or more DBR files; the results are written into DBR files or (optionally) into a dump text output file.
Main.cpp and main.h files in "standalone/node-gem/" folder shows how the coupling with a "transport part" having several nodes can be arranged, where the chemical data and parameters are kept in dynamic arrays and exchanged with the GEMS3K in computer memory. The modification to a “real” mass transport method is straightforward, although this would make the example program too complicated.
In these examples, some options are shown that may significantly improve performance and accuracy of GEMS3K.As the mass transport part has no direct access to internal IPM work data, the calculations of equilibria occur as if they were done in “real” amounts consistent with the current size of the node. The use of internal re-scaling mechanism leads, in general, to a much better convergence of GEMS3K and to higher precision of its results.
The current chemical speciation (xDC vector, amounts of Dependent Components) is usually the output from GEM IPM runs (can also be used as input initial approximation in SIA mode). The ultimate input of GEMS3K is the bulk chemical composition vector (bIC), temperature and pressure. So, at time steps, the mass transport part should change the bIC vector and, maybe, T and/or P in every node, and this must be sufficient.
However, in many existing advection/diffusion transport codes (e.g. OpenGeoSys), not the total amounts of independent components, but rather the amounts of aqueous species are moved between nodes. Conversely, after the mass transport step, only the speciation vector xDC is changed in each node, and an additional step would be needed to compute a corrected bIC vector (i.e. elemental bulk composition) in the mass transport part, which will be the GEM input for computing the new speciation xDC. In turn, this would require access to several data objects (such as the A stoichiometry matrix) in the DATACH structure.
In order to make this unnecessary, the TNode class provides a variant of GEM_from_MT() function call where the re-calculation of input (changed) xDC speciation into new bIC bulk composition vector is performed inside of GEMS3K before the internal mass scaling and calculation of new equilibrium state. Actually, total amounts of Independent Components are calculated from stoichiometries of Dependent Components (as provided in the A matrix) and their amounts (in the xDC vector passed through the p_xDC parameter), and then added to the amounts given in the bIC vector (passed through p_bIC parameter). Here, two cases are possible. One case is when the mass transport step puts in the p_xDC parameter the differences between the new and the old speciation in the node chemical system. In this case, the old p_bIC content (e.g. unchanged since the last GEM calculation for this node) must be provided and it will be automatically corrected. Another (more usual) case is when the mass-transport part delivers in p_xDC the complete new speciation (absolute amounts) after the transport time step. In this case, the p_bIC vector must be filled with zeros because its new content will be completely retrieved from the p_xDC vector.
Concerning the “Smart Initial Approximation” (SIA) mode, the TNode class contains a parameter in GEM_run() call and yet another variant of GEM_from_MT() function that control what initial approximation of GEM speciation vector will be used. Basically, regardless of what the mass transport part does, the GEMS3K can use previous contents of speciation and activity coefficients vectors in MULTI structure as an initial approximation for the next GEM calculation – even if the bulk composition b or temperature T have had changed. This mode is set by the default uPrimalSol parameter value (false).
Another, more efficient variant (uPrimalSol = true) can be used when the mass transport part directly changes the chemical speciation in the node while taking care of material balance, and the bIC input vector must be restored anyway from the new xDC vector contents (provided together with activity coefficients stored for this node at previous time step). In this case, the SIA mode is very efficient because the mass balance is already ensured, and this may lead to a very serious reduction of the number of GEM iterations and total computing time. Note that this SIA mode works only when the DATACH structure contains a full list of Dependent Components used in GEM calculations.
The TNode class provides dozens of access methods to single items of chemical data. Usage of these public methods is strongly preferred over direct access to fields of the DATABR and DATACH memory structures. The provided methods use either DATABR or DATACH indexation.
In DBR files and DATABR structure, the results for volumes, masses and bulk compositions are only provided for multi-component phases, because for pure phases they can be easily retrieved. Nevertheless, TNode public functions such as Ph_Volume(), Ph_Mass(), Ph_BC() and Ph_SatInd() work both for multi- and single-component phases. This is done in order to make access to such data from the side of mass-transport part of the program as simple as possible.Programming example fragment (see full example in "standalone/gemcalc/main.cpp"):
| #include "node.h" int main( int argc, char* argv[] ) { long nRecipes = 0; char (*recipes)[fileNameLength] = 0; // Analyzing command line arguments // Default arguments char input_system_file_list_name[256] = "system-dat.lst"; char input_recipes_file_list_name[256] = "more_recipes.lst"; if (argc >= 2 ) strncpy( input_system_file_list_name, argv[1], 256); // list of DCH, IPM and DBR input files for initializing GEMS3K // Creates TNode structure instance accessible trough the "node" pointer TNode* node = new TNode(); // (1) Initialization of GEMS3K internal data by reading files // whose names are given in the input_system_file_list_name if( node->GEM_init( input_system_file_list_name ) ) { // error occured during reading the files return 1; } // Getting direct access to work node DATABR structure which exchanges the // data with GEM IPM3 (already filled out by reading the DBR input file) DATABR* dBR = node->pCNode(); // . . . . . . . . . . . . . . . // Asking GEM to run with automatic initial approximation dBR->NodeStatusCH = NEED_GEM_AIA; // (2) re-calculating equilibrium by calling GEMS3K, getting the status back int NodeStatusCH = node->GEM_run( false ); if( NodeStatusCH == OK_GEM_AIA || NodeStatusCH == OK_GEM_SIA ) { // (3) Writing results in default DBR file node->GEM_write_dbr( NULL, false, true, false ); // node->GEM_print_ipm( NULL ); // possible debugging printout } else { // (4) possible return status analysis, error message node->GEM_print_ipm( NULL ); // possible debugging printout return 5; // GEM IPM did not converge properly - error message needed } // . . . . . . . . . . . . . . . // Here a possible loop on more input recipes begins if (argc >= 3 ) { char NextRecipeFileName[256]; char NextRecipeOutFileName[300]; char input_recipes_file_list_path[256-fileNameLength] = ""; strncpy( input_recipes_file_list_name, argv[2], 256); // list of additional recipes (dbr.dat files) e.g. for simulation // of a titration or another irreversible process // Reading list of recipes names from file recipes = f_getfiles( input_recipes_file_list_name, input_recipes_file_list_path, nRecipes, ','); for(int cRecipe=0; cRecipe < nRecipes; cRecipe++ ) { // Trying to read the next file name sprintf(NextRecipeFileName , "%s%s", input_recipes_file_list_path, recipes[cRecipe] ); // (5) Reading the next DBR file with different input composition or temperature node->GEM_read_dbr( NextRecipeFileName ); // Asking GEM IPM2 to run (faster) with smart initial approximation dBR->NodeStatusCH = NEED_GEM_SIA; NodeStatusCH = node->GEM_run( false ); if( NodeStatusCH == OK_GEM_AIA || NodeStatusCH == OK_GEM_SIA ) { sprintf(NextRecipeOutFileName , "%s.out", NextRecipeFileName ); node->GEM_write_dbr( NextRecipeOutFileName, false, false, false ); sprintf(NextRecipeOutFileName , "%s.Dump.out", NextRecipeFileName ); node->GEM_print_ipm( NextRecipeOutFileName ); } else { // error message, debugging printout sprintf(NextRecipeOutFileName , "%s.Dump.out", NextRecipeFileName ); node->GEM_print_ipm( NextRecipeOutFileName ); // return 5; // GEM IPM did not converge properly - error message needed } } } // end of possible loop on input recipes delete node; if( recipes ) delete recipes; // End of example return 0; } |
This case is
typical for coupling GEMS3K with an existing reactive
mass transport code that has its own data structure for keeping
chemical information (e.g. speciation, bulk composition of the system)
for each node. As in the previous example, it uses one work instance of
the DATABR memory
structure as data bridge for the exchange with internal GEM IPM data
structures.
This example code consists of two C++ files: "main.h" and "main.cpp",
located in the "standalone/node-gem/" folder. The former file
contains a definition of TMyTransport
class representing a typical dynamic memory structure used for
keeping current chemical composition and speciation for all nodes
involved in the reactive mass transport system. The latter file
contains
an implementation of some public methods of TMyTransport class, as
well as the main function which shows the sequence of calls to TNode
class functions, necessary to perform coupled reactive transport
calculations based on the "operator splitting" approach.
This example
is designed to help learning how to couple the existing transport code
with GEMS3K in a simple way, and it does not contain a 'real' transport
code which would make the example more difficult to observe. The
transport model itself is located in the TMyTransport::OneTimeStepRun()
function and can be improved there to any desired extent.
Here is a code snippet with the TMyTransport class definition (see complete code in standalone/node-gem/main.h):
| #include "node.h" class TMyTransport { public: long int nNodes, // Number of mass transport nodes nTimes, // Number of time steps nIC, // Number of chemical independent components nDC, // Number of chemical dependent components nPH, // Number of chemical phases nPS, // Number of chemical phases-solutions nRecipes; // Number of different input node recipes to set boundary conditions long int *aNodeHandle, // Node identification handles *aNodeStatusCH, // Node status codes (changed after GEM calculation) *aIterDone; // Number of GEM IPM iterations performed for each node // at the last time step double *aT, // Array of node temperatures T, Kelvin *aP, // Array of node pressures P, Pa *aVs, // Array of node volume V of reactive subsystem, m3 *aMs, // Array of node mass of reactive subsystem, kg *aGs, // Array of node total Gibbs energy of reactive subsystems, J *aHs, // Array of node total enthalpy of reactive subsystems, J (reserved) *aIC, // Array of node effective aqueous ionic strengths, molal *apH, // Array of node pH of aqueous solutions *ape, // Array of node pe of aqueous solutions *aEh; // Array of node Eh of aqueous solution, V double **axDC, // Array of node mole amounts of dependent components (speciation) **agam, // Array of node activity coefficients of dependent components **axPH, // Array of node total mole amounts of all reactive phases **aaPH, // Array of node specific surface areas of phases, m2/kg **avPS, // Array of node total volumes of multicomponent phases, m3 **amPS, // Array of node total masses of multicomponent phases,kg **abPS, // Array of node bulk compositions of multicomponent phases, moles **axPA, // Array of node amount of carrier in asymmetric phases, moles **aaPh, // Array of node surface areas of phases, m2 **adul, // Array of node upper restrictions to amounts of dependent components **adll, // Array of node lower restrictions to amounts of dependent components **abIC, // Array of node bulk mole amounts of independent components **arMB, // Array of node mole balance residuals for independent components **auIC, // Array of node chemical potentials of independent components (norm.) **abSP; // Array for bulk composition of solid part of equilibrated sub-system TMyTransport( // Constructor (dynamic memory allocation) long int p_nNod, // Number of nodes long int p_nTim, // Number of time steps long int p_nIC, // Number of chemical independent components long int p_nDC, // Number of chemical dependent components long int p_nPH, // Number of chemical phases long int p_nPS, // Number of chemical phases - solutions long int p_nRcps // Number of different input node recipes to set boundary conditions ); ~TMyTransport(); // Destructor of dynamic memory void OneTimeStepRun( // Placeholder function for one transport time step double *stoich, // Stoichiometry coefficients long int *ICndx, // Indexes of mobile independent components long int nICndx // Number of mobile independent components ); }; |
All data and methods in this class are declared as 'public' to make the example program more straightforward. The three functions belonging to TMyTransport class are implemented in the "standalone/node-gem/main.cpp" file as follows:
| #include "node.h" #include "main.h" TMyTransport::TMyTransport( long int p_nNod, long int p_nTim, long int p_nIC, long int p_nDC, long int p_nPH, long int p_nPS, long int p_nRcps ) { nNodes = p_nNod; nTimes = p_nTim; nIC = p_nIC; nDC = p_nDC; nPH = p_nPH; nPS = p_nPS; nRecipes = p_nRcps; aNodeHandle = new long int [nNodes]; aNodeStatusCH = new long int [nNodes]; aIterDone = new long int [nNodes]; aT = new double [nNodes]; aP = new double [nNodes]; aVs = new double [nNodes]; aMs = new double [nNodes]; aGs = new double [nNodes]; aHs = new double [nNodes]; aIC = new double [nNodes]; apH = new double [nNodes]; ape = new double [nNodes]; aEh = new double [nNodes]; axDC = new double *[nNodes]; agam = new double *[nNodes]; axPH = new double *[nNodes]; aaPH = new double *[nNodes]; avPS = new double *[nNodes]; amPS = new double *[nNodes]; abPS = new double *[nNodes]; axPA = new double *[nNodes]; aaPh = new double *[nNodes]; adul = new double *[nNodes]; adll = new double *[nNodes]; abIC = new double *[nNodes]; arMB = new double *[nNodes]; auIC = new double *[nNodes]; abSP = new double *[nNodes]; for (long int in=0; in<nNodes; in++) { abIC[in] = new double [nIC]; arMB[in] = new double [nIC]; auIC[in] = new double [nIC]; axDC[in] = new double [nDC]; agam[in] = new double [nDC]; adul[in] = new double [nDC]; adll[in] = new double [nDC]; aaPH[in] = new double [nPH]; axPH[in] = new double [nPH]; avPS[in] = new double [nPS]; amPS[in] = new double [nPS]; axPA[in] = new double [nPS]; aaPh[in] = new double [nPH]; abPS[in] = new double [nIC*nPS]; abSP[in] = new double [nIC]; } } TMyTransport::~TMyTransport() { // Deleting chemical data arrays for nodes for (long int in=0; in<nNodes; in++) { delete[]abIC[in]; delete[]arMB[in]; delete[]auIC[in]; delete[]axDC[in]; delete[]agam[in]; delete[]adul[in]; delete[]adll[in]; delete[]aaPH[in]; delete[]axPH[in]; delete[]avPS[in]; delete[]amPS[in]; delete[]axPA[in]; delete[]aaPh[in]; delete[]abPS[in]; delete[]abSP[in]; } delete[]axDC; delete[]agam; delete[]axPH; delete[]aaPH; delete[]avPS; delete[]amPS; delete[]abPS; delete[]axPA; delete[]aaPh; delete[]adul; delete[]adll; delete[]abIC; delete[]arMB; delete[]auIC; delete[]abSP; delete[]aNodeHandle; delete[]aNodeStatusCH; delete[]aIterDone; delete[]aT; delete[]aP; delete[]aVs; delete[]aMs; delete[]aGs; delete[]aHs; delete[]aIC; delete[]apH; delete[]ape; delete[]aEh; } // A very simple example of transport algorithm void TMyTransport::OneTimeStepRun( double *stoicf, long int *ICndx, long int nICndx ) { double parcel[nICndx]; long int in; for( in=1; in< nNodes; in++ ) { // some operators that change in some nodes some amounts of some migrating // chemical species (axDC array or some amounts of migrating chemical // elements in abIC array), possibly using data from other arrays in such // a way that the mass conservation within the whole array of nodes is retained for (int i=0; i<nICndx; i++) { parcel[i] = 0.01* stoicf[i]*abIC[in-1][ICndx[i]]; abIC[in][ICndx[i]] += parcel[i]; // if( in > 1 ) // abIC[in-1][ICndx[i]] -= parcel[i]; } // The above example loop implements a zero-order flux of MgCl2 in one direction. // Real advective/diffusive transport models are much more complex, but essentially // do similar things } } |
#include "node.h" #include "main.h" //The case of data exchange in computer memory int main( int argc, char* argv[] ) { char (*recipes)[fileNameLength] = 0; // Analyzing command line arguments ( Default arguments) char input_system_file_list_name[256] = "system-dat.lst"; char input_recipes_file_list_name[256] = "more_recipes.lst"; if (argc >= 2 ) strncpy( input_system_file_list_name, argv[1], 256); // list of DCH, IPM and DBR input files for initializing GEMS3K if (argc >= 3 ) // list of DBR files for setting boundary conditions strncpy( input_recipes_file_list_name, argv[2], 256); // Creates TNode structure instance accessible trough the "node" pointer TNode* node = new TNode(); // (1) Initialization of GEMS3K internal data by reading files // whose names are given in the ipm_input_system_file_list_name if( node->GEM_init( input_system_file_list_name ) ) { cout << "Error occured during reading the files" ; return 1; } // Getting direct access to work node DATABR structure which exchanges the // data with GEMS3K (already filled out by reading the DBR input file) DATABR* dBR = node->pCNode(); // Getting direct access to DataCH structure in GEMS3K instance memory DATACH* dCH = node->pCSD(); // Creating memory for mass transport nodes // 11 nodes, 99 time steps TMyTransport mt( 11, 100, dCH->nICb, dCH->nDCb, dCH->nPHb, dCH->nPSb, 1 ); |
| // Initialization of GEMS3K and chemical information for nodes kept in the MT part long int in; for( in=0; in< mt.nNodes; in++ ) { // Asking GEM IPM to run with automatic initial approximation dBR->NodeStatusCH = NEED_GEM_AIA; // (2) re-calculating equilibrium by calling GEMIPM2K, getting the status back mt.aNodeStatusCH[in] = node->GEM_run( false); if( !( mt.aNodeStatusCH[in] == OK_GEM_AIA || mt.aNodeStatusCH[in] == OK_GEM_SIA ) ) { cout << "Error occured during re-calculating equilibrium" ; return 5; } // Extracting GEM IPM input data to mass-transport program arrays node->GEM_restore_MT( mt.aNodeHandle[in], mt.aNodeStatusCH[in], mt.aT[in], mt.aP[in], mt.aVs[in], mt.aMs[in], mt.abIC[in], mt.adul[in], mt.adll[in], mt.aaPH[in] ); // Extracting GEM IPM output data to mass-transport program arrays node->GEM_to_MT( mt.aNodeHandle[in], mt.aNodeStatusCH[in], mt.aIterDone[in], mt.aVs[in], mt.aMs[in], mt.aGs[in], mt.aHs[in], mt.aIC[in], mt.apH[in], mt.ape[in], mt.aEh[in], mt.arMB[in], mt.auIC[in], mt.axDC[in], mt.agam[in], mt.axPH[in], mt.avPS[in], mt.amPS[in], mt.abPS[in], mt.axPA[in], mt.aaPh[in], mt.abSP[in] ); // Here the setup of initial differences between node compositions, // temperatures, etc. can be implemented // // Here the file output for the initial conditions can be implemented } if (argc >= 3 ) // Read DATABR structure from text file { char NextRecipeFileName[256]; char NextRecipeOutFileName[300]; char input_recipes_file_list_path[256-fileNameLength] = ""; // Reading list of recipes names from file recipes = f_getfiles( input_recipes_file_list_name, input_recipes_file_list_path, mt.nRecipes, ','); // in this example, nRecipes = 1 (one additional input recipe) for( in=0; in<min(mt.nRecipes, mt.nNodes); in++ ) { // Trying to read the next DBR file name sprintf(NextRecipeFileName , "%s%s", input_recipes_file_list_path, recipes[in] ); // (5) Reading the next DBR file (boundary condition on the left) node->GEM_read_dbr( NextRecipeFileName ); // Asking GEM IPM to run with automatic initial approximation dBR->NodeStatusCH = NEED_GEM_AIA; // (2) Re-calculating chemical equilibrium by calling GEM mt.aNodeStatusCH[in] = node->GEM_run( false ); if( !( mt.aNodeStatusCH[in] == OK_GEM_AIA || mt.aNodeStatusCH[in] == OK_GEM_SIA ) ) { cout << "Error occured during re-calculating chemical equilibrium" ; return 5; } // (6) Extracting GEMIPM input data to mass-transport program arrays node->GEM_restore_MT( mt.aNodeHandle[in], mt.aNodeStatusCH[in], mt.aT[in], mt.aP[in], mt.aVs[in], mt.aMs[in], mt.abIC[in], mt.adul[in], mt.adll[in], mt.aaPH[in] ); // (7) Extracting GEMIPM output data to mass-transport program arrays node->GEM_to_MT( mt.aNodeHandle[in], mt.aNodeStatusCH[in], mt.aIterDone[in], mt.aVs[in], mt.aMs[in], mt.aGs[in], mt.aHs[in], mt.aIC[in], mt.apH[in], mt.ape[in], mt.aEh[in], mt.arMB[in], mt.auIC[in], mt.axDC[in], mt.agam[in], mt.axPH[in], mt.avPS[in], mt.amPS[in], mt.abPS[in], mt.axPA[in], mt.aaPh[in], mt.abSP[in] ); // Here the setup of initial differences between node compositions, // temperatures, etc. can be implemented // // Here the file output for the initial conditions can be implemented } // end loop on in } |
The last code snippet below shows how the main loop of coupled
calculations is organized.
Firstly, the use of access methods for extracting indexes of
independent chemical components and phases is shown. Then, the
stoichiometry of matter to be transported is defined using the stoich[5] work array.
Next, the 'time loop' begins with the 'mass transport' as such
represented as a single mt.OneTimeStepRun() function
call; in principle, in this function, also temperature of each node may
change.
After that, the nested 'chemical equilibriation loop' performs
re-calculation of equilibria in all nodes whose bulk composition had
just been modified at the 'mass transport' stage at this time step.
This loop also demonstrates how to use the 'smart initial
approximation' (SIA) mode of GEM IPM. After some output indicating
changes in chemical speciation at this time step in each node,
the chemical equilibration loop is closed, and control goes to
the mass transport time loop. Note that at this stage, various
back-coupling calculations, i.e. how precipitation/dissolution of some
phases has affected the porosity and transport properties, can be
implemented (in a real coupling, of course).
| // Main loop - iterations over nTimes time steps int xCalcite = node->Ph_name_to_xDB("Calcite"); int xDolomite = node->Ph_name_to_xDB("Dolomite-dis"); int xAq_gen = node->Ph_name_to_xDB("aq_gen"); long int ICndx[5]; ICndx[0] = node->IC_name_to_xDB("Ca"); ICndx[1] = node->IC_name_to_xDB("C"); ICndx[2] = node->IC_name_to_xDB("O"); ICndx[3] = node->IC_name_to_xDB("Mg"); ICndx[4] = node->IC_name_to_xDB("Cl"); // Checking indexes cout << "xCa= " << ICndx[0] << " xC=" << ICndx[1] << " xO=" << ICndx[2] << " xMg=" << ICndx[3] << " xCl=" << ICndx[4] << endl << " xCalcite=" << xCalcite << " xDolomite=" << xDolomite << " xAq_gen=" << xAq_gen << endl << endl; double stoich[5] = { 0., 0., 0., 1., 2. }; // defines what is 'transported' long int it; for( it=0; it< mt.nTimes; it++ ) { cout << "Time step " << it << endl; // Mass transport loop over nodes (not a real transport model) mt.OneTimeStepRun( stoich, ICndx, 5 ); // Chemical equilibration loop over nodes for( in=0; in< mt.nNodes; in++ ) { mt.aNodeHandle[in] = in; mt.aNodeStatusCH[in] = NEED_GEM_SIA; // (8) Setting input data for GEM IPM to use available node speciation as // initial approximation node->GEM_from_MT( mt.aNodeHandle[in], mt.aNodeStatusCH[in], mt.aT[in], mt.aP[in], mt.aVs[in], mt.aMs[in], mt.abIC[in], mt.adul[in], mt.adll[in], mt.aaPH[in], mt.axDC[in], mt.agam[in] ); // (9) Passing current FMT iteration information into the work DATABR structure node->GEM_set_MT( (double)it, 1. ); // Calling GEMIPM calculation mt.aNodeStatusCH[in] = node->GEM_run( true ); if( ( mt.aNodeStatusCH[in] == ERR_GEM_AIA || mt.aNodeStatusCH[in] == ERR_GEM_SIA || mt.aNodeStatusCH[in] == T_ERROR_GEM ) ) { cout << "Error: GEM calculation results are not retrieved. Time step" << it << " node " << in << endl; } else { if( ( mt.aNodeStatusCH[in] == BAD_GEM_AIA || mt.aNodeStatusCH[in] == BAD_GEM_SIA ) ) { cout << "Insufficient quality of GEM solution, but GEM results are retrieved" << it << " node " << in << endl; } else // (7) Extracting GEMIPM output data to FMT part node->GEM_to_MT( mt.aNodeHandle[in], mt.aNodeStatusCH[in], mt.aIterDone[in], mt.aVs[in], mt.aMs[in], mt.aGs[in], mt.aHs[in], mt.aIC[in], mt.apH[in], mt.ape[in], mt.aEh[in],mt.arMB[in], mt.auIC[in], mt.axDC[in], mt.agam[in], mt.axPH[in], mt.avPS[in], mt.amPS[in], mt.abPS[in], mt.axPA[in], mt.aaPh[in], mt.abSP[in] ); } // Here, the output upon completion of the time step is usually implemented // to monitor the coupled simulation or collect results cout << " Node " << in ; cout << ": Aq= " << mt.axPH[in][xAq_gen] << " pH= " << mt.apH[in] << " Calcite= " << mt.axPH[in][xCalcite] << endl; } } // Calculations finished - end time reached // Final output e.g. of total simulation time or of the final distribution of // components and phases in all nodes can be implemented here // deleting GEM IPM and data exchange memory structures delete node; // end of example return 0; } |
<under construction>
Last change: 16.10.2012 DK
Copyright (c) 2012 GEMS Development Team