From 46e5797c7e1edf52371dbb1946735ea8de73f4c5 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Wed, 6 Oct 2021 18:25:39 -0400 Subject: [PATCH 1/4] An effort to embed a Python script in MOM6 via C++ - Credits for the general procedure goes to https://github.com/wangsl/python-embedding - Also of interest is https://docs.python.org/3/extending/embedding.html https://www.noahbrenowitz.com/post/calling-fortran-from-python/ --- .../MOM_state_initialization.F90 | 12 ++ src/initialization/pyext.C | 137 ++++++++++++++++++ src/initialization/runpy.F | 48 ++++++ src/initialization/test.py | 55 +++++++ 4 files changed, 252 insertions(+) create mode 100644 src/initialization/pyext.C create mode 100644 src/initialization/runpy.F create mode 100644 src/initialization/test.py diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 1907a75c74..865aff53ce 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -184,6 +184,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & # include "version_variable.h" integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB + Integer Nx, Ny + Parameter(Nx = 15) + Parameter(Ny = 15) + Real*8 X(Nx), Y(Ny) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -643,7 +647,15 @@ 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 ...") + Do I = 1, Nx + X(I) = DBLE(I) + Y(I) = 0.0D0 + EndDo + Call RunPy('test.py', 'my_test', X, Nx, Y, Ny) !This worked + !Call RunPy_3darray('test.py', 'my_test', h, is,ie,js,je,1,nz) + if (is_root_pe()) call MOM_mesg("MOM_initialize_state: Stop to debug!!"); stop end subroutine MOM_initialize_state !> Reads the layer thicknesses or interface heights from a file. diff --git a/src/initialization/pyext.C b/src/initialization/pyext.C new file mode 100644 index 0000000000..25e1ec32f7 --- /dev/null +++ b/src/initialization/pyext.C @@ -0,0 +1,137 @@ +#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 py_test(const char *py_script_, const char *py_function_, + double *x, const int &nx, + double *y, const int &ny) +{ + if(!Py_IsInitialized()) Py_Initialize(); + + assert(Py_IsInitialized()); + assert(init_numpy()); + + 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. + //int x_in = 1; + //int y_in = 1; + //assert(PyObject_CallFunctionObjArgs(py_function, x_in, y_in, NULL)); //Segmentation fault - invalid memory reference + //assert(PyObject_CallFunctionObjArgs(py_function, x_in, NULL)); //Segmentation fault - invalid memory reference + //assert(PyObject_CallFunctionObjArgs(py_function, &x_in, NULL)); //Works but cannot access value. Ahaa, input needs to be a poniter + //assert(PyObject_CallFunctionObjArgs(py_function, dims, NULL)); //Works but cannot access value. invalid memory reference + /* The following worked fine. + long dims[2]; + dims[0] = 4; dims[1] = 4; + PyObject *pinValue, *pArgs,*pValue ; + pArgs = PyTuple_New(2); + for (int i = 0; i < 2; ++i) { + pinValue = PyLong_FromLong(dims[i]); + PyTuple_SetItem(pArgs, i, pinValue); + } + pValue = PyObject_CallFunctionObjArgs(py_function, pArgs, NULL); + if (pValue != NULL) { + printf("Result of call: %ld\n", PyLong_AsLong(pValue)); + Py_DECREF(pValue); + } + + double ddims[2]; + ddims[0] = 4.5; ddims[1] = 5.5; + pArgs = PyTuple_New(2); + for (int i = 0; i < 2; ++i) { + pinValue = PyFloat_FromDouble(ddims[i]); + PyTuple_SetItem(pArgs, i, pinValue); + } + pValue = PyObject_CallFunctionObjArgs(py_function, pArgs, NULL); + if (pValue != NULL) { + printf("Result of call: %f\n", PyFloat_AsDouble(pValue)); + Py_DECREF(pValue); + } + */ + const npy_intp dim_x [] = { nx }; + PyObject *x_py = PyArray_SimpleNewFromData(1, dim_x, NPY_DOUBLE, x); + const npy_intp dim_y [] = { ny }; + PyObject *y_py = PyArray_SimpleNewFromData(1, dim_y, NPY_DOUBLE, y); + 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; + Py_DECREF(y_py); y_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: PyTest +extern "C" void FORT(pytest)(const char *py_script, const int &len_py_script, + const char *py_function, const int &len_py_function, + double *x, const int &nx, + double *y, const int &ny) +{ + 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'; + + py_test(py_script_, py_function_, x, nx, y, ny); + + if(py_script_) { delete [] py_script_; py_script_ = 0; } + if(py_function_) { delete [] py_function_; py_function_ = 0; } +} + +// Fortran interface: PyFinalize +extern "C" void FORT(pyfinalize)() +{ + py_finalize(); +} + diff --git a/src/initialization/runpy.F b/src/initialization/runpy.F new file mode 100644 index 0000000000..d543f92fb6 --- /dev/null +++ b/src/initialization/runpy.F @@ -0,0 +1,48 @@ + Subroutine RunPy(PyScript, PyFunction, X, Nx, Y, Ny) + Implicit None + Character*(*) PyScript, PyFunction + Real*8 X(*), Y(*) + Integer Nx, Ny + + Integer I, Len, LinEnd + External LinEnd + Integer LenPyScript, LenPyFunction + + LenPyScript = LinEnd(PyScript) + If(PyScript(LenPyScript-2:LenPyScript) .eq. '.py') + $ LenPyScript = LenPyScript-3 + + LenPyFunction = LinEnd(PyFunction) + + Call Flush(6) + + Call PyTest(PyScript, LenPyScript, + $ PyFunction, LenPyFunction, + $ X, Nx, Y, Ny) + + Return + End + + integer function linend(cline) + implicit real*8(a-h,o-z) +c +c function which returns the length of a character string, +c excluding trailing blanks, tabs and nulls. +c + character*(*) cline, blank + character*1 null, tab + parameter (blank=' ') +c + null = char(0) + tab = char(9) + do 5 i = len(cline), 1, -1 + if(cline(i:i).ne.blank.and. + $ cline(i:i).ne.null.and. + $ cline(i:i).ne.tab) goto 10 + 5 continue + linend = 0 + return + 10 linend = i + return + end + diff --git a/src/initialization/test.py b/src/initialization/test.py new file mode 100644 index 0000000000..b6277a96a5 --- /dev/null +++ b/src/initialization/test.py @@ -0,0 +1,55 @@ +#!/bin/env python + +import sys +from math import sin +#import torch +#import numpy #This gives all kinds of library errors at runtime + +def my_test0(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() + +def my_test(x) : + print(' **** From my_test1 ****',x) + sys.stdout.flush() + a=sum(x) + return a + +def my_test2(x,y) : + print(' **** From my_test2 ****',x,y) + sys.stdout.flush() + a=sum(x) + b=sum(y) + return a,b + +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) + + From 72b8a0d83923d01adeffb2e053331b3ad6d4d581 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Wed, 6 Oct 2021 18:55:38 -0400 Subject: [PATCH 2/4] An effort to embed a Python script in MOM6 via C++ - Credits for the general procedure goes to https://github.com/wangsl/python-embedding - Also of interest is https://docs.python.org/3/extending/embedding.html https://www.noahbrenowitz.com/post/calling-fortran-from-python/ --- src/initialization/pyext.C | 1 + src/initialization/runpy.F | 2 ++ src/initialization/test.py | 2 +- 3 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/initialization/pyext.C b/src/initialization/pyext.C index 25e1ec32f7..2b08c26a38 100644 --- a/src/initialization/pyext.C +++ b/src/initialization/pyext.C @@ -1,3 +1,4 @@ +//Initially copied from https://github.com/wangsl/python-embedding #include #include #include diff --git a/src/initialization/runpy.F b/src/initialization/runpy.F index d543f92fb6..405763467b 100644 --- a/src/initialization/runpy.F +++ b/src/initialization/runpy.F @@ -1,3 +1,5 @@ +!Initially copied from https://github.com/wangsl/python-embedding +! Subroutine RunPy(PyScript, PyFunction, X, Nx, Y, Ny) Implicit None Character*(*) PyScript, PyFunction diff --git a/src/initialization/test.py b/src/initialization/test.py index b6277a96a5..8b7a263bc6 100644 --- a/src/initialization/test.py +++ b/src/initialization/test.py @@ -1,5 +1,5 @@ #!/bin/env python - +#Initially copied from https://github.com/wangsl/python-embedding import sys from math import sin #import torch From c7d65e96413d2d2f86b206d999c383c822f2e734 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 7 Oct 2021 17:57:09 -0400 Subject: [PATCH 3/4] Operate with a Python tool on a given MOM6 array - This update shows an example to use Python tool in model runtime to operate on a MOM6 array and calculate/plot something. - This works for a demo in single column model that can be obtained from https://github.com/nikizadehgfdl/MOM6_OBGC_examples --- .../MOM_state_initialization.F90 | 7 +- src/initialization/pyext.C | 138 ---------------- src/initialization/runpy.F | 50 ------ src/initialization/runpyC.C | 154 ++++++++++++++++++ src/initialization/runpyF.F90 | 26 +++ src/initialization/test.py | 46 ++++-- 6 files changed, 216 insertions(+), 205 deletions(-) delete mode 100644 src/initialization/pyext.C delete mode 100644 src/initialization/runpy.F create mode 100644 src/initialization/runpyC.C create mode 100644 src/initialization/runpyF.F90 diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 865aff53ce..97ca2b9e4b 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 @@ -653,9 +654,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & Y(I) = 0.0D0 EndDo - Call RunPy('test.py', 'my_test', X, Nx, Y, Ny) !This worked - !Call RunPy_3darray('test.py', 'my_test', h, is,ie,js,je,1,nz) - if (is_root_pe()) call MOM_mesg("MOM_initialize_state: Stop to debug!!"); stop + !Call RunPy('test.py', 'my_test', X, Nx, Y, Ny) !This worked + 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/pyext.C b/src/initialization/pyext.C deleted file mode 100644 index 2b08c26a38..0000000000 --- a/src/initialization/pyext.C +++ /dev/null @@ -1,138 +0,0 @@ -//Initially copied from https://github.com/wangsl/python-embedding -#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 py_test(const char *py_script_, const char *py_function_, - double *x, const int &nx, - double *y, const int &ny) -{ - if(!Py_IsInitialized()) Py_Initialize(); - - assert(Py_IsInitialized()); - assert(init_numpy()); - - 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. - //int x_in = 1; - //int y_in = 1; - //assert(PyObject_CallFunctionObjArgs(py_function, x_in, y_in, NULL)); //Segmentation fault - invalid memory reference - //assert(PyObject_CallFunctionObjArgs(py_function, x_in, NULL)); //Segmentation fault - invalid memory reference - //assert(PyObject_CallFunctionObjArgs(py_function, &x_in, NULL)); //Works but cannot access value. Ahaa, input needs to be a poniter - //assert(PyObject_CallFunctionObjArgs(py_function, dims, NULL)); //Works but cannot access value. invalid memory reference - /* The following worked fine. - long dims[2]; - dims[0] = 4; dims[1] = 4; - PyObject *pinValue, *pArgs,*pValue ; - pArgs = PyTuple_New(2); - for (int i = 0; i < 2; ++i) { - pinValue = PyLong_FromLong(dims[i]); - PyTuple_SetItem(pArgs, i, pinValue); - } - pValue = PyObject_CallFunctionObjArgs(py_function, pArgs, NULL); - if (pValue != NULL) { - printf("Result of call: %ld\n", PyLong_AsLong(pValue)); - Py_DECREF(pValue); - } - - double ddims[2]; - ddims[0] = 4.5; ddims[1] = 5.5; - pArgs = PyTuple_New(2); - for (int i = 0; i < 2; ++i) { - pinValue = PyFloat_FromDouble(ddims[i]); - PyTuple_SetItem(pArgs, i, pinValue); - } - pValue = PyObject_CallFunctionObjArgs(py_function, pArgs, NULL); - if (pValue != NULL) { - printf("Result of call: %f\n", PyFloat_AsDouble(pValue)); - Py_DECREF(pValue); - } - */ - const npy_intp dim_x [] = { nx }; - PyObject *x_py = PyArray_SimpleNewFromData(1, dim_x, NPY_DOUBLE, x); - const npy_intp dim_y [] = { ny }; - PyObject *y_py = PyArray_SimpleNewFromData(1, dim_y, NPY_DOUBLE, y); - 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; - Py_DECREF(y_py); y_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: PyTest -extern "C" void FORT(pytest)(const char *py_script, const int &len_py_script, - const char *py_function, const int &len_py_function, - double *x, const int &nx, - double *y, const int &ny) -{ - 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'; - - py_test(py_script_, py_function_, x, nx, y, ny); - - if(py_script_) { delete [] py_script_; py_script_ = 0; } - if(py_function_) { delete [] py_function_; py_function_ = 0; } -} - -// Fortran interface: PyFinalize -extern "C" void FORT(pyfinalize)() -{ - py_finalize(); -} - diff --git a/src/initialization/runpy.F b/src/initialization/runpy.F deleted file mode 100644 index 405763467b..0000000000 --- a/src/initialization/runpy.F +++ /dev/null @@ -1,50 +0,0 @@ -!Initially copied from https://github.com/wangsl/python-embedding -! - Subroutine RunPy(PyScript, PyFunction, X, Nx, Y, Ny) - Implicit None - Character*(*) PyScript, PyFunction - Real*8 X(*), Y(*) - Integer Nx, Ny - - Integer I, Len, LinEnd - External LinEnd - Integer LenPyScript, LenPyFunction - - LenPyScript = LinEnd(PyScript) - If(PyScript(LenPyScript-2:LenPyScript) .eq. '.py') - $ LenPyScript = LenPyScript-3 - - LenPyFunction = LinEnd(PyFunction) - - Call Flush(6) - - Call PyTest(PyScript, LenPyScript, - $ PyFunction, LenPyFunction, - $ X, Nx, Y, Ny) - - Return - End - - integer function linend(cline) - implicit real*8(a-h,o-z) -c -c function which returns the length of a character string, -c excluding trailing blanks, tabs and nulls. -c - character*(*) cline, blank - character*1 null, tab - parameter (blank=' ') -c - null = char(0) - tab = char(9) - do 5 i = len(cline), 1, -1 - if(cline(i:i).ne.blank.and. - $ cline(i:i).ne.null.and. - $ cline(i:i).ne.tab) goto 10 - 5 continue - linend = 0 - return - 10 linend = i - return - end - 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 index 8b7a263bc6..e69f0acf81 100644 --- a/src/initialization/test.py +++ b/src/initialization/test.py @@ -1,11 +1,41 @@ #!/bin/env python -#Initially copied from https://github.com/wangsl/python-embedding +#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_test0(x, y) : +def my_test_torch(x, y) : print(' From Python test: {}'.format(x.size)) for i in range(x.size) : x[i] += 1.0 @@ -27,18 +57,6 @@ def my_test0(x, y) : sys.stdout.flush() -def my_test(x) : - print(' **** From my_test1 ****',x) - sys.stdout.flush() - a=sum(x) - return a - -def my_test2(x,y) : - print(' **** From my_test2 ****',x,y) - sys.stdout.flush() - a=sum(x) - b=sum(y) - return a,b if __name__ == '__main__' : import numpy as np From c96a37c16b2e17620b49c4c217cae463aa849409 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Thu, 7 Oct 2021 18:04:27 -0400 Subject: [PATCH 4/4] Operate with a Python tool on a given MOM6 array - This update shows an example to use Python tool in model runtime to operate on a MOM6 array and calculate/plot something. - This works for a demo in single column model that can be obtained from https://github.com/nikizadehgfdl/MOM6_OBGC_examples --- src/initialization/MOM_state_initialization.F90 | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 97ca2b9e4b..26077a2d77 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -185,10 +185,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & # include "version_variable.h" integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - Integer Nx, Ny - Parameter(Nx = 15) - Parameter(Ny = 15) - Real*8 X(Nx), Y(Ny) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB @@ -649,12 +645,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & endif if (is_root_pe()) call MOM_mesg("MOM_initialize_state: Testing Python embedding ...") - Do I = 1, Nx - X(I) = DBLE(I) - Y(I) = 0.0D0 - EndDo - - !Call RunPy('test.py', 'my_test', X, Nx, Y, Ny) !This worked 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