Solubility.py

import os
from modena import *
from modena import ForwardMappingModel,BackwardMappingModel,SurrogateModel,\
    CFunction,IndexSet
import modena.Strategy as Strategy
from fireworks.user_objects.firetasks.script_task import FireTaskBase, ScriptTask
from fireworks import Firework, Workflow, FWAction
from fireworks.utilities.fw_utilities import explicit_serialize
from blessings import Terminal
from jinja2 import Template

# Create terminal for colour output
term = Terminal()

species = IndexSet(
    name= 'solubility_pol_species',
    names= [ 'Air', 'CO2', 'CyP', 'O2', 'N2' ]
)
system = IndexSet(
    name= 'solubility_num_of_components',
    names= [ '2', '3' ]
)

@explicit_serialize

class SolubilityExactSim(ModenaFireTask):
    def task(self, fw_spec):
        # Write input for detailed model
        ff = open('in.txt', 'w')
        Tstr = str(self['point']['T'])
        ff.write('%s \n' %(Tstr))
        if (self['indices']['B']=='2'):
            ff.write('2 \n')       #number of components in system
            if (self['indices']['A']=='CO2'):
                ff.write('co2 \n')
            elif (self['indices']['A']=='Air'):
                ff.write('air \n')
            elif (self['indices']['A']=='CyP'):
                ff.write('cyclopentane \n')
            elif (self['indices']['A']=='O2'):
                ff.write('o2 \n')
            elif (self['indices']['A']=='N2'):
                ff.write('n2 \n')
            #pass molar liquid composition
            x1l_str = str(self['point']['xl1'])
            x2l_str = str(self['point']['xl2'])
            ff.write('%s \n' %(x1l_str))
            ff.write('%s \n' %(x2l_str))
        elif (self['indices']['B']=='3'):
            ff.write('3 \n')       #number of components in system
            if (self['indices']['A']=='CO2'):
                ff.write('co2 \n')
            elif (self['indices']['A']=='Air'):
                ff.write('air \n')
            elif (self['indices']['A']=='CyP'):
                ff.write('cyclopentane \n')
            elif (self['indices']['A']=='O2'):
                ff.write('o2 \n')
            elif (self['indices']['A']=='N2'):
                ff.write('n2 \n')
            #pass molar liquid composition
            x1l_str = str(self['point']['xl1'])
            x2l_str = str(self['point']['xl2'])
            x3l_str = str(self['point']['xl3'])
            ff.write('%s \n' %(x1l_str))
            ff.write('%s \n' %(x2l_str))
            ff.write('%s \n' %(x3l_str))
        ff.close()

        #create output file for detailed code
        fff = open('out.txt', 'w+')
        fff.close()

        # Execute the detailed model
        # path to **this** file + /src/...
        # will break if distributed computing
        ret = os.system(os.path.dirname(os.path.abspath(__file__))+\
            '/src/pcsaft')
        # This call enables backward mapping capabilities
        self.handleReturnCode(ret)
        # Analyse output
        # os.getcwd() returns the path to the "launcher" directory
        try:
            FILE = open(os.getcwd()+'/out.txt','r')        except IOError:
            raise IOError("File not found")
        self['point']['H'] = float(FILE.readline())
        FILE.close()


Ccode2='''
#include "modena.h"
#include "math.h"

void surroSolubility2
(
const modena_model_t* model,
const double* inputs,
double *outputs
)
{
{% block variables %}{% endblock %}

const double P0 = parameters[0];
const double P1 = parameters[1];
const double P2 = parameters[2];

const double term1 = P1*(1/T - 1/P2);
const double term2 = exp(term1);

outputs[0] = P0*term2;

outputs[0] = P0 + T*P1 + P2*T*T;
}
'''
Ccode3='''
#include "modena.h"
#include "math.h"

void surroSolubility3
(
const modena_model_t* model,
const double* inputs,
double *outputs
)
{
{% block variables %}{% endblock %}

const double P0 = parameters[0];
const double P1 = parameters[1];
const double P2 = parameters[2];
const double P3 = parameters[3];

const double term1 = P1*(1/T - 1/P2);
const double term2 = exp(term1);

outputs[0] = P0*term2;

outputs[0] = P0 + T*P1 + P2*T*T;

outputs[0] = P0+P1*exp(-pow((T-P2),2)/P3);
}
'''
outputs={
    'H': { 'min': 9e99, 'max': -9e99, 'argPos': 0 }
}
parameters={
    'param0': { 'min': 1e-12, 'max': 1E1, 'argPos': 0 },    #check if boundaries are reasonable!!!
    'param1': { 'min': -1e-1, 'max': 1e-1, 'argPos': 1 },
    'param2': { 'min': -1e-1, 'max': 1e-1, 'argPos': 2 },
}
parameters4={
    'param0': { 'min': -9e9, 'max': 9e9+2*1.42885440e-04, 'argPos': 0 },    #check if boundaries are reasonable!!!
    'param1': { 'min': -9e9, 'max': 9e9+2*1.22132172e+00, 'argPos': 1 },
    'param2': { 'min': -9e9, 'max': 9e9+2*-7.97789449e+02, 'argPos': 2 },
    'param3': { 'min': -9e9, 'max': 9e9+2*1.33835999e+05, 'argPos': 3 },
}
indices={
    'A': species,
    'B': system
}
inputs2={
    'T': { 'min': 200.0, 'max': 550.0},        #check if boundaries reasonable, from this range, the random values for the DOE are chosen!
    'xl1': { 'min': 0.0, 'max': 1.0 },
    'xl2': { 'min': 0.0, 'max': 1.0 },
}
inputs3={
    'T': { 'min': 200.0, 'max': 550.0},        #check if boundaries reasonable, from this range, the random values for the DOE are chosen!
    'xl1': { 'min': 0.0, 'max': 1.0 },
    'xl2': { 'min': 0.0, 'max': 1.0 },
    'xl3': { 'min': 0.0, 'max': 1.0 },
}
f2 = CFunction(Ccode=Ccode2,
    inputs=inputs2,
    outputs=outputs,
    parameters=parameters,
    indices=indices
)
f3 = CFunction(Ccode=Ccode3,
    inputs=inputs3,
    outputs=outputs,
    parameters=parameters4,
    indices=indices
)

outOfBoundsStrategy=Strategy.ExtendSpaceStochasticSampling(
    nNewPoints=4
)
parameterFittingStrategy=Strategy.NonLinFitWithErrorContol(
    testDataPercentage=0.2,
    maxError=0.1,
    improveErrorStrategy=Strategy.StochasticSampling(
        nNewPoints=2
    ),
    maxIterations=5  # Currently not used
)
m_solubilityCO2PU = BackwardMappingModel(
    _id='Solubility[A=CO2,B=2]',
    surrogateFunction=f2,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 320, 350, 380],
            'xl1': [1.1e-3, 1.0e-3, 1.0e-3, 1.0e-4],
            'xl2': [0.9989, 0.999, 0.999, 0.9999],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
m_solubilityAirPU = BackwardMappingModel(
    _id='Solubility[A=Air,B=2]',
    surrogateFunction=f2,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 320, 350, 380],
            'xl1': [1.1e-3, 1.0e-3, 1.0e-3, 1.0e-4],
            'xl2': [0.9989, 0.999, 0.999, 0.9999],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
m_solubilityCyclopentanePU = BackwardMappingModel(
    _id='Solubility[A=CyP,B=2]',
    surrogateFunction=f2,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 320, 350, 380],
            'xl1': [1.1e-3, 1.0e-3, 1.0e-3, 1.0e-4],
            'xl2': [0.9989, 0.999, 0.999, 0.9999],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
m_solubilityO2PU = BackwardMappingModel(
    _id='Solubility[A=O2,B=2]',
    surrogateFunction=f2,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 320, 350, 380],
            'xl1': [1.1e-3, 1.0e-3, 1.0e-3, 1.0e-4],
            'xl2': [0.9989, 0.999, 0.999, 0.9999],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
m_solubilityN2PU = BackwardMappingModel(
    _id='Solubility[A=N2,B=2]',
    surrogateFunction=f2,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 320, 350, 380],
            'xl1': [1.1e-3, 1.0e-3, 1.0e-3, 1.0e-4],
            'xl2': [0.9989, 0.999, 0.999, 0.9999],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
# 3 component models
m_solubilityCO2 = BackwardMappingModel(
    _id='Solubility[A=CO2,B=3]',
    surrogateFunction=f3,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 350, 450, 550],
            'xl1': [0.9e-4, 1.1e-4, 1e-4, 1e-4],
            'xl2': [0.51, 0.49, 0.5, 0.5],
            'xl3': [0.49, 0.51, 0.5, 0.5],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
m_solubilityAir = BackwardMappingModel(
    _id='Solubility[A=Air,B=3]',
    surrogateFunction=f3,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 350, 450, 550],
            'xl1': [0.9e-4, 1.1e-4, 1e-4, 1e-4],
            'xl2': [0.51, 0.49, 0.5, 0.5],
            'xl3': [0.49, 0.51, 0.5, 0.5],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
m_solubilityCyclopentane = BackwardMappingModel(
    _id='Solubility[A=CyP,B=3]',
    surrogateFunction=f3,
    exactTask=SolubilityExactSim(),
    substituteModels=[],
    initialisationStrategy=Strategy.InitialPoints(
        initialPoints={
            'T': [290, 350, 450, 550],
            'xl1': [0.9e-4, 1.1e-4, 1e-4, 1e-4],
            'xl2': [0.51, 0.49, 0.5, 0.5],
            'xl3': [0.49, 0.51, 0.5, 0.5],
        },
    ),
    outOfBoundsStrategy=outOfBoundsStrategy,
    parameterFittingStrategy=parameterFittingStrategy
)
# below are experimental solubility models
# they are needed, because we want substitute model for bubble growth model
inputsExp={
    'T': { 'min': 200.0, 'max': 550.0}
}
CcodeR11Baser='''
#include "modena.h"
#include "math.h"

void surroSolubilityR11Baser
(
const modena_model_t* model,
const double* inputs,
double *outputs
)
{
{% block variables %}{% endblock %}

const double P0 = parameters[0];
const double P1 = parameters[1];
const double P2 = parameters[2];
const double P3 = parameters[3];

outputs[0] = (P0 + P1*exp(-pow(T-P2,2)/(2*P3*P3)))*1100.0/137.37e-3/101e3;
}
'''
parameters4={
    'param0': { 'min': 1e-12, 'max': 1E1, 'argPos': 0 },
    'param1': { 'min': -1e-1, 'max': 1e-1, 'argPos': 1 },
    'param2': { 'min': -1e-1, 'max': 1e-1, 'argPos': 2 },
    'param3': { 'min': -1e-1, 'max': 1e-1, 'argPos': 3 },
}
fR11Baser = CFunction(Ccode=CcodeR11Baser,
    inputs=inputsExp,
    outputs=outputs,
    parameters=parameters4
)

parR11Baser = [1e-7, 4.2934, 203.3556, 40.016]
m_solubilityR11Baser = ForwardMappingModel(
    _id='SolubilityR11Baser',
    surrogateFunction=fR11Baser,
    substituteModels=[],
    parameters=parR11Baser,
    inputs=inputsExp,
    outputs=outputs,
)
CcodeCO2Baser='''
#include "modena.h"
#include "math.h"

void surroSolubilityCO2Baser
(
const modena_model_t* model,
const double* inputs,
double *outputs
)
{
{% block variables %}{% endblock %}

const double P0 = parameters[0];

outputs[0] = P0;
}
'''
parameters1={
    'param0': { 'min': 1e-12, 'max': 1E1, 'argPos': 0 },
}
if 'argPos' in inputsExp['T']: #added by previous CFunction
    inputsExp['T'].pop('argPos')
fCO2Baser = CFunction(Ccode=CcodeCO2Baser,
    inputs=inputsExp,
    outputs=outputs,
    parameters=parameters1
)

parCO2Baser = [1.1e-4]
m_solubilityCO2Baser = ForwardMappingModel(
    _id='SolubilityCO2Baser',
    surrogateFunction=fCO2Baser,
    substituteModels=[],
    parameters=parCO2Baser,
    inputs=inputsExp,
    outputs=outputs,
)
CcodePentGupta='''
#include "modena.h"
#include "math.h"

void surroSolubilityPentGupta
(
const modena_model_t* model,
const double* inputs,
double *outputs
)
{
{% block variables %}{% endblock %}

const double P0 = parameters[0];
const double P1 = parameters[1];
const double P2 = parameters[2];
const double P3 = parameters[3];
const double P4 = parameters[4];

outputs[0] = P0/(exp((P1-P2*T)/(P3-T))-P4)*1100.0/72.15e-3/101e3;
}
'''
parameters5={
    'param0': { 'min': 1e-12, 'max': 1E1, 'argPos': 0 },
    'param1': { 'min': -1e-1, 'max': 1e-1, 'argPos': 1 },
    'param2': { 'min': -1e-1, 'max': 1e-1, 'argPos': 2 },
    'param3': { 'min': -1e-1, 'max': 1e-1, 'argPos': 3 },
    'param4': { 'min': -1e-1, 'max': 1e-1, 'argPos': 4 },
}
if 'argPos' in inputsExp['T']: #added by previous CFunction
    inputsExp['T'].pop('argPos')
fPentGupta = CFunction(Ccode=CcodePentGupta,
    inputs=inputsExp,
    outputs=outputs,
    parameters=parameters5
)

parPentGupta = [-3.3e-4, 2.09e4, 67.5, 8.69e4, 1.01]
m_solubilityPentGupta = ForwardMappingModel(
    _id='SolubilityPentGupta',
    surrogateFunction=fPentGupta,
    substituteModels=[],
    parameters=parPentGupta,
    inputs=inputsExp,
    outputs=outputs,
)
CcodePentWinkler='''
#include "modena.h"
#include "math.h"

void surroSolubilityPentWinkler
(
const modena_model_t* model,
const double* inputs,
double *outputs
)
{
{% block variables %}{% endblock %}

const double P0 = parameters[0];
const double P1 = parameters[1];
const double P2 = parameters[2];
const double P3 = parameters[3];

outputs[0] = (P0 + P1*exp(-pow(T-P2,2)/(2*P3*P3)))*1100.0/72.15e-3/101e3;
}
'''
if 'argPos' in inputsExp['T']: #added by previous CFunction
    inputsExp['T'].pop('argPos')
fPentWinkler = CFunction(Ccode=CcodePentWinkler,
    inputs=inputsExp,
    outputs=outputs,
    parameters=parameters5
)

parPentWinkler = [0.0064,0.0551,298.0,17.8]
m_solubilityPentWinkler = ForwardMappingModel(
    _id='SolubilityPentWinkler',
    surrogateFunction=fPentWinkler,
    substituteModels=[],
    parameters=parPentWinkler,
    inputs=inputsExp,
    outputs=outputs,
)