ansys / pymapdl

Pythonic interface to MAPDL
https://mapdl.docs.pyansys.com
MIT License
428 stars 120 forks source link

Retrieving results with get_float #242

Closed Erwt10 closed 4 years ago

Erwt10 commented 4 years ago

Hi,

I am attempting to improve our flow of retrieving results from pipeline analyses with the use of pyansys. I am running into trouble with retrieving values present in ansys to variables in python. In /POST1 I am able to retrieve the results I want with: (based on issue 157/152)

ansys.run("/POST1")

print("start creating output")
#results at last (=current) timestep via POST1
ansys.run("ESEL,S,ENAME,,59")   #select all elements of pipe59 type
ansys.etable(lab="FX", item="SMISC",comp="1")       #axial force
ansys.etable(lab="SDIR", item="SMISC",comp="13")    #direct axial stress
ansys.etable(lab="SB", item="NMISC",comp="88")       #bending stress
ansys.etable(lab="SEQVL", item="NMISC",comp="87")       #equivalent stress
ansys.etable(lab="MY", item="SMISC",comp="5")       #bending moment
ansys.etable(lab="MZ", item="SMISC",comp="6")       #bending moment

maxelem=100
python_FX = np.zeros(maxelem)
python_SDIR = np.zeros(maxelem)
python_SB = np.zeros(maxelem)
python_SEQVL = np.zeros(maxelem)
python_MY = np.zeros(maxelem)
python_MZ = np.zeros(maxelem)

for i in range(maxelem):
    python_FX[i] = ansys.get_float("ELEM", i+1, "ETAB", "FX")
    python_SDIR[i] = ansys.get_float("ELEM", i+1, "ETAB", "SDIR")
    python_SB[i] = ansys.get_float("ELEM", i+1, "ETAB", "SB")
    python_SEQVL[i] = ansys.get_float("ELEM", i+1, "ETAB", "SEQVL")
    python_MY[i] = ansys.get_float("ELEM", i+1, "ETAB", "MY")
    python_MZ[i] = ansys.get_float("ELEM", i+1, "ETAB", "MZ")
print("POST1 elements completed")

There are also results I want that are retrieved from POST26, however I run into problems with assigning the internal result of ansys to a python variable

#for fatigue extract time-dependant longitudinal stress 
#axial stress per location, via ESOL, which requires element number and node number
#ESOL is only valid in POST26
ansys.set_log_level('INFO')
ansys.post26()   
python_SD_LOC1 = np.zeros(maxelem) 
for i in range(maxelem):    
    print(i)
    ansys.esol(nvar='2',elem=str(i+1),node=str(i+1),item='LS',comp='1') #ls,1 is longitudinal stress at loc 1 (0 deg location)
    ansys.vget(par='SD_LOC1',ir='2',tstrt='1') #store in a variable with a name
    python_SD_LOC1[i] = ansys.get_float("ELEM", i+1, "VGET", "SD_LOC1") # <- I do not understand of use get_float()
    #get_float(entity='', entnum='', item1='', it1num='', item2='', it2num='', **kwargs)

The documentation on this function does not provide enough information to help me further.

Could you elaborate on how to use get_float in this situation, or other means of retrieving data?

example file.zip

akaszynski commented 4 years ago

Nice work on your script!

You're quite right, the documentation isn't great for those methods; I've improved them as well as added an example for them at read_float_parameter

Here's how you would read in those element results from MAPDL from the vget array. You'll need pyansys==0.41.4, which should be out shortly on PyPi.

#for fatigue extract time-dependant longitudinal stress 
#axial stress per location, via ESOL, which requires element number and node number
#ESOL is only valid in POST26
ansys.set_log_level('INFO')
ansys.post26()
python_SD_LOC1 = np.empty(maxelem)
for i in range(maxelem):
    # ls,1 is longitudinal stress at loc 1 (0 deg location)
    ansys.esol(nvar='2', elem=str(i+1), node=str(i+1), item='LS', comp='1')
    ansys.vget(par='SD_LOC1', ir='2', tstrt='1') # store in a variable with a name

    # read the first value in the array
    python_SD_LOC1[i] = ansys.read_float_parameter('SD_LOC1(1)')

However, it's more efficient to read the results directly from the binary file (using pyansys). Here's how you'd go about doing that:

# Read in the element stress for all elements directly from the file
result = ansys.result
estress, enum, enode = result.element_stress(0)

# grab the stress from the first value for each element
python_SD_LOC1_from_rst = np.array([elem[0, 0] for elem in estress[:maxelem]])

Prove they are identical:

>>> np.allclose(python_SD_LOC1_from_rst, python_SD_LOC1)
True

>>> python_SD_LOC1_from_rst
array([-4266.18994141, -4271.24462891, -4276.28662109, -4281.30566406,
       -4286.29101562, -4291.23388672, -4296.12597656, -4300.95898438,
       -4305.72558594, -4310.41992188, -4315.03613281, -4319.56884766,
       ...
       -4421.67919922, -4421.47167969, -4421.26220703, -4421.05126953,
       -4420.83935547, -4420.62695312, -4420.4140625 , -4420.20214844,
       -4419.99072266, -4419.78125   , -4419.57275391, -4419.36669922])
Erwt10 commented 4 years ago

Hi,

Thank you for your response, this helps me make a step forward.

The solution with result.element_stress() is not exactly what I want though. This gives a single stress for a node, but I want to investigate stress at multiple locations around the pipe, which I can extract (in APDL script) with ESOL, i, ENUM(i-1), NNUM(i-1), LS,X where X is 1, 5, ..., 29 for the 8 circumferential positions of PIPE59.

I did some digging and I think the solution could lie in: enum, edata = result.element_solution_data(0, datatype='ENS') this gives 64 values for a PIPE59 element and 192 for a PIPE288 element (while exploring pyansys I am trying to update our analyses to use the more recent element). Am I correct to assume that this corresponds to indices CI and CJ as here ?

This function also gives me some further questions. With datatype= 'EMN' or 'EMS' you can get the 'NMISC' or 'SMISC' components. The retrieved data has 40 or 164 values for these options, respectively. The table only has NMISC components up to 28 and SMISC up to 67. Where are the other values coming from?

akaszynski commented 4 years ago

When you get a chance, please upload your PIPE288 version.

Looking at this with PIPE59 you're correct that "ENS" corresponds to the "LS" ESOL data:

# compare "LS" ESOL results between pyansys and ANSYS
enum, edata, enode = result.element_solution_data(0, datatype='ENS')

# just here to compare getting data using esol vs element_solution_data
def tmp_get_value(comp):
    ansys.esol(nvar='2', elem=enum[0], node=enode[0][0], item='LS', comp=comp)
    ansys.vget(par='SD_LOC1', ir='2', tstrt='1') # store in a variable with a name
    return ansys.read_float_parameter('SD_LOC1(1)')

# show these are identical
ansys.post26()
ansys.set_log_level('ERROR')
from_ansys = [tmp_get_value(i) for i in range(1, 65)]
assert np.allclose(from_ansys, edata[0])

Am I correct to assume that this corresponds to indices CI and CJ as here ?

I can also confirm that there's a correspondence between the table values in the item and sequence numbers and the values of the results from ENS. Here's the table from the 12.1 Element Reference documentation:

image

This function also gives me some further questions. With datatype= 'EMN' or 'EMS' you can get the 'NMISC' or 'SMISC' components. The retrieved data has 40 or 164 values for these options, respectively. The table only has NMISC components up to 28 and SMISC up to 67. Where are the other values coming from?

Without diving into the ANSYS source, I can only speculate that there's extra undocumented data written to the result file that is either used internally within ANSYS to avoid recalculating some results on the fly, or the data being written to file isn't necessary and is a hold-over from some feature.

I appreciate your digging, and it's given me some ideas about how to expose additional data to the user through the binary interface. There's so much data available in these result files, it's quite difficult to provide all of it as a nicely packaged and labeled result. Ideally, the results from an "ENS" query would return a dictionary for each element, but I'm afraid that's out of scope at the moment of this project.

Erwt10 commented 4 years ago

Hi,

It has been a while (sidetracked with other projects), but your response has helped me to complete my task.

I will close this issue as my original question has been resolved. Thank you.