Tea: different behavior when computing properties together

[seb007]

Hi,

I encountered a strange behavior when calling Tea: It provides dirrerent values if I compute one or several properties simultaneously. Here is a "minimal" example (interesting part is the end after #ifdef):

Code: Select all

	HRESULT hr = CoInitialize(NULL);
	if (FAILED(hr)) messages::get().error("Cannot initialize COM: Error %i", hr);

	_mat = NULL;
	hr = ::CoCreateInstance(__uuidof(fpCapeOpenCOM),
		NULL,
		CLSCTX_INPROC_SERVER,
		__uuidof(ICapeThermoMaterial),
		(void**)&_mat); 

	if (FAILED(hr)) messages::get().error("Cannot instantiate CAPE-OPEN material object: Error %i", hr);

	GUID ppmId;
	hr = CLSIDFromString(OLESTR("COCO_TEA.PropertyPackManager.1"), (LPCLSID)&ppmId);

	//hr = CLSIDFromString(OLESTR("DWSIM.CAPEOPENPropertyPackageManager"), (LPCLSID)&teaId);
	if (FAILED(hr)) messages::get().error("Cannot find property package id: Error %i", hr);

	hr = ::CoCreateInstance(ppmId,
		NULL,
		CLSCTX_INPROC_SERVER,
		__uuidof(ICapeThermoPropertyPackageManager),
		(void**)&_ppManager);
	if (FAILED(hr)) messages::get().error("Cannot instantiate CAPE-OPEN property package manager and get interface: Error %i", hr);


	/*CComVariant propPackages;
	hr = _ppManager->GetPropertyPackageList(&propPackages);
	if (FAILED(hr)) messages::get().error("Cannot get property packages list: Error %i", hr);
	CComSafeArray <BSTR> propPackagesArray;
	propPackagesArray.Attach(propPackages.parray);
	messages::get().info("Number of available property packages %i", propPackagesArray.GetCount());
	for (int i = 0; i < propPackagesArray.GetCount(); i++){
		CComBSTR mElement(propPackagesArray.GetAt(i));
		messages::get().info("     ---- (%i): %s", i, _com_util::ConvertBSTRToString(mElement.m_str));
	}
	propPackagesArray.Detach();*/

	hr = _ppManager->GetPropertyPackage(L"Longpipe", (IDispatch**)&_pp);
	//hr = _ppManager->GetPropertyPackage(L"RapidDecompression", (IDispatch**)&_pp);

	if (FAILED(hr)) messages::get().error("Cannot get property package : Error %i", hr);


	hr = _pp->QueryInterface(__uuidof(ICapeThermoMaterialContext), (LPVOID*)&_matContext);
	if (FAILED(hr)) messages::get().error("Cannot get material context from property package : Error %i", hr);

	hr = _matContext->SetMaterial(_mat);
	if (FAILED(hr)) messages::get().error("Cannot set material to material context : Error %i", hr);

	hr = _pp->QueryInterface(__uuidof(ICapeThermoPropertyRoutine), (LPVOID*)&_thermoPropRoutine);
	if (FAILED(hr)) messages::get().error("Cannot get property routine from property package : Error %i", hr);

	//Setting phases
	CComSafeArray <BSTR> phases(1);
	phases[0] = L"Vapor";
	CComVariant phasesVar = phases;
	CComSafeArray <int> status(1);
	status[0] = CAPE_UNKNOWNPHASESTATUS;
	CComVariant statusVar = status;
	hr = _mat->SetPresentPhases(phasesVar, statusVar);
	if (FAILED(hr)) messages::get().error("Cannot set phase");

	//Setting temperature
	CComSafeArray <double> currentTArr(1);
	CComSafeArray <double> currentPArr(1);
	currentTArr[0] = 270;
	currentPArr[0] = 1.0e6;
	CComVariant currentTVar = currentTArr;
	CComVariant currentPVar = currentPArr;
	hr = _mat->SetSinglePhaseProp(L"temperature", L"Vapor", NULL, currentTVar);
	if (FAILED(hr)) messages::get().error("Cannot set temperature");
	hr = _mat->SetSinglePhaseProp(L"pressure", L"Vapor", NULL, currentPVar);
	if (FAILED(hr)) messages::get().error("Cannot set pressure");

	CComSafeArray <double> fraction(10);
	fraction[0] = 0.0137;
	fraction[1] = 0.0138;
	fraction[2] = 0.8742;
	fraction[3] = 0.0697;
	fraction[4] = 0.0213;
	fraction[5] = 0.003;
	fraction[6] = 0.0025;
	fraction[7] = 0.0006;
	fraction[8] = 0.0003;
	fraction[9] = 0.0009;
	CComVariant fractionVar = fraction;
	hr = _mat->SetSinglePhaseProp(L"fraction", L"Vapor", L"Mole", fractionVar);

	if (FAILED(hr)) messages::get().error("Cannot set mole fraction");

#if TEST==1
	printf("test1\n");
	int nbProps = 5;
	CComSafeArray <BSTR> props(nbProps);
	props[0] = L"density";
	props[1] = L"internalEnergy";
	props[2] = L"internalEnergy.Dtemperature";
	props[3] = L"internalEnergy.Dpressure";
	props[4] = L"molecularWeight";
	CComVariant propsVar = props;
#elif TEST==2
	printf("test2\n");
	int nbProps = 3;
	CComSafeArray <BSTR> props(nbProps);
	props[0] = L"density";
	props[1] = L"internalEnergy";
	props[2] = L"molecularWeight"; 
	CComVariant propsVar = props;
#elif TEST==3
	printf("test3\n");
	int nbProps = 2;
	CComSafeArray <BSTR> props(nbProps);
	props[0] = L"density";
	props[1] = L"internalEnergy";
	CComVariant propsVar = props;
#endif

	hr = _thermoPropRoutine->CalcSinglePhaseProp(propsVar, L"Vapor");
	if (FAILED(hr)) messages::get().error("Cannot compute properties: Error %i",hr);
	CComVariant propGet;
	CComSafeArray <double> propGetArr;
	for (int iProp = 0; iProp < nbProps; iProp++){
		hr = _mat->GetSinglePhaseProp(props[iProp], L"Vapor", NULL, &propGet);
		if (FAILED(hr)) messages::get().error("Cannot get property %ws", props[iProp]);
		propGetArr.Attach(propGet.parray);
		printf("Property %ws is %.14e\n",props[iProp],propGetArr[0]);
		propGetArr.Detach();
	}
	_thermoPropRoutine->Release();
	_matContext->Release();
	_pp->Release();
	_ppManager->Release();
	_mat->Release();
	CoUninitialize();

The three different output while changing TEST define value (look at the value of "internalEnergy"):

Code: Select all

test1
Property density is 4.63068260659736e+002
Property internalEnergy is -3.48112448296645e+003
Property internalEnergy.Dtemperature is 2.96241259172201e+001
Property internalEnergy.Dpressure is -1.80752517698902e-004
Property molecularWeight is 1.85132287300110e+001

Code: Select all

test2
Property density is 4.63068260659736e+002
Property internalEnergy is -3.48112533568150e+003
Property molecularWeight is 1.85132287300110e+001

Code: Select all

test3
ERROR : Cannot compute properties: Error -2147220223
ERROR -99

I use the standard Peng-Robinson template with (ordered) Nitrogen, Carbon Dioxide, Methane, Ethane, Propane, Isobutane, N-butane, Isopentane, N-pentane, N-hexane. It may come from my implementation of material object, but it is rather simple. Also, I thought of a tolerance issue during solver iterations, but there is no iteration for this kind of properties, right (i.e. analytical properties AND derivatives)? And the test3 results is even more strange.

Any idea where it can come from? Is it possible to test multiple properties computation in COFE calculator (I did not find a way).

Thanks

[jasper]

You do not need to set the phases for single phase property calculations. Only prior to phase equilibrium calculations do you want to indicate which phases are allowed to take part in the equilibrium. Unless of course your material object internally requires the SetPhases call....

When you get back a CAPE-OPEN error code (any of below), you can QI for ECapeRoot and get a textual error description from it.

Code: Select all

enum eCapeErrorInterfaceHR_tag
{
    ECapeUnknownHR = -2147220223,
    ECapeDataHR = -2147220222,
    ECapeLicenceErrorHR = -2147220221,
    ECapeBadCOParameterHR = -2147220220,
    ECapeBadArgumentHR = -2147220219,
    ECapeInvalidArgumentHR = -2147220218,
    ECapeOutOfBoundsHR = -2147220217,
    ECapeImplementationHR = -2147220216,
    ECapeNoImplHR = -2147220215,
    ECapeLimitedImplHR = -2147220214,
    ECapeComputationHR = -2147220213,
    ECapeOutOfResourcesHR = -2147220212,
    ECapeNoMemoryHR = -2147220211,
    ECapeTimeOutHR = -2147220210,
    ECapeFailedInitialisationHR = -2147220209,
    ECapeSolvingErrorHR = -2147220208,
    ECapeBadInvOrderHR = -2147220207,
    ECapeInvalidOperationHR = -2147220206,
    ECapePersistenceHR = -2147220205,
    ECapeIllegalAccessHR = -2147220204,
    ECapePersistenceNotFoundHR = -2147220203,
    ECapePersistenceSystemErrorHR = -2147220202,
    ECapePersistenceOverflowHR = -2147220201,
    ECapeOutsideSolverScopeHR = -2147220200,
    ECapeHessianInfoNotAvailableHR = -2147220199,
    ECapeThrmPropertyNotAvailableHR = -2147220198
};

Otherwise GetLastError will probably get you a long way (often this requires loading NTDLL.dll, which contains most of the system error messages).

Your strings must be BSTR values, you are assigning const wchar_t * values. This is not correct and will crash under various conditions. You can use CComBSTR to allocate one, or manually use SysAllocString and friends (don't forget to SysFreeString).

The relative error between the two internal energy values is 2e-7, I expect this to be the convergence tolerance of the EOS. The EOS is solved iteratively (the analytic solution can be fraught with very large numerical error).

The COFE calculator does not allow for simultaneous property calculation. However, any custom unit operation will allow you do this. As does for example Matlab or Scilab thermo import.

I do not see the definition of your COM pointers, are they smart pointers? If so, why the Release() at the end? Note that you must Terminate() PMC objects, as well as Initialize() them via ICapeUtilities.

[seb007]

Thanks for these very useful comments. I fixed the BSTR issues (by always using CComBSTR), and called Initialize and Terminate functions.

I do not use smart pointers for COM objects, but normal pointers initialized to NULL in the class constructor and released in the destructor, because they are class members declared in a .h file in which I do not want to include any COM stuff.

I agree that the difference can result from the solver tolerance. Note that the relative error can be as high as 1e-4 for certain values of p and T:

Computing internal energy and its derivatives:
T 2.704082e+002 p 6.816327e+006 U -2.533797e+005

Computing only internal energy:
T 2.704082e+002 p 6.816327e+006 U -2.534055e+005

[jasper]

I agree that the difference can result from the solver tolerance. Note that the relative error can be as high as 1e-4 for certain values of p and T

This is not impossible, the tolerance is on Z. You can test by tightening the tolerance.

Did you get an error message for case 3?

[jasper]

Note that Terminate and Release must be done on 'primary PMC objects'. These are the ones you create (e.g. with CoCreateInstance), but also the ones that are created from a 'factory' object. Particularly property packages created by thermo systems or property package managers. Derived objects (parameters and ports of unit operations for example) do not need to be initialized and terminated.

[seb007]

It looks like the error on TEST3 has been fixed after I made the changes you proposed.

I played with TEA properties: it is not the Newton tolerance, but the value of the temperature and pressure perturbation, used to compute derivatives I suppose.

The strange thing is that I can set them to 1e-14 and get machine precision match between tests 1 and 2. I expected a value of about 1e-7 to be optimal.

[jasper]

The perturbation is used for enthalpy in case the approx enthalpy derivative model is not selected. The analytic derivatives for enthalpy require first order derivatives of (d ln Phi / d T) or (d ln gamma / d T) with respect to T, P and X, so second order derivatives of ln Phi and ln gamma. I have implementations that provide this, but not in TEA; TEA does these by perturbation. Internal energy is calculated from enthalpy, volume and pressure.

This should not affect the precision of the internal energy value itself though - merely its derivatives.

[seb007]

In my case, the internal energy is clearly affected by this parameter, and only IF derivatives are computed (compare with the values in my first post).

Here is what I obtain when I keep the default newton tolerance but change the rel. perturbation of t and p to 0.1:

Code: Select all

XXXXXXXXXXXXXXXXXXXXXXXX - test1
Property density is 4.63068260659736e+002
Property internalEnergy is -3.47260004151000e+003
Property internalEnergy.Dtemperature is 2.96630141387408e+001
Property internalEnergy.Dpressure is -1.81207174324797e-004
Property molecularWeight is 1.85132287300110e+001
XXXXXXXXXXXXXXXXXXXXXXXXX - test2
Property density is 4.63068260659736e+002
Property internalEnergy is -3.48112533568150e+003
Property molecularWeight is 1.85132287300110e+001
XXXXXXXXXXXXXXXXXXXXXXXXX - test3
Property density is 4.63068260659736e+002
Property internalEnergy is -3.48112533568150e+003

Is it possible that, in Tea, when computing the derivative (and then a perturbed value for internal energy), this perturbed variable is then not properly set back to its unperturbed value?

[seb007]

From the following example, it seems that we have a difference when derivative regards to p is computed:

Code: Select all

	HRESULT hr = CoInitialize(NULL);
	if (FAILED(hr)) messages::get().error("Cannot initialize COM: Error %i", hr);

	hr = ::CoCreateInstance(__uuidof(fpCapeOpenCOM),
		NULL,
		CLSCTX_INPROC_SERVER,
		__uuidof(ICapeThermoMaterial),
		(void**)&_mat);
	if (FAILED(hr)) messages::get().error("Cannot instantiate CAPE-OPEN material object: Error %i", hr);

	GUID ppmId;
	hr = CLSIDFromString(OLESTR("COCO_TEA.PropertyPackManager.1"), (LPCLSID)&ppmId);

	//hr = CLSIDFromString(OLESTR("DWSIM.CAPEOPENPropertyPackageManager"), (LPCLSID)&teaId);
	if (FAILED(hr)) messages::get().error("Cannot find property package id: Error %i", hr);

	hr = ::CoCreateInstance(ppmId,
		NULL,
		CLSCTX_INPROC_SERVER,
		__uuidof(ICapeThermoPropertyPackageManager),
		(void**)&_ppManager);
	if (FAILED(hr)) messages::get().error("Cannot instantiate CAPE-OPEN property package manager and get interface: Error %i", hr);

	hr = _ppManager->QueryInterface(__uuidof(ICapeUtilities), (LPVOID*)&_ppMUtil);
	if (FAILED(hr)) messages::get().error("Cannot get ICapeUtilities interface: Error %i", hr);
	hr = _ppMUtil->Initialize();
	if (FAILED(hr)) messages::get().error("Cannot initialize property package manager: Error %i", hr);

	CComBSTR longpipeStr("Longpipe");
	hr = _ppManager->GetPropertyPackage(longpipeStr, (IDispatch**)&_pp);
	//hr = _ppManager->GetPropertyPackage(L"RapidDecompression", (IDispatch**)&_pp);
	if (FAILED(hr)) messages::get().error("Cannot get property package : Error %i", hr);

	hr = _pp->QueryInterface(__uuidof(ICapeUtilities), (LPVOID*)&_ppUtil);
	if (FAILED(hr)) messages::get().error("Cannot get ICapeUtilities interface: Error %i", hr);
	hr = _ppUtil->Initialize();
	if (FAILED(hr)) messages::get().error("Cannot initialize property package: Error %i", hr);

	hr = _pp->QueryInterface(__uuidof(ICapeThermoMaterialContext), (LPVOID*)&_matContext);
	if (FAILED(hr)) messages::get().error("Cannot get material context from property package : Error %i", hr);

	hr = _matContext->SetMaterial(_mat);
	if (FAILED(hr)) messages::get().error("Cannot set material to material context : Error %i", hr);

	hr = _pp->QueryInterface(__uuidof(ICapeThermoPropertyRoutine), (LPVOID*)&_thermoPropRoutine);
	if (FAILED(hr)) messages::get().error("Cannot get property routine from property package : Error %i", hr);



	int nbProps[4] = { 3, 2, 2, 1 };
	CComBSTR vapStr("Vapor");
	CComBSTR tStr("temperature");
	CComBSTR pStr("pressure");
	for (int test = 0; test < 4; test++){
		_mat->ClearAllProps();

		//Setting phases
		CComSafeArray <BSTR> phases(1);
		phases.SetAt(0,CComBSTR("Vapor"), true);
		CComVariant phasesVar = phases;
		CComSafeArray <int> status(1);
		status.SetAt(0, CAPE_UNKNOWNPHASESTATUS, true);
		CComVariant statusVar = status;
		hr = _mat->SetPresentPhases(phasesVar, statusVar);
		if (FAILED(hr)) messages::get().error("Cannot set phase");

		//Setting temperature
		CComSafeArray <double> currentTArr(1);
		CComSafeArray <double> currentPArr(1);
		currentTArr.SetAt(0, 270, true);
		currentPArr.SetAt(0, 1.0e6, true);
		CComVariant currentTVar = currentTArr;
		CComVariant currentPVar = currentPArr;
		hr = _mat->SetSinglePhaseProp(tStr, vapStr, NULL, currentTVar);
		if (FAILED(hr)) messages::get().error("Cannot set temperature");
		hr = _mat->SetSinglePhaseProp(pStr, vapStr, NULL, currentPVar);
		if (FAILED(hr)) messages::get().error("Cannot set pressure");

		CComSafeArray <double> fraction(10);
		fraction[0] = 0.0137;
		fraction[1] = 0.0138;
		fraction[2] = 0.8742;
		fraction[3] = 0.0697;
		fraction[4] = 0.0213;
		fraction[5] = 0.003;
		fraction[6] = 0.0025;
		fraction[7] = 0.0006;
		fraction[8] = 0.0003;
		fraction[9] = 0.0009;
		CComVariant fractionVar = fraction;
		hr = _mat->SetSinglePhaseProp(L"fraction", L"Vapor", L"Mole", fractionVar);
		if (FAILED(hr)) messages::get().error("Cannot set mole fraction");

		CComSafeArray <BSTR> props(nbProps[test]);
		switch (test){
		case 0:
			printf("XXXXXXXXXXXXXXXXXXXXXXXX - test0\n");
			props.SetAt(0, CComBSTR("internalEnergy"), true);
			props.SetAt(1, CComBSTR("internalEnergy.Dtemperature"), true);
			props.SetAt(2, CComBSTR("internalEnergy.Dpressure"), true);
			break;
		case 1:
			printf("XXXXXXXXXXXXXXXXXXXXXXXX - test1\n");
			props.SetAt(0, CComBSTR("internalEnergy"), true);
			props.SetAt(1, CComBSTR("internalEnergy.Dtemperature"), true);
			break;
		case 2:
			printf("XXXXXXXXXXXXXXXXXXXXXXXXX - test2\n");
			props.SetAt(0, CComBSTR("internalEnergy"), true);
			props.SetAt(1, CComBSTR("internalEnergy.Dpressure"), true);
			break;
		case 3:
			printf("XXXXXXXXXXXXXXXXXXXXXXXXX - test3\n");
			props.SetAt(0, CComBSTR("internalEnergy"), true);
			break;
		}

		CComVariant propsVar = props;
		hr = _thermoPropRoutine->CalcSinglePhaseProp(propsVar, vapStr);
		if (FAILED(hr)) messages::get().error("Cannot compute properties: Error %i", hr);
		CComVariant propGet;
		CComSafeArray <double> propGetArr;
		for (int iProp = 0; iProp < nbProps[test]; iProp++){
			hr = _mat->GetSinglePhaseProp(props[iProp], vapStr, NULL, &propGet);
			if (FAILED(hr)) messages::get().error("Cannot get property %ws", props[iProp]);
			propGetArr.Attach(propGet.parray);
			printf("Property %ws is %.14e\n", props[iProp], propGetArr[0]);
			propGetArr.Detach();
		}
	}
	hr = _ppMUtil->Terminate();
	if (FAILED(hr)) messages::get().error("Cannot terminate property package manager: Error %i", hr);
	hr = _ppUtil->Terminate();
	if (FAILED(hr)) messages::get().error("Cannot terminate property package: Error %i", hr);
	_thermoPropRoutine->Release();
	_matContext->Release();
	_pp->Release();
	_ppManager->Release();
	_mat->Release();
	_ppUtil->Release();
	_ppMUtil->Release();
	CoUninitialize();

producing:

Code: Select all

XXXXXXXXXXXXXXXXXXXXXXXX - test0
Property internalEnergy is -3.47260004151000e+003
Property internalEnergy.Dtemperature is 2.96630141387408e+001
Property internalEnergy.Dpressure is -1.81207174324797e-004
XXXXXXXXXXXXXXXXXXXXXXXX - test1
Property internalEnergy is -3.48112533568150e+003
Property internalEnergy.Dtemperature is 2.96630141387408e+001
XXXXXXXXXXXXXXXXXXXXXXXXX - test2
Property internalEnergy is -3.47260004151000e+003
Property internalEnergy.Dpressure is -1.81207174324797e-004
XXXXXXXXXXXXXXXXXXXXXXXXX - test3
Property internalEnergy is -3.48112533568150e+003

[jasper]

You are correct that there was a book-keeping problem there (in volume values, when volume is evaluated from the EOS calculation routine in combination with a perturbation, e.g. enthalpy.Dpressure, such as is done by the internal energy routine). Corrected, and the correction is available via CUP. Thank you for reporting this issue.