exothermic crash

[p.m.berger]

Hello,

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?

Thanks,

Peter

 

#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
scope
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;
}

[Jia Wang]

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

[p.m.berger]

Hello,

 

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.

 

Thanks,

Peter

[Jia Wang]

Hello Peter,

Thank you for the additional information. We are looking into this issue now.

Best regards,
Jia

[dan saalfeld]

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,

Dan