diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 1907a75c74..26077a2d77 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -97,6 +97,7 @@ module MOM_state_initialization use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field use MOM_oda_incupd, only: calc_oda_increments, output_oda_incupd_inc +use MOM_python_embedding, only: runpyF_1_3d implicit none ; private @@ -643,7 +644,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & PF, oda_incupd_CSp, restart_CS, Time) endif - + if (is_root_pe()) call MOM_mesg("MOM_initialize_state: Testing Python embedding ...") + Call runpyF_1_3d('test.py', 'py_test1', 'tv%T', tv%T) + if (is_root_pe()) call MOM_mesg("MOM_initialize_state: Stop here to debug!!"); stop end subroutine MOM_initialize_state !> Reads the layer thicknesses or interface heights from a file. diff --git a/src/initialization/runpyC.C b/src/initialization/runpyC.C new file mode 100644 index 0000000000..40d7777e95 --- /dev/null +++ b/src/initialization/runpyC.C @@ -0,0 +1,154 @@ +//Inspired from https://github.com/wangsl/python-embedding +//#include "runpyC.h" +#include +#include +#include + +#include +#include + +// https://docs.python.org/3/extending/embedding.html +// https://github.com/dusty-nv/jetson-utils/issues/44 +// export PYTHONPATH=.:$PYTHONPATH + +#define FORT(x) x##_ + +static PyObject *py_module = 0; +static PyObject *py_function = 0; + +inline int init_numpy() +{ + if(!PyArray_API) import_array(); + return PyArray_API ? 1 : 0; +} + +static void runpyC_1_1d(const char *py_script_, const char *py_function_, + double *x, const int &n1) +{ + if(!Py_IsInitialized()) Py_Initialize(); + + assert(Py_IsInitialized()); + assert(init_numpy()); + printf("runpyC_1_1d: n1,%d \n", n1); + printf("From runpyC_1_1d: Trying to call %s in %s \n",py_function_,py_script_); + if(!py_module) { + PyObject *py_script = PyUnicode_DecodeFSDefault(py_script_); + assert(py_script); + py_module = PyImport_Import(py_script); +//Interesting: The following assertion fails if +// the py_script has a python bug even outside the function that is called! +// the py_script has import torch + assert(py_module); + Py_DECREF(py_script); + } + if(!py_function) { + assert(py_module); + py_function = PyObject_GetAttrString(py_module, py_function_); + assert(py_function);// This assertation fails if py_function is not found in py_script + assert(PyCallable_Check(py_function)); + } + //Simle tests + //assert(PyObject_CallFunctionObjArgs(py_function, NULL)); //This worked. Printed from the function. + + const npy_intp dim_x [] = { n1 }; + PyObject *x_py = PyArray_SimpleNewFromData(1, dim_x, NPY_DOUBLE, x); + PyObject *pValue; + assert(pValue = PyObject_CallFunctionObjArgs(py_function, x_py, NULL)); + if (pValue != NULL) { + printf("Result of call: %f\n", PyFloat_AsDouble(pValue)); + Py_DECREF(pValue); + } + Py_DECREF(x_py); x_py = 0; + + // PyRun_SimpleString("import sys; sys.stdout.flush()"); + + std::cout.flush(); +} + +static void runpyC_1_3d(const char *py_script_, const char *py_function_, + const char *varname_, + double *x, const int &n1, const int &n2, const int &n3 ) +{ + if(!Py_IsInitialized()) Py_Initialize(); + + assert(Py_IsInitialized()); + assert(init_numpy()); + printf("runpyC_1_3d: varname %s , dims %d,%d,%d \n", varname_, n1,n2,n3); + printf("From runpyC_1_3d: Trying to call %s in %s.py \n",py_function_,py_script_); + if(!py_module) { + PyObject *py_script = PyUnicode_DecodeFSDefault(py_script_); + assert(py_script); + py_module = PyImport_Import(py_script); +//Interesting: The following assertion fails if +// the py_script has a python bug even outside the function that is called! +// the py_script has import torch + assert(py_module && " This assertation fails if the py_script has a python bug.."); + Py_DECREF(py_script); + } + if(!py_function) { + assert(py_module); + py_function = PyObject_GetAttrString(py_module, py_function_); + assert(py_function && " This assertation fails if the function is not found in the script."); //This assertation fails if py_function is not found in py_script + assert(PyCallable_Check(py_function)); + } + //Simle tests + //assert(PyObject_CallFunctionObjArgs(py_function, NULL)); //This worked. Printed from the function. + + const npy_intp dim_x [] = { n1,n2,n3 }; + PyObject *x_py = PyArray_SimpleNewFromData(3, dim_x, NPY_DOUBLE, x); + PyObject *pValue; + assert(pValue = PyObject_CallFunctionObjArgs(py_function, x_py, NULL)); + if (pValue != NULL) { + printf("Result of call: %f\n", PyFloat_AsDouble(pValue)); + Py_DECREF(pValue); + } + Py_DECREF(x_py); x_py = 0; + + // PyRun_SimpleString("import sys; sys.stdout.flush()"); + + std::cout.flush(); +} + +static void py_finalize() +{ + if(py_function) { Py_DECREF(py_function); py_function = 0; } + if(py_module) { Py_DECREF(py_module); py_module = 0; } + if(Py_IsInitialized()) assert(!Py_FinalizeEx()); + std::cout.flush(); +} + +// Fortran interface: NOTE: the name of function should be al lowercase! +extern "C" void FORT(runpyc_1array3d) (const char *py_script, const int &len_py_script, + const char *py_function, const int &len_py_function, + const char *varname, const int &len_varname, + double *x, const int &n1, const int &n2, const int &n3 ) + +{ + char *py_script_ = new char [len_py_script+1]; + assert(py_script); + memcpy(py_script_, py_script, len_py_script*sizeof(char)); + py_script_[len_py_script] = '\0'; + + char *py_function_ = new char [len_py_function+1]; + assert(py_function); + memcpy(py_function_, py_function, len_py_function*sizeof(char)); + py_function_[len_py_function] = '\0'; + + char *varname_ = new char [len_varname+1]; + assert(varname); + memcpy(varname_, varname, len_varname*sizeof(char)); + varname_[len_varname] = '\0'; + + runpyC_1_3d(py_script_, py_function_, varname_, x, n1,n2,n3); + + if(py_script_) { delete [] py_script_; py_script_ = 0; } + if(py_function_) { delete [] py_function_; py_function_ = 0; } + if(varname_) { delete [] varname_; varname_ = 0; } +} + +// Fortran interface: PyFinalize +extern "C" void FORT(pyfinalize)() +{ + py_finalize(); +} + diff --git a/src/initialization/runpyF.F90 b/src/initialization/runpyF.F90 new file mode 100644 index 0000000000..8e1c3883cb --- /dev/null +++ b/src/initialization/runpyF.F90 @@ -0,0 +1,26 @@ +!Inspired by https://github.com/wangsl/python-embedding +! +module MOM_python_embedding +implicit none +public runPyF_1_3d +contains + +subroutine runpyF_1_3d(pyScript, pyFunction, varname, x) + character(len=*), intent(in) :: pyScript, pyFunction,varname + real, dimension(:,:,:), intent(inout) :: x + integer :: n1,n2,n3, len_pyScript, len_pyFunction , len_varname + + len_pyScript = Len_Trim(pyScript) + if(pyScript(len_pyScript-2:len_pyScript) .eq. '.py') len_pyScript = len_pyScript-3 + len_pyFunction = Len_Trim(pyFunction) + len_varname = Len_Trim(varname) + + n1 = size(x,1) + n2 = size(x,2) + n3 = size(x,3) + call runpyc_1array3d(pyScript, len_pyScript, pyFunction, len_pyFunction,& + varname, len_varname, x, n1,n2,n3) + +end subroutine runpyF_1_3d + +end module MOM_python_embedding diff --git a/src/initialization/test.py b/src/initialization/test.py new file mode 100644 index 0000000000..e69f0acf81 --- /dev/null +++ b/src/initialization/test.py @@ -0,0 +1,73 @@ +#!/bin/env python +#Inspired from https://github.com/wangsl/python-embedding +import sys +import numpy +import matplotlib.pyplot as plt +#A test function with one argument +def py_test1(x) : + print(' **** From py_test1 ****') + print("py_test1: Shape of the input array x : ", x.shape) + print("py_test1: x[4,4,0] x[4,4,-1] : ", x[4,4,0],x[4,4,-1]) + + plt.plot(x[4,4,:]);plt.show(); + + a=numpy.max(x) + print("py_test1: Returning the numpy.max(x): ", a) + + sys.stdout.flush() + return a + +#A test function with two arguments +def py_test2(x,y) : + print(' **** From py_test2 ****',x,y) + sys.stdout.flush() + a=sum(x) + b=sum(y) + return a,b + +#A test function with no arguments +def py_test0(): + print(' **** From py_test0 ****') + sys.stdout.flush() + return + +from math import sin +#import torch +#import numpy #This gives all kinds of library errors at runtime + +def my_test_torch(x, y) : + print(' From Python test: {}'.format(x.size)) + for i in range(x.size) : + x[i] += 1.0 + + a = torch.from_numpy(x) + print(a) + + b = torch.from_numpy(x.astype(numpy.float32)) + print(b) + + if torch.cuda.is_available() : + b_dev = b.cuda() + print(b_dev) + + for i in range(y.size) : + y[i] = 2*numpy.float64(b[i]) + 0.1 + + print(y) + + sys.stdout.flush() + + +if __name__ == '__main__' : + import numpy as np + + x = np.arange(10, dtype=numpy.float64) + print(x) + + y = np.arange(10, dtype=numpy.float64) + print(y) + my_test(x, y) + + print(x) + +