RadialSolver solution gives unallocated memory for Love number results

[jrenaud90]

Users have reported that successful runs of RadialSolver.radial_solver give unallocated values when accessing the solutions.k, .h, .l etc.

# Problem
def run():
    solution = radial_solver(*input)
    k = solution.k
    return k

k2 = run()
print(k2)  # Random value. Usually something super small like ~1e-300

# Workaround
import numpy as np
def run():
    solution = radial_solver(*input)
    k = solution.k[0]
    # or 
    k = np.copy(solution.k)
    return k
    # or return the solution instance: `return solution`

k2 = run()
print(k2)  # Correct value!

This has occurred with, at least, TidalPy v0.5.4 but as discussed below, it is not a TidalPy problem per-say.

The problem is that the RadialSolverSolution class performs a heap allocation to memory in order to store the Love number results. It is responsible for freeing this memory when the class instance is destroyed. This particular problem comes up when the solution is destroyed too early.

To Reproduce The common situation where this problem arises is when an outer python function calls radial solver and gets a solution but does not return the solution. This can be reproduced with pure Cython and no use of TidalPy:

To reproduce make a cython class and functions like...

## Cython
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free

import numpy as np
cimport numpy as np

cdef class Solution:

    cdef size_t n
    cdef double complex* result_ptr
    cdef double complex[::1] result_view

    def __init__(self, size_t n):

        self.n = n

        self.result_ptr = <double complex*> PyMem_Malloc(self.n * sizeof(double complex))
        if not (self.result_ptr is NULL):
            self.result_view = <double complex[:self.n]>self.result_ptr

    @property
    def result(self):
        return np.ascontiguousarray(self.result_view, dtype=np.complex128)

    def __dealloc__(self):
        if not (self.result_ptr is NULL):
            PyMem_Free(self.result_ptr)

cdef Solution cy_run():

    cdef size_t n = 3
    cdef size_t i
    cdef Solution solution = Solution(n)

    for i in range(n):
        solution.result_ptr[i] = i * (10. + 1.j)

    return solution

def run():
    return cy_run()

Call this function from another python wrapping function...

## Python
import numpy as np
def run_outer():

    solution = run()
    k = solution.result

    return k

print(run_outer())

You will see unallocated values are printed.

Workaround We can get the correct values by grabbing a specific value of the array, or by making a copy of the array, then the correct values are printed...

import numpy as np
def run_outer():

    result = run()
    k = np.copy(result.result)
    # or
    k = result.result[0]

    return k

print(run_outer())

Alternatively, and probably the better solution, is to return the solution instance itself so that it is no longer tied to the scope of run_outer

Future Fix It is hard to know the solution to this problem that is not one of the workarounds mentioned above. We want RadialSolverSolution to own the heap memory otherwise its owner may never free it properly. The downside is a situation like this where it is deallocated too early.

Certain open to suggestions but for now this probably just requires updates to the documentation.

[jrenaud90]

Okay found a solution to this but not sure if I will implement. If in the RadialSolverSolution method we change...

@property
def k(self):
    return np.ascontiguousarray(self.complex_love_view[0::3])
    # Replace with  
    return np.copy(np.ascontiguousarray(self.complex_love_view[0::3]))

Then the problem is fixed. But that extra copy doubles the overhead of calling that property attribute.