jranalli/pyromat_plot_demo.py

## pyromat_plot_demo.py
import pyromat as pm
import numpy as np

import matplotlib.pyplot as plt


def ismultiphase(subst):
    """Test whether the PYroMat substance instance is a multi-phase model
    :param subst: A PYroMat substance instance
    :return True/False:
"""
    # Explicitly list the supported classes
    # This was my favorite of a few candidate algorithms.  It is quick,
    # it ensures that the instance will have the correct property
    # methods (which is really what you want to test for), and it will
    # be immune from creating strange new substance models in the future.
    testfor = ['mp1']
    for thisclass in testfor:
        if isinstance(subst, pm.reg.registry[thisclass]):
            return True
    return False
    # There are a few other candidate algorithms we could consider:
    # 1) We could test the substance's collection from its id string
    #    This has advantages if we ever add ideal liquid, solid, or
    #    other substance collections.
    # 2) We could test for a complete list of mandatory property methods
    #    This would ensure compatibility with the algorithm, but it
    #    would be slower.
    # 3) We could test for a single candidate property method (like ps)
    #    This isn't bad, but it's prone to problems if we aren't careful
    #    down the line.  However unlikely, it's possible that we'll make
    #    a new substance model with ps() that has a different meaning.


def get_practical_limits(subst):
    """Return practical temperature and pressure boundaries for a substance
    Tmin,pmin,Tmax,pmax = get_practical_limits( subst )

Temeprature limits are always set based on the substance's Tlim() method,
but ideal gas substance models have no explicit pressure limits, and
most multi-phase minimum pressures are zero.  If the substance is a
multiphase instance with zero minimum pressure, 10% of the triple point
is used.  Ideal gas models use the pressure at T=Tmin,d=0.01 and T=Tmax,
d=1000.
"""
    if ismultiphase(subst):
        pmin, pmax = subst.plim()
        Tt, pt = subst.triple()
        if pmin == 0:
            pmin = 0.1 * pt
        Tmin, Tmax = subst.Tlim()
    else:
        Tmin, Tmax = subst.Tlim()
        pmin, pmax = np.array(
            [subst.p(T=Tmin, d=0.01), subst.p(T=Tmax, d=1000)]).flatten()

    peps_hi, peps_lo = 0.01 * pmax, pmin
    Teps = 0.001 * (Tmax - Tmin)

    return Tmin + Teps, pmin + peps_lo, Tmax - Teps, pmax - peps_hi


def get_default_lines(subst, prop):
    """
    Get a default set of isolines for a given property.
    :param subst: A pyromat substance object
    :param prop: A string representing the requested property
    :return: vals - a np array of suitable default values.
    """
    if not hasattr(subst, prop):
        raise pm.utility.PMParamError(f"{subst} has no such property: {prop}")

    vals = None

    multiphase = hasattr(subst, 'ps')
    Tmin, pmin, Tmax, pmax = get_practical_limits(subst)
    if multiphase:
        Tt, pt = subst.triple()

    if prop == 'p':
        vals = np.logspace(np.log10(pmin), np.log10(pmax), 10)
    elif prop == 'T':
        vals = np.linspace(Tmin, Tmax, 10)
    elif prop == 'x':
        vals = np.linspace(0.1, 0.9, 9)
    elif prop in ['v', 'd', 'h', 'e', 's']:
        pfn = getattr(subst, prop)
        try:  # Finding the low value can be flaky for some substances
            propmin = pfn(T=Tmin, p=pmax)
        except pm.utility.PMParamError:
            if multiphase:
                pfns = getattr(subst, prop + 's')
                propmin = pfns(T=Tt)[0]
            else:
                propmin = pfn(T=Tmin, p=pmin)
        propmax = pfn(T=Tmax, p=pmin)
        if prop == 'v':
            vals = np.logspace(np.log10(propmin), np.log10(propmax), 10)
        elif prop == 'd':  # Inverse of v means max and min flip
            vals = np.logspace(np.log10(propmax), np.log10(propmin), 10)
        else:
            vals = np.linspace(propmin, propmax, 10)
    else:
        raise pm.utility.PMParamError(f'Default Lines Undefined for {prop}.')

    return vals

def compute_iso_line(subst, n=25, scaling='linear', **kwargs):
    """
    Compute a constant line for a given property at a given value
    :param subst: a pyromat substance object
    :param n: The number of points to compute to define the line
    :param scaling: Should point spacing be 'linear' or 'log'
    :param kwargs: A property specified by name. If 'default' is specified in
                    kwargs, the value of the prop will be ignored and a set
                    of default lines for that prop will be computed (see
                    get_default_lines()).
    :return: A dict containing arrays of properties. If 'default' flag is set
                the response will be an array of dicts representing all the
                individual lines.
    """

    # Keep track of whether the users want the defaults
    default_mode = False
    if 'default' in kwargs:
        kwargs.pop('default')  # remove from prop list
        default_mode = True

    # Test the legality of inputs
    if len(kwargs) != 1:
        raise pm.utility.PMParamError("Specify exactly one property "
                                      "for an isoline")

    # We now definitely have only a single property, so pull its name
    prop = list(kwargs.keys())[0]

    # If the user wants multiple lines, we'll behave differently
    multiline = None
    if default_mode:
        # In default_mode, we'll get the default lines for the substance
        multiline = get_default_lines(subst, prop)
    elif hasattr(kwargs[prop], "__iter__") and np.size(kwargs[prop]) > 1:
        # The user might have also submitted an array to ask for multiple lines
        multiline = kwargs[prop]

    # If we're looking for multiples, compute each by recursion
    if multiline is not None:
        lines = []
        for val in multiline:
            arg = {prop: val}  # Build an argument
            try:
                lines.append(compute_iso_line(subst, n, scaling, **arg))
            except pm.utility.PMParamError:
                # This may error if stuff is out of bounds, just skip that line
                pass
        return lines

    # We were asked for a single line, so begin computations for that line

    # Set P&T limits, including special cases for multiphase
    multiphase = hasattr(subst, 'Ts')
    Tmin, pmin, Tmax, pmax = get_practical_limits(subst)
    if multiphase:
        Tc, pc = subst.critical()
        Tt, pt = subst.triple()

    # The props for which we will plot against a T list
    if prop in ['p', 'd', 'v', 's', 'x']:

        # If quality, we stop at the crit pt
        if 'x' in kwargs:
            if multiphase:
                Tmax = Tc
            else:
                raise pm.utility.PMParamError('x cannot be computed for non-'
                                              'multiphase substances.')

        line_T = np.linspace(Tmin, Tmax, n).flatten()

        # We can insert the phase change points
        if multiphase and 'p' in kwargs and pc > kwargs['p'] > pt:
            Tsat = subst.Ts(p=kwargs['p'])
            i_insert = np.argmax(line_T > Tsat)

            line_T = np.insert(line_T, i_insert,
                               np.array([Tsat, Tsat]).flatten())
            x = -np.ones_like(line_T)
            x[line_T == Tsat] = np.array([0, 1])
            kwargs['x'] = x

        kwargs['T'] = line_T

    elif prop in ['h', 'e']:
        try:  # Finding the low value can be flaky for some substances
            dmax = subst.d(T=Tmin, p=pmax)
        except pm.utility.PMParamError:
            if multiphase:
                dmax = subst.ds(T=Tt)[0]
            else:
                dmax = subst.d(T=Tmin, p=pmin)
        dmin = subst.d(T=Tmax, p=pmin)
        line_d = np.logspace(np.log10(dmin), np.log10(dmax), n).flatten()
        # line_d = np.linspace(dmin, dmax, n)
        kwargs['d'] = line_d

    elif prop in ['T', 'h', 'e']:
        # ph & pe are going to be really slow, but what's better?

        line_p = np.logspace(np.log10(pmin), np.log10(pmax), n).flatten()

        # We can insert the phase change points
        if multiphase and 'T' in kwargs and Tc > kwargs['T'] > Tt:
            psat = subst.ps(T=kwargs['T'])
            i_insert = np.argmax(line_p > psat)

            line_p = np.insert(line_p, i_insert,
                               np.array([psat, psat]).flatten())
            x = -np.ones_like(line_p)
            x[line_p == psat] = np.array([1, 0])
            kwargs['x'] = x

        kwargs['p'] = line_p

    else:  # Should never arrive here without error
        raise pm.utility.PMParamError('property invalid')

    states = subst.state(**kwargs)

    return states

def compute_steam_dome(subst):
    Tc, pc, dc = subst.critical(density=True)
    Tt, pt = subst.triple()
    # Use an epsilon, we'll manually add the critical point later
    ep = (Tc - Tt) * .01
    Ts = np.linspace(Tt + ep, Tc - ep, 31)

    data = {}

    # OK, we've got Ts - go calculate the state
    dsL, dsV = subst.ds(T=Ts)
    data['liquid'] = subst.state(T=Ts, d=dsL)
    data['vapor'] = subst.state(T=Ts, d=dsV)
    # Throw away the liquid pressure - it is not numerically correct.
    data['liquid']['p'] = data['vapor']['p']


    crit_state = subst.state(p=pc, d=dc)
    for prop in crit_state:
        data['liquid'][prop] = \
            np.append(data['liquid'][prop], crit_state[prop])
        data['vapor'][prop] = \
            np.append(data['vapor'][prop], crit_state[prop])

    return data


if __name__ == "__main__":
    water = pm.get("mp.H2O")
    sd_data = compute_steam_dome(water)

    plt.figure()
    plt.xlabel('s (kJ/kg-K)')
    plt.ylabel('T (K)')

    # Steam Dome
    for line in ['liquid', 'vapor']:
        plt.plot(sd_data[line]['s'], sd_data[line]['T'], 'k-')

    # P isolines
    for line in compute_iso_line(water, p=0, n=25, default=True):
        plt.plot(line['s'], line['T'], 'r-')

    # v isolines
    for line in compute_iso_line(water, v=0, n=25, default=True):
        plt.plot(line['s'], line['T'], 'b-')

    # h isolines
    for line in compute_iso_line(water, h=0, n=25, default=True):
        plt.plot(line['s'], line['T'], 'g-')

    plt.show()
	import pyromat as pm
	import numpy as np

	import matplotlib.pyplot as plt


	def ismultiphase(subst):
	"""Test whether the PYroMat substance instance is a multi-phase model
	:param subst: A PYroMat substance instance
	:return True/False:
	"""
	# Explicitly list the supported classes
	# This was my favorite of a few candidate algorithms. It is quick,
	# it ensures that the instance will have the correct property
	# methods (which is really what you want to test for), and it will
	# be immune from creating strange new substance models in the future.
	testfor = ['mp1']
	for thisclass in testfor:
	if isinstance(subst, pm.reg.registry[thisclass]):
	return True
	return False
	# There are a few other candidate algorithms we could consider:
	# 1) We could test the substance's collection from its id string
	# This has advantages if we ever add ideal liquid, solid, or
	# other substance collections.
	# 2) We could test for a complete list of mandatory property methods
	# This would ensure compatibility with the algorithm, but it
	# would be slower.
	# 3) We could test for a single candidate property method (like ps)
	# This isn't bad, but it's prone to problems if we aren't careful
	# down the line. However unlikely, it's possible that we'll make
	# a new substance model with ps() that has a different meaning.



	def get_practical_limits(subst):
	"""Return practical temperature and pressure boundaries for a substance
	Tmin,pmin,Tmax,pmax = get_practical_limits( subst )

	Temeprature limits are always set based on the substance's Tlim() method,
	but ideal gas substance models have no explicit pressure limits, and
	most multi-phase minimum pressures are zero. If the substance is a
	multiphase instance with zero minimum pressure, 10% of the triple point
	is used. Ideal gas models use the pressure at T=Tmin,d=0.01 and T=Tmax,
	d=1000.
	"""
	if ismultiphase(subst):
	pmin, pmax = subst.plim()
	Tt, pt = subst.triple()
	if pmin == 0:
	pmin = 0.1 * pt
	Tmin, Tmax = subst.Tlim()
	else:
	Tmin, Tmax = subst.Tlim()
	pmin, pmax = np.array(
	[subst.p(T=Tmin, d=0.01), subst.p(T=Tmax, d=1000)]).flatten()

	peps_hi, peps_lo = 0.01 * pmax, pmin
	Teps = 0.001 * (Tmax - Tmin)

	return Tmin + Teps, pmin + peps_lo, Tmax - Teps, pmax - peps_hi


	def get_default_lines(subst, prop):
	"""
	Get a default set of isolines for a given property.
	:param subst: A pyromat substance object
	:param prop: A string representing the requested property
	:return: vals - a np array of suitable default values.
	"""
	if not hasattr(subst, prop):
	raise pm.utility.PMParamError(f"{subst} has no such property: {prop}")

	vals = None

	multiphase = hasattr(subst, 'ps')
	Tmin, pmin, Tmax, pmax = get_practical_limits(subst)
	if multiphase:
	Tt, pt = subst.triple()

	if prop == 'p':
	vals = np.logspace(np.log10(pmin), np.log10(pmax), 10)
	elif prop == 'T':
	vals = np.linspace(Tmin, Tmax, 10)
	elif prop == 'x':
	vals = np.linspace(0.1, 0.9, 9)
	elif prop in ['v', 'd', 'h', 'e', 's']:
	pfn = getattr(subst, prop)
	try: # Finding the low value can be flaky for some substances
	propmin = pfn(T=Tmin, p=pmax)
	except pm.utility.PMParamError:
	if multiphase:
	pfns = getattr(subst, prop + 's')
	propmin = pfns(T=Tt)[0]
	else:
	propmin = pfn(T=Tmin, p=pmin)
	propmax = pfn(T=Tmax, p=pmin)
	if prop == 'v':
	vals = np.logspace(np.log10(propmin), np.log10(propmax), 10)
	elif prop == 'd': # Inverse of v means max and min flip
	vals = np.logspace(np.log10(propmax), np.log10(propmin), 10)
	else:
	vals = np.linspace(propmin, propmax, 10)
	else:
	raise pm.utility.PMParamError(f'Default Lines Undefined for {prop}.')

	return vals

	def compute_iso_line(subst, n=25, scaling='linear', **kwargs):
	"""
	Compute a constant line for a given property at a given value
	:param subst: a pyromat substance object
	:param n: The number of points to compute to define the line
	:param scaling: Should point spacing be 'linear' or 'log'
	:param kwargs: A property specified by name. If 'default' is specified in
	kwargs, the value of the prop will be ignored and a set
	of default lines for that prop will be computed (see
	get_default_lines()).
	:return: A dict containing arrays of properties. If 'default' flag is set
	the response will be an array of dicts representing all the
	individual lines.
	"""

	# Keep track of whether the users want the defaults
	default_mode = False
	if 'default' in kwargs:
	kwargs.pop('default') # remove from prop list
	default_mode = True

	# Test the legality of inputs
	if len(kwargs) != 1:
	raise pm.utility.PMParamError("Specify exactly one property "
	"for an isoline")

	# We now definitely have only a single property, so pull its name
	prop = list(kwargs.keys())[0]

	# If the user wants multiple lines, we'll behave differently
	multiline = None
	if default_mode:
	# In default_mode, we'll get the default lines for the substance
	multiline = get_default_lines(subst, prop)
	elif hasattr(kwargs[prop], "__iter__") and np.size(kwargs[prop]) > 1:
	# The user might have also submitted an array to ask for multiple lines
	multiline = kwargs[prop]

	# If we're looking for multiples, compute each by recursion
	if multiline is not None:
	lines = []
	for val in multiline:
	arg = {prop: val} # Build an argument
	try:
	lines.append(compute_iso_line(subst, n, scaling, **arg))
	except pm.utility.PMParamError:
	# This may error if stuff is out of bounds, just skip that line
	pass
	return lines

	# We were asked for a single line, so begin computations for that line

	# Set P&T limits, including special cases for multiphase
	multiphase = hasattr(subst, 'Ts')
	Tmin, pmin, Tmax, pmax = get_practical_limits(subst)
	if multiphase:
	Tc, pc = subst.critical()
	Tt, pt = subst.triple()

	# The props for which we will plot against a T list
	if prop in ['p', 'd', 'v', 's', 'x']:

	# If quality, we stop at the crit pt
	if 'x' in kwargs:
	if multiphase:
	Tmax = Tc
	else:
	raise pm.utility.PMParamError('x cannot be computed for non-'
	'multiphase substances.')

	line_T = np.linspace(Tmin, Tmax, n).flatten()

	# We can insert the phase change points
	if multiphase and 'p' in kwargs and pc > kwargs['p'] > pt:
	Tsat = subst.Ts(p=kwargs['p'])
	i_insert = np.argmax(line_T > Tsat)

	line_T = np.insert(line_T, i_insert,
	np.array([Tsat, Tsat]).flatten())
	x = -np.ones_like(line_T)
	x[line_T == Tsat] = np.array([0, 1])
	kwargs['x'] = x

	kwargs['T'] = line_T

	elif prop in ['h', 'e']:
	try: # Finding the low value can be flaky for some substances
	dmax = subst.d(T=Tmin, p=pmax)
	except pm.utility.PMParamError:
	if multiphase:
	dmax = subst.ds(T=Tt)[0]
	else:
	dmax = subst.d(T=Tmin, p=pmin)
	dmin = subst.d(T=Tmax, p=pmin)
	line_d = np.logspace(np.log10(dmin), np.log10(dmax), n).flatten()
	# line_d = np.linspace(dmin, dmax, n)
	kwargs['d'] = line_d

	elif prop in ['T', 'h', 'e']:
	# ph & pe are going to be really slow, but what's better?

	line_p = np.logspace(np.log10(pmin), np.log10(pmax), n).flatten()

	# We can insert the phase change points
	if multiphase and 'T' in kwargs and Tc > kwargs['T'] > Tt:
	psat = subst.ps(T=kwargs['T'])
	i_insert = np.argmax(line_p > psat)

	line_p = np.insert(line_p, i_insert,
	np.array([psat, psat]).flatten())
	x = -np.ones_like(line_p)
	x[line_p == psat] = np.array([1, 0])
	kwargs['x'] = x

	kwargs['p'] = line_p

	else: # Should never arrive here without error
	raise pm.utility.PMParamError('property invalid')

	states = subst.state(**kwargs)

	return states

	def compute_steam_dome(subst):
	Tc, pc, dc = subst.critical(density=True)
	Tt, pt = subst.triple()
	# Use an epsilon, we'll manually add the critical point later
	ep = (Tc - Tt) * .01
	Ts = np.linspace(Tt + ep, Tc - ep, 31)

	data = {}

	# OK, we've got Ts - go calculate the state
	dsL, dsV = subst.ds(T=Ts)
	data['liquid'] = subst.state(T=Ts, d=dsL)
	data['vapor'] = subst.state(T=Ts, d=dsV)
	# Throw away the liquid pressure - it is not numerically correct.
	data['liquid']['p'] = data['vapor']['p']


	crit_state = subst.state(p=pc, d=dc)
	for prop in crit_state:
	data['liquid'][prop] = \
	np.append(data['liquid'][prop], crit_state[prop])
	data['vapor'][prop] = \
	np.append(data['vapor'][prop], crit_state[prop])

	return data


	if __name__ == "__main__":
	water = pm.get("mp.H2O")
	sd_data = compute_steam_dome(water)

	plt.figure()
	plt.xlabel('s (kJ/kg-K)')
	plt.ylabel('T (K)')

	# Steam Dome
	for line in ['liquid', 'vapor']:
	plt.plot(sd_data[line]['s'], sd_data[line]['T'], 'k-')

	# P isolines
	for line in compute_iso_line(water, p=0, n=25, default=True):
	plt.plot(line['s'], line['T'], 'r-')

	# v isolines
	for line in compute_iso_line(water, v=0, n=25, default=True):
	plt.plot(line['s'], line['T'], 'b-')

	# h isolines
	for line in compute_iso_line(water, h=0, n=25, default=True):
	plt.plot(line['s'], line['T'], 'g-')

	plt.show()