Jump to content
Geochemist's Workbench Support Forum

exothermic crash


Recommended Posts


I am trying to link the internal heat source to the reaction rate of pyrite to simulate an exothermic reaction. I created a 64 bit library and was able load the function when into X1t. When the simulation starts, the program instantly closes. Any suggestions about how I broke this?




#X1t input file

data = thermo.tdat verify
conductivity = conductivity-USGS.dat
suppress ALL
unsuppress  Hematite K-feldspar Pyrite
interval start at 0 day, fluid = fluid_1
time end at .25 m.y.
length = 250 m
width = 1 m
height = 1 m
Nx = 25
discharge start = .307 m3/m2/yr
diffusion_coef = 1e-7 cm2/s
temperature = 100 C
heat_source = file D:\Model\Teena\Heat\heat.dll:heat J/cm3/yr transient, Tmax = 300 C
scope = initial
   H2O          = 1 free kg
   Na+          = 1 g/kg
   balance on Cl-
   Fe++         = 1 mg/kg
   swap Pyrite for O2(aq)
   Pyrite       = 10 free volume%
   SO4--        = 20 mg/kg
   pH           = 6.5
scope = fluid_1
   H2O          = 1 free kg
   Na+          = 1 g/kg
   balance on Cl-
   Fe++         = 1 mg/kg
   swap e- for O2(aq)
   Eh           = .2 V
   SO4--        = 20 mg/kg
   pH           = 6
kinetic Pyrite pre-exp = 3.02e-12 act_eng = 56900 surface = 1e4


// Heatlib.cpp : Defines the functions for the library.
#define EXPORT extern "C" __declspec(dllexport)
#define NOMINMAX

#include <stdio.h>
#include <math.h>
#include "gwb_context.h"

// heat:  Examine the amount of pyrite oxidized and return a heat source
//      proportional to the pyrite reaction rate and the time step.

EXPORT double heat(int i, Nodal_block &n, GWBcontext &c)
    double pre, act_eng, qnk, surf, temp;
    pre = n.rct[0].preexp;
    act_eng = n.rct[0].acteng;
    qnk =  n.rct[0].QoverK;
    surf = n.rct[0].surfac;
    temp = n.tempc + 273.15;
    if (qnk > 1.0)
        return 0.0;
    //This is the reaction rate in mol/s which needs to be in J/cm3/s
    //so divide by the domain size and number of nodal blocks
    //which gets mol/cm3/s and then multiply by the energy of the reaction
    //in J/mol which I took from https://doi.org/10.3390/min11060565
    return pre*surf*(1.0-qnk)*exp(-act_eng/8.314/temp)/c.Length/c.Width/c.Height/c.Nnode*1409000.0;

Link to comment
Share on other sites

Hello Peter,

Thank you for the script. To troubleshoot, could you double check that the header files used to compile the script are the same version as the version of the software installed?

The dll should be recompiled every time the code and function is updated. If the above doesn't help to resolve the issue, you can try writing a function that returns a value of 0 to check if the dll is compiling correctly for the X1t to read in.

Hope this helps,
Jia Wang
Aqueous Solutions LLC

Link to comment
Share on other sites



This computer has never had anything but the GWB 17 installed on it, the current version is 17.0.2. Running the dll that instantly returns 0 works fine. If I put the return 0 after:

 pre = n.rct[0].preexp;

it crashes. It seems to work fine if I put:

   if (c.Xi == 0)
      return 0;

at the beginning so it doesn't make a calculation for the first time step. I assume that means that the first heat transport step is done before the reaction properties are initialized.




Link to comment
Share on other sites

Hi Peter,

Since your function is called during the domain set up at the start of the run before the reactants are laid out in memory you will need to check either c.Xi as you did or you could also check n.rct.size() before attempting to access a particular reactant.

Once time stepping begins your heat source function will be called before heat transport is done.


Hope this helps,




Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

  • Create New...