From 4d0d1a9677d78a8e1b75416be3299fa291565432 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Mon, 2 Aug 2021 17:03:13 -0700 Subject: [PATCH 1/8] First stab at Fortran-friendly named BC functions. --- demo/transient/transient_snes_mpfaof90.F90 | 12 +- include/tdycore.h | 6 + src/f90-mod/tdycoremod.F90 | 125 +++++++++++++++++++++ src/tdyconditions.c | 86 ++++++++++++++ src/tdycore.c | 16 +++ 5 files changed, 242 insertions(+), 3 deletions(-) diff --git a/demo/transient/transient_snes_mpfaof90.F90 b/demo/transient/transient_snes_mpfaof90.F90 index 127dcea0..ebdd6c2d 100644 --- a/demo/transient/transient_snes_mpfaof90.F90 +++ b/demo/transient/transient_snes_mpfaof90.F90 @@ -84,13 +84,13 @@ subroutine MaterialPropAlpha_PFLOTRAN(alpha) PetscReal, intent(out) :: alpha alpha = 1.d-4 end subroutine MaterialPropAlpha_PFLOTRAN - + subroutine MaterialPropM_PFLOTRAN(m) implicit none PetscReal, intent(out) :: m m = 0.3d0 end subroutine MaterialPropM_PFLOTRAN - + subroutine ResidualSat_PFLOTRAN(resSat) implicit none PetscReal resSat @@ -160,6 +160,11 @@ program main call TDyInit(ierr); CHKERRA(ierr); + + ! Register some functions. + call TDyRegisterPressureFn("p0", PressureFunction, ierr) + CHKERRA(ierr); + call TDyCreate(tdy, ierr); CHKERRA(ierr); call TDySetDiscretizationMethod(tdy,MPFA_O,ierr); @@ -340,7 +345,8 @@ program main end if if (bc_type == MPFAO_DIRICHLET_BC .OR. bc_type == MPFAO_SEEPAGE_BC ) then - call TDySetBoundaryPressureFn(tdy,PressureFunction,0,ierr); +! call TDySetBoundaryPressureFn(tdy,PressureFunction,0,ierr); + call TDySelectBoundaryPressureFn(tdy,"p0",ierr); CHKERRQ(ierr) endif diff --git a/include/tdycore.h b/include/tdycore.h index fcfbca07..491063b4 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -102,8 +102,14 @@ PETSC_EXTERN PetscErrorCode TDySetSoilDensityFunction(TDy,PetscErrorCode(*)(TDy, PETSC_EXTERN PetscErrorCode TDySetSoilSpecificHeatFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); // Set boundary and source-sink: via PETSc operations +PetscErrorCode TDyRegisterPressureFn(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*)); +PetscErrorCode TDyRegisterTemperatureFn(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*)); +PetscErrorCode TDyRegisterVelocityFn(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*)); PETSC_EXTERN PetscErrorCode TDySetForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); PETSC_EXTERN PetscErrorCode TDySetEnergyForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); +PETSC_EXTERN PetscErrorCode TDySelectBoundaryPressureFn(TDy,const char*,void*); +PETSC_EXTERN PetscErrorCode TDySelectBoundaryTemperatureFn(TDy,const char*,void*); +PETSC_EXTERN PetscErrorCode TDySelectBoundaryVelocityFn(TDy,const char*,void*); PETSC_EXTERN PetscErrorCode TDySetBoundaryPressureFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); PETSC_EXTERN PetscErrorCode TDySetBoundaryTemperatureFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); PETSC_EXTERN PetscErrorCode TDySetBoundaryVelocityFn(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index a5b060c1..bd600219 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -450,4 +450,129 @@ subroutine TDyInit(ierr) call TDyInitNoArguments(ierr) end subroutine TDyInit + subroutine TDyRegisterPressureFn(func_name, func, ierr) + use, intrinsic :: iso_c_binding + implicit none + character(len=*), intent(in) :: func_name + interface + subroutine func(tdy, x, pressure, dummy, ierr) + use tdycoredef + use petscvec + TDy :: tdy + PetscReal, intent(in) :: x(3) + PetscReal, intent(out) :: pressure + integer :: dummy(*) + PetscErrorCode :: ierr + end subroutine + + function RegisterPressureFn(name, func) bind (c, name="TDyRegisterPressureFn") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: name + type(c_funptr), value :: func + integer(c_int) :: ierr + end function + end interface + integer, intent(out) :: ierr + + ierr = RegisterPressureFn(f_to_c_string(func_name), c_funloc(func)) + end subroutine + + subroutine TDyRegisterTemperatureFn(func_name, func, ierr) + use, intrinsic :: iso_c_binding + implicit none + character(len=*), intent(in) :: func_name + interface + subroutine func(tdy, x, temperature, dummy, ierr) + use tdycoredef + use petscvec + TDy :: tdy + PetscReal, intent(in) :: x(3) + PetscReal, intent(out) :: temperature + integer :: dummy(*) + PetscErrorCode :: ierr + end subroutine + + function RegisterTemperatureFn(name, func) bind (c, name="TDyRegisterTemperatureFn") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: name + type(c_funptr), value :: func + integer(c_int) :: ierr + end function + end interface + integer, intent(out) :: ierr + + ierr = RegisterTemperatureFn(f_to_c_string(func_name), c_funloc(func)) + end subroutine + + subroutine TDyRegisterVelocityFn(func_name, func, ierr) + use, intrinsic :: iso_c_binding + implicit none + character(len=*), intent(in) :: func_name + interface + subroutine func(tdy, x, velocity, dummy, ierr) + use tdycoredef + use petscvec + TDy :: tdy + PetscReal, intent(in) :: x(3) + PetscReal, intent(out) :: velocity + integer :: dummy(*) + PetscErrorCode :: ierr + end subroutine + + function RegisterVelocityFn(name, func) bind (c, name="TDyRegisterVelocityFn") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: name + type(c_funptr), value :: func + integer(c_int) :: ierr + end function + end interface + integer, intent(out) :: ierr + + ierr = RegisterVelocityFn(f_to_c_string(func_name), c_funloc(func)) + end subroutine + + subroutine TDySelectBoundaryPressureFn(tdy, func_name, ierr) + use, intrinsic :: iso_c_binding + use tdycoredef + implicit none + TDy, target :: tdy + character(len=*), intent(in) :: func_name + integer, intent(out) :: ierr + interface + function SelectBoundaryPressureFn(tdy, name, ctx) bind (c, name="TDySelectBoundaryPressureFn") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: tdy + type(c_ptr), value :: name + type(c_ptr), value :: ctx + integer(c_int) :: ierr + end function + end interface + + ierr = SelectBoundaryPressureFn(c_loc(tdy), f_to_c_string(func_name), c_null_ptr) + end subroutine + + function f_to_c_string(f_string) result(c_string) + use, intrinsic :: iso_c_binding + implicit none + character(len=*), target :: f_string + character(len=:), pointer :: f_ptr + type(c_ptr) :: c_string + + interface + function new_c_string(f_str_ptr, f_str_len) bind (c) result(c_string) + use, intrinsic :: iso_c_binding + type(c_ptr), value :: f_str_ptr + integer(c_int), value :: f_str_len + type(c_ptr) :: c_string + end function new_c_string + end interface + + f_ptr => f_string + c_string = new_c_string(c_loc(f_ptr), len(f_string)) + end function f_to_c_string + end module tdycore diff --git a/src/tdyconditions.c b/src/tdyconditions.c index b10cf179..8193c2f6 100644 --- a/src/tdyconditions.c +++ b/src/tdyconditions.c @@ -1,9 +1,53 @@ #include +#include /* Boundary and source-sink conditions are set by PETSc operations */ +// Here's a registry of pressure, temperature, and velocity functions. +typedef PetscErrorCode(*BCFunction)(TDy, PetscReal*, PetscReal*, void*); +KHASH_MAP_INIT_STR(TDY_BCFUNC_MAP, BCFunction) +static khash_t(TDY_BCFUNC_MAP)* pressure_funcs_ = NULL; +static khash_t(TDY_BCFUNC_MAP)* temperature_funcs_ = NULL; +static khash_t(TDY_BCFUNC_MAP)* velocity_funcs_ = NULL; + +PetscErrorCode TDyRegisterPressureFn(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { + PetscFunctionBegin; + if (pressure_funcs_ == NULL) { + pressure_funcs_ = kh_init(TDY_BCFUNC_MAP); + } + + int retval; + khiter_t iter = kh_put(TDY_BCFUNC_MAP, pressure_funcs_, name, &retval); + kh_val(pressure_funcs_, iter) = f; + PetscFunctionReturn(0); +} + +PetscErrorCode TDyRegisterTemperatureFn(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { + PetscFunctionBegin; + if (temperature_funcs_ == NULL) { + temperature_funcs_ = kh_init(TDY_BCFUNC_MAP); + } + + int retval; + khiter_t iter = kh_put(TDY_BCFUNC_MAP, temperature_funcs_, name, &retval); + kh_val(temperature_funcs_, iter) = f; + PetscFunctionReturn(0); +} + +PetscErrorCode TDyRegisterVelocityFn(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { + PetscFunctionBegin; + if (velocity_funcs_ == NULL) { + velocity_funcs_ = kh_init(TDY_BCFUNC_MAP); + } + + int retval; + khiter_t iter = kh_put(TDY_BCFUNC_MAP, velocity_funcs_, name, &retval); + kh_val(velocity_funcs_, iter) = f; + PetscFunctionReturn(0); +} + PetscErrorCode TDySetForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*),void *ctx) { PetscFunctionBegin; if (f) tdy->ops->computeforcing = f; @@ -18,6 +62,48 @@ PetscErrorCode TDySetEnergyForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,Petsc PetscFunctionReturn(0); } +PetscErrorCode TDySelectBoundaryPressureFn(TDy tdy, const char* name, void* ctx) { + PetscFunctionBegin; + int ierr; + + khiter_t iter = kh_get(TDY_BCFUNC_MAP, pressure_funcs_, name); + if (iter != kh_end(pressure_funcs_)) { // found it! + BCFunction f = kh_val(pressure_funcs_, iter); + ierr = TDySetBoundaryPressureFn(tdy, f, ctx); CHKERRQ(ierr); + } else { + printf("Uh oh!\n"); // TODO: handle error condition + } + PetscFunctionReturn(0); +} + +PetscErrorCode TDySelectBoundaryTemperatureFn(TDy tdy, const char* name, void* ctx) { + PetscFunctionBegin; + int ierr; + + khiter_t iter = kh_get(TDY_BCFUNC_MAP, temperature_funcs_, name); + if (iter != kh_end(temperature_funcs_)) { // found it! + BCFunction f = kh_val(temperature_funcs_, iter); + ierr = TDySetBoundaryTemperatureFn(tdy, f, ctx); CHKERRQ(ierr); + } else { + printf("Uh oh!\n"); // TODO: handle error condition + } + PetscFunctionReturn(0); +} + +PetscErrorCode TDySelectBoundaryVelocityFn(TDy tdy, const char* name, void* ctx) { + PetscFunctionBegin; + int ierr; + + khiter_t iter = kh_get(TDY_BCFUNC_MAP, velocity_funcs_, name); + if (iter != kh_end(velocity_funcs_)) { // found it! + BCFunction f = kh_val(velocity_funcs_, iter); + ierr = TDySetBoundaryVelocityFn(tdy, f, ctx); CHKERRQ(ierr); + } else { + printf("Uh oh!\n"); // TODO: handle error condition + } + PetscFunctionReturn(0); +} + PetscErrorCode TDySetBoundaryTemperatureFn(TDy tdy, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*),void *ctx) { PetscFunctionBegin; if (f) tdy->ops->compute_boundary_temperature = f; diff --git a/src/tdycore.c b/src/tdycore.c index ba66fd12..746f468b 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -1460,3 +1460,19 @@ PetscErrorCode TDyCreateJacobian(TDy tdy) { PetscFunctionReturn(0); } +// Here's a registry of C-backed strings created from Fortran. +// Currently we limit one to 1024 such strings. +static char* fortran_strings_[1024]; +static int num_fortran_strings_ = 0; + +// Returns a newly-allocated C string for the given Fortran string pointer with +// the given length. Resources for this string are managed by the running model. +const char* new_c_string(char* f_str_ptr, int f_str_len) { + char* new_string = malloc(sizeof(char) * (f_str_len+1)); + memcpy(new_string, f_str_ptr, sizeof(char) * f_str_len); + new_string[f_str_len] = '\0'; + fortran_strings_[num_fortran_strings_] = new_string; + ++num_fortran_strings_; + return (const char*)new_string; +} + From 4453215e40f3a52c87c504ec845c554d3e4e5609 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 3 Aug 2021 09:51:51 -0700 Subject: [PATCH 2/8] Added a registry for shutdown functions and for Fortran strings. --- include/tdycore.h | 3 +- src/f90-mod/tdycoremod.F90 | 18 +++++----- src/tdycore.c | 73 ++++++++++++++++++++++++++++++-------- src/tdytimers.c | 11 ++++++ 4 files changed, 81 insertions(+), 24 deletions(-) diff --git a/include/tdycore.h b/include/tdycore.h index 491063b4..6f1e851c 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -74,8 +74,9 @@ PETSC_EXTERN PetscClassId TDY_CLASSID; PETSC_EXTERN PetscLogEvent TDy_ComputeSystem; /* ---------------------------------------------------------------- */ -PETSC_EXTERN PetscErrorCode TDyInit(int argc, char* argv[]); +PETSC_EXTERN PetscErrorCode TDyInit(int, char*[]); PETSC_EXTERN PetscErrorCode TDyInitNoArguments(void); +PETSC_EXTERN PetscErrorCode TDyOnFinalize(void (*)(void)); PETSC_EXTERN PetscErrorCode TDyFinalize(void); PETSC_EXTERN PetscErrorCode TDyCreate(TDy*); diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index bd600219..6ae83bf9 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -475,7 +475,7 @@ function RegisterPressureFn(name, func) bind (c, name="TDyRegisterPressureFn") r end interface integer, intent(out) :: ierr - ierr = RegisterPressureFn(f_to_c_string(func_name), c_funloc(func)) + ierr = RegisterPressureFn(FtoCString(func_name), c_funloc(func)) end subroutine subroutine TDyRegisterTemperatureFn(func_name, func, ierr) @@ -503,7 +503,7 @@ function RegisterTemperatureFn(name, func) bind (c, name="TDyRegisterTemperature end interface integer, intent(out) :: ierr - ierr = RegisterTemperatureFn(f_to_c_string(func_name), c_funloc(func)) + ierr = RegisterTemperatureFn(FtoCString(func_name), c_funloc(func)) end subroutine subroutine TDyRegisterVelocityFn(func_name, func, ierr) @@ -531,7 +531,7 @@ function RegisterVelocityFn(name, func) bind (c, name="TDyRegisterVelocityFn") r end interface integer, intent(out) :: ierr - ierr = RegisterVelocityFn(f_to_c_string(func_name), c_funloc(func)) + ierr = RegisterVelocityFn(FtoCString(func_name), c_funloc(func)) end subroutine subroutine TDySelectBoundaryPressureFn(tdy, func_name, ierr) @@ -552,10 +552,10 @@ function SelectBoundaryPressureFn(tdy, name, ctx) bind (c, name="TDySelectBounda end function end interface - ierr = SelectBoundaryPressureFn(c_loc(tdy), f_to_c_string(func_name), c_null_ptr) + ierr = SelectBoundaryPressureFn(c_loc(tdy), FtoCString(func_name), c_null_ptr) end subroutine - function f_to_c_string(f_string) result(c_string) + function FtoCString(f_string) result(c_string) use, intrinsic :: iso_c_binding implicit none character(len=*), target :: f_string @@ -563,16 +563,16 @@ function f_to_c_string(f_string) result(c_string) type(c_ptr) :: c_string interface - function new_c_string(f_str_ptr, f_str_len) bind (c) result(c_string) + function NewCString(f_str_ptr, f_str_len) bind (c, name="NewCString") result(c_string) use, intrinsic :: iso_c_binding type(c_ptr), value :: f_str_ptr integer(c_int), value :: f_str_len type(c_ptr) :: c_string - end function new_c_string + end function NewCString end interface f_ptr => f_string - c_string = new_c_string(c_loc(f_ptr), len(f_string)) - end function f_to_c_string + c_string = NewCString(c_loc(f_ptr), len(f_string)) + end function FtoCString end module tdycore diff --git a/src/tdycore.c b/src/tdycore.c index 746f468b..3620109e 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -140,13 +140,38 @@ PetscBool TDyInitialized(void) { return TDyPackageInitialized; } +// A registry of functions called at shutdown. This can be used by all +// subsystems via TDyOnFinalize(). +typedef void (*ShutdownFunc)(void); +static ShutdownFunc *shutdown_funcs_ = NULL; +static int num_shutdown_funcs_ = 0; +static int shutdown_funcs_cap_ = 0; + +/// Call this to register a shutdown function that is called during TDyFinalize. +PetscErrorCode TDyOnFinalize(void (*shutdown_func)(void)) { + if (shutdown_funcs_ == NULL) { + shutdown_funcs_cap_ = 32; + shutdown_funcs_ = malloc(sizeof(ShutdownFunc) * shutdown_funcs_cap_); + } else if (num_shutdown_funcs_ == shutdown_funcs_cap_) { // need more space! + shutdown_funcs_cap_ *= 2; + shutdown_funcs_ = realloc(shutdown_funcs_, sizeof(ShutdownFunc) * shutdown_funcs_cap_); + } + shutdown_funcs_[num_shutdown_funcs_] = shutdown_func; + ++num_shutdown_funcs_; + + PetscFunctionReturn(0); +} + PetscErrorCode TDyFinalize() { PetscFunctionBegin; - // Dump timing information before we leave. - TDyWriteTimingProfile("tdycore_profile.csv"); - - TDyDestroyTimers(); + // Call shutdown functions and destroy the list. + if (shutdown_funcs_ != NULL) { + for (int i = 0; i < num_shutdown_funcs_; ++i) { + shutdown_funcs_[i](); + } + free(shutdown_funcs_); + } // Finalize PETSc. PetscFinalize(); @@ -1461,18 +1486,38 @@ PetscErrorCode TDyCreateJacobian(TDy tdy) { } // Here's a registry of C-backed strings created from Fortran. -// Currently we limit one to 1024 such strings. -static char* fortran_strings_[1024]; -static int num_fortran_strings_ = 0; +KHASH_SET_INIT_STR(TDY_STRING_SET) +static khash_t(TDY_STRING_SET)* fortran_strings_ = NULL; + +// This function is called on finalization to destroy the Fortran Ń•tring +// registry. +static void DestroyFortranStrings() { + kh_destroy(TDY_STRING_SET, fortran_strings_); +} // Returns a newly-allocated C string for the given Fortran string pointer with // the given length. Resources for this string are managed by the running model. -const char* new_c_string(char* f_str_ptr, int f_str_len) { - char* new_string = malloc(sizeof(char) * (f_str_len+1)); - memcpy(new_string, f_str_ptr, sizeof(char) * f_str_len); - new_string[f_str_len] = '\0'; - fortran_strings_[num_fortran_strings_] = new_string; - ++num_fortran_strings_; - return (const char*)new_string; +const char* NewCString(char* f_str_ptr, int f_str_len) { + if (fortran_strings_ == NULL) { + fortran_strings_ = kh_init(TDY_STRING_SET); + TDyOnFinalize(DestroyFortranStrings); + } + + // Copy the Fortran character array to a C string. + char c_str[f_str_len+1]; + memcpy(c_str, f_str_ptr, sizeof(char) * f_str_len); + c_str[f_str_len] = '\0'; + + // Does this string already exist? + khiter_t iter = kh_get(TDY_STRING_SET, fortran_strings_, c_str); + if (iter != kh_end(fortran_strings_)) { // yep + return kh_key(fortran_strings_, iter); + } else { // nope + int retval; + const char* str = malloc(sizeof(char) * (f_str_len+1)); + strcpy((char*)str, c_str); + iter = kh_put(TDY_STRING_SET, fortran_strings_, str, &retval); + return str; + } } diff --git a/src/tdytimers.c b/src/tdytimers.c index 3a2649b5..cc62451a 100644 --- a/src/tdytimers.c +++ b/src/tdytimers.c @@ -24,6 +24,14 @@ khash_t(TDY_PROFILING_MD_MAP)* TDY_PROFILING_METADATA = NULL; // Are timers enabled? static PetscBool timersEnabled_ = PETSC_FALSE; +// This function is called at finalization time. +static void ShutDownTimers() { + // Dump timing information before we leave. + TDyWriteTimingProfile("tdycore_profile.csv"); + + TDyDestroyTimers(); +} + PetscErrorCode TDyInitTimers() { // Register timers table. if (TDY_TIMERS == NULL) { @@ -46,6 +54,9 @@ PetscErrorCode TDyInitTimers() { ierr = TDyAddProfilingStage("TDycore Stepping"); CHKERRQ(ierr); ierr = TDyAddProfilingStage("TDycore I/O"); CHKERRQ(ierr); + // Register our shutdown function. + TDyOnFinalize(ShutDownTimers); + return ierr; } From c8fb2862a731523c2f1c479d099edb41a5492da5 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Tue, 3 Aug 2021 16:44:22 -0700 Subject: [PATCH 3/8] Switched to PETSc one-off Fortran bindings. --- src/f90-mod/tdycoremod.F90 | 125 ----------------------------------- src/interface/ftn/tdycoref.c | 12 ++++ 2 files changed, 12 insertions(+), 125 deletions(-) diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index 6ae83bf9..a5b060c1 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -450,129 +450,4 @@ subroutine TDyInit(ierr) call TDyInitNoArguments(ierr) end subroutine TDyInit - subroutine TDyRegisterPressureFn(func_name, func, ierr) - use, intrinsic :: iso_c_binding - implicit none - character(len=*), intent(in) :: func_name - interface - subroutine func(tdy, x, pressure, dummy, ierr) - use tdycoredef - use petscvec - TDy :: tdy - PetscReal, intent(in) :: x(3) - PetscReal, intent(out) :: pressure - integer :: dummy(*) - PetscErrorCode :: ierr - end subroutine - - function RegisterPressureFn(name, func) bind (c, name="TDyRegisterPressureFn") result(ierr) - use, intrinsic :: iso_c_binding - implicit none - type(c_ptr), value :: name - type(c_funptr), value :: func - integer(c_int) :: ierr - end function - end interface - integer, intent(out) :: ierr - - ierr = RegisterPressureFn(FtoCString(func_name), c_funloc(func)) - end subroutine - - subroutine TDyRegisterTemperatureFn(func_name, func, ierr) - use, intrinsic :: iso_c_binding - implicit none - character(len=*), intent(in) :: func_name - interface - subroutine func(tdy, x, temperature, dummy, ierr) - use tdycoredef - use petscvec - TDy :: tdy - PetscReal, intent(in) :: x(3) - PetscReal, intent(out) :: temperature - integer :: dummy(*) - PetscErrorCode :: ierr - end subroutine - - function RegisterTemperatureFn(name, func) bind (c, name="TDyRegisterTemperatureFn") result(ierr) - use, intrinsic :: iso_c_binding - implicit none - type(c_ptr), value :: name - type(c_funptr), value :: func - integer(c_int) :: ierr - end function - end interface - integer, intent(out) :: ierr - - ierr = RegisterTemperatureFn(FtoCString(func_name), c_funloc(func)) - end subroutine - - subroutine TDyRegisterVelocityFn(func_name, func, ierr) - use, intrinsic :: iso_c_binding - implicit none - character(len=*), intent(in) :: func_name - interface - subroutine func(tdy, x, velocity, dummy, ierr) - use tdycoredef - use petscvec - TDy :: tdy - PetscReal, intent(in) :: x(3) - PetscReal, intent(out) :: velocity - integer :: dummy(*) - PetscErrorCode :: ierr - end subroutine - - function RegisterVelocityFn(name, func) bind (c, name="TDyRegisterVelocityFn") result(ierr) - use, intrinsic :: iso_c_binding - implicit none - type(c_ptr), value :: name - type(c_funptr), value :: func - integer(c_int) :: ierr - end function - end interface - integer, intent(out) :: ierr - - ierr = RegisterVelocityFn(FtoCString(func_name), c_funloc(func)) - end subroutine - - subroutine TDySelectBoundaryPressureFn(tdy, func_name, ierr) - use, intrinsic :: iso_c_binding - use tdycoredef - implicit none - TDy, target :: tdy - character(len=*), intent(in) :: func_name - integer, intent(out) :: ierr - interface - function SelectBoundaryPressureFn(tdy, name, ctx) bind (c, name="TDySelectBoundaryPressureFn") result(ierr) - use, intrinsic :: iso_c_binding - implicit none - type(c_ptr), value :: tdy - type(c_ptr), value :: name - type(c_ptr), value :: ctx - integer(c_int) :: ierr - end function - end interface - - ierr = SelectBoundaryPressureFn(c_loc(tdy), FtoCString(func_name), c_null_ptr) - end subroutine - - function FtoCString(f_string) result(c_string) - use, intrinsic :: iso_c_binding - implicit none - character(len=*), target :: f_string - character(len=:), pointer :: f_ptr - type(c_ptr) :: c_string - - interface - function NewCString(f_str_ptr, f_str_len) bind (c, name="NewCString") result(c_string) - use, intrinsic :: iso_c_binding - type(c_ptr), value :: f_str_ptr - integer(c_int), value :: f_str_len - type(c_ptr) :: c_string - end function NewCString - end interface - - f_ptr => f_string - c_string = NewCString(c_loc(f_ptr), len(f_string)) - end function FtoCString - end module tdycore diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index 97c591cd..bcf78c52 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -39,6 +39,12 @@ #define tdysetpermeabilityfunction_ TDYSETPERMEABILITYFUNCTION #define tdysetresidualsaturationfunction_ TDYSETRESIDUALSATURATIONFUNCTION #define tdysetforcingfunction_ TDYSETFORCINGFUNCTION +#define tdyregisterpressurefn_ TDYREGISTERPRESSUREFN +#define tdyregistertemperaturefn_ TDYREGISTERTEMPERATUREFN +#define tdyregistervelocityfn_ TDYREGISTERVELOCITYFN +#define tdyselectboundarypressurefn_ TDYSELECTBOUNDARYPRESSUREFN +#define tdyselectboundarytemperaturefn_ TDYSELECTBOUNDARYTEMPERATUREFN +#define tdyselectboundaryvelocityfn_ TDYSELECTBOUNDARYVELOCITYFN #define tdysetboundarypressurefn_ TDYSETBOUNDARYPRESSUREFN #define tdysetboundaryvelocityfn_ TDYSETBOUNDARYVELOCITYFN #define tdysetporosityvalueslocal0_ TDYSETPOROSITYVALUESLOCAL0 @@ -101,6 +107,12 @@ #define tdysetpermeabilityfunction_ tdysetpermeabilityfunction #define tdysetresidualsaturationfunction_ tdysetresidualsaturationfunction #define tdysetforcingfunction_ tdysetforcingfunction +#define tdyregisterpressurefn_ tdyregisterpressurefn +#define tdyregistertemperaturefn_ tdyregistertemperaturefn +#define tdyregistervelocityfn_ tdyregistervelocityfn +#define tdyselectboundarypressurefn_ tdyselectboundarypressurefn +#define tdyselectboundarytemperaturefn_ tdyselectboundarytemperaturefn +#define tdyselectboundaryvelocityfn_ tdyselectboundaryvelocityfn #define tdysetboundarypressurefn_ tdysetboundarypressurefn #define tdysetboundaryvelocityfn_ tdysetboundaryvelocityfn #define tdysetporosityvalueslocal0_ tdysetporosityvalueslocal0 From 1b8377e9e7104f81506fc45c3cf78048e0bea351 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Wed, 4 Aug 2021 10:39:23 -0700 Subject: [PATCH 4/8] "Saving" my "game." --- src/interface/ftn/tdycoref.c | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index bcf78c52..a8e67da0 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -142,6 +142,7 @@ #define tdydestroy_ tdydestroy #endif +// Legacy Fortran function pointer registry static struct { PetscFortranCallbackId porosity; PetscFortranCallbackId permeability; @@ -154,6 +155,29 @@ static struct { #endif } _cb; +// This struct and wrapper allow us to call a Fortran function (with an +// additional ierr argument) from C. This uses our dynamic C function registry +// instead of _cb above. +typedef struct { + // Fortran function pointer + PetscErrorCode (*fort_func)(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode*),void *ctx,PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptr)) + // Fortran function context + void* fort_ctx; + // Double-pointer (PGI?) compiler hack +#if defined(PETSC_HAVE_F90_2PTR_ARG) + PetscFortranCallbackId fort_pgiptr; +#endif +} TDyFortFunc; +static TDyFortFunc* createfortfunc(PetscErrorCode (*func)(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode*),void *ctx,PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptr)) +static PetscErrorCode tdyfunc_wrapper(TDy tdy,PetscReal *x,PetscReal *f,void *ctx) +{ + TDyFortFunc* ffunc = (TDyFortFunc*)ctx; +#if defined(PETSC_HAVE_F90_2PTR_ARG) + void* ptr; + PetscObjectGetFortranCallback((PetscObject)tdy,PETSC_FORTRAN_CALLBACK_CLASS,ffunc->fort_pgiptr,NULL,&ptr); +#endif + PetscObjectUseFortranCallback(tdy,ffunc->fort_func,(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode* PETSC_F90_2PTR_PROTO_NOVAR),(&tdy,x,f,ffunc->fort_ctx,&ierr PETSC_F90_2PTR_PARAM(ptr))); +} #if defined(__cplusplus) extern "C" { @@ -175,6 +199,16 @@ PETSC_EXTERN void tdyfinalize_(int *__ierr){ } #endif +#if defined(__cplusplus) +extern "C" { +#endif +PETSC_EXTERN void tdyregisterpressurefn(char* name, PetscErrorCode (*func)(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode*),void *ctx,PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptr)) +*__ierr = TDyRegisterPressureFn((TDy)PetscToPointer((tdy))); +} +#if defined(__cplusplus) +} +#endif + #if defined(__cplusplus) extern "C" { #endif From adcb822c6a75a0c2c123c089b924d59bb8331bd3 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Wed, 4 Aug 2021 17:02:37 -0700 Subject: [PATCH 5/8] Trying to work around C/Fortran interoperability issues. --- demo/transient/transient_snes_mpfaof90.F90 | 2 +- include/tdycore.h | 4 +- src/f90-mod/tdycoremod.F90 | 80 ++++++++++++++++++++++ src/interface/ftn/tdycoref.c | 40 ----------- src/tdyconditions.c | 67 +++++++----------- 5 files changed, 106 insertions(+), 87 deletions(-) diff --git a/demo/transient/transient_snes_mpfaof90.F90 b/demo/transient/transient_snes_mpfaof90.F90 index ebdd6c2d..a61d7b60 100644 --- a/demo/transient/transient_snes_mpfaof90.F90 +++ b/demo/transient/transient_snes_mpfaof90.F90 @@ -162,7 +162,7 @@ program main CHKERRA(ierr); ! Register some functions. - call TDyRegisterPressureFn("p0", PressureFunction, ierr) + call TDyRegisterFunction("p0", PressureFunction, ierr) CHKERRA(ierr); call TDyCreate(tdy, ierr); diff --git a/include/tdycore.h b/include/tdycore.h index 6f1e851c..e5d0a8f2 100644 --- a/include/tdycore.h +++ b/include/tdycore.h @@ -103,9 +103,7 @@ PETSC_EXTERN PetscErrorCode TDySetSoilDensityFunction(TDy,PetscErrorCode(*)(TDy, PETSC_EXTERN PetscErrorCode TDySetSoilSpecificHeatFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); // Set boundary and source-sink: via PETSc operations -PetscErrorCode TDyRegisterPressureFn(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*)); -PetscErrorCode TDyRegisterTemperatureFn(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*)); -PetscErrorCode TDyRegisterVelocityFn(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*)); +PetscErrorCode TDyRegisterFunction(const char*, PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*)); PETSC_EXTERN PetscErrorCode TDySetForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); PETSC_EXTERN PetscErrorCode TDySetEnergyForcingFunction(TDy,PetscErrorCode(*)(TDy,PetscReal*,PetscReal*,void*),void*); PETSC_EXTERN PetscErrorCode TDySelectBoundaryPressureFn(TDy,const char*,void*); diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index a5b060c1..af49edce 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -437,6 +437,17 @@ subroutine TDyPostSolveSNESSolver(a,b,z) end subroutine TDyPostSolveSNESSolver end interface + abstract interface + subroutine TDyFunction(tdy, x, f, dummy, ierr) + use tdycoredef + TDy :: tdy + PetscReal, intent(in) :: x(3) + PetscReal, intent(out) :: f + integer :: dummy(*) + PetscErrorCode :: ierr + end subroutine + end interface + contains subroutine TDyInit(ierr) @@ -450,4 +461,73 @@ subroutine TDyInit(ierr) call TDyInitNoArguments(ierr) end subroutine TDyInit + subroutine TDyRegisterFunction(name, func, ierr) + use, intrinsic :: iso_c_binding + implicit none + character(len=*), intent(in) :: name + procedure(TDyFunction) :: func + PetscErrorCode :: ierr + + interface + function RegisterFn(name, func) bind (c, name="TDyRegisterFunction") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: name + type(c_funptr), value :: func + integer(c_int) :: ierr + end function + end interface + + ierr = RegisterFn(FtoCString(name), c_funloc(func)) + end subroutine + + subroutine TDySelectBoundaryPressureFn(tdy, name, ierr) + use, intrinsic :: iso_c_binding + use tdycoredef + implicit none + TDy :: tdy + character(len=*), intent(in) :: name + PetscErrorCode :: ierr + + type(c_funptr) :: c_func + procedure(TDyFunction), pointer :: f_func + + interface + function GetFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: name + type(c_funptr) :: c_func + integer(c_int) :: ierr + end function + end interface + + ierr = GetFn(FtoCString(name), c_func) + call c_f_procpointer(c_func, f_func) + call TDySetBoundaryPressureFn(tdy, f_func, c_null_ptr) + end subroutine + + ! Here's a function that converts a Fortran string to a C string and + ! stashes it in TDycore's Fortran string registry. This allows us to + ! create more expressive (and standard) Fortran interfaces. + function FtoCString(f_string) result(c_string) + use, intrinsic :: iso_c_binding + implicit none + character(len=*), target :: f_string + character(len=:), pointer :: f_ptr + type(c_ptr) :: c_string + + interface + function NewCString(f_str_ptr, f_str_len) bind (c, name="NewCString") result(c_string) + use, intrinsic :: iso_c_binding + type(c_ptr), value :: f_str_ptr + integer(c_int), value :: f_str_len + type(c_ptr) :: c_string + end function NewCString + end interface + + f_ptr => f_string + c_string = NewCString(c_loc(f_ptr), len(f_string)) + end function FtoCString + end module tdycore diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index a8e67da0..38efc9cd 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -39,9 +39,6 @@ #define tdysetpermeabilityfunction_ TDYSETPERMEABILITYFUNCTION #define tdysetresidualsaturationfunction_ TDYSETRESIDUALSATURATIONFUNCTION #define tdysetforcingfunction_ TDYSETFORCINGFUNCTION -#define tdyregisterpressurefn_ TDYREGISTERPRESSUREFN -#define tdyregistertemperaturefn_ TDYREGISTERTEMPERATUREFN -#define tdyregistervelocityfn_ TDYREGISTERVELOCITYFN #define tdyselectboundarypressurefn_ TDYSELECTBOUNDARYPRESSUREFN #define tdyselectboundarytemperaturefn_ TDYSELECTBOUNDARYTEMPERATUREFN #define tdyselectboundaryvelocityfn_ TDYSELECTBOUNDARYVELOCITYFN @@ -107,9 +104,6 @@ #define tdysetpermeabilityfunction_ tdysetpermeabilityfunction #define tdysetresidualsaturationfunction_ tdysetresidualsaturationfunction #define tdysetforcingfunction_ tdysetforcingfunction -#define tdyregisterpressurefn_ tdyregisterpressurefn -#define tdyregistertemperaturefn_ tdyregistertemperaturefn -#define tdyregistervelocityfn_ tdyregistervelocityfn #define tdyselectboundarypressurefn_ tdyselectboundarypressurefn #define tdyselectboundarytemperaturefn_ tdyselectboundarytemperaturefn #define tdyselectboundaryvelocityfn_ tdyselectboundaryvelocityfn @@ -155,30 +149,6 @@ static struct { #endif } _cb; -// This struct and wrapper allow us to call a Fortran function (with an -// additional ierr argument) from C. This uses our dynamic C function registry -// instead of _cb above. -typedef struct { - // Fortran function pointer - PetscErrorCode (*fort_func)(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode*),void *ctx,PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptr)) - // Fortran function context - void* fort_ctx; - // Double-pointer (PGI?) compiler hack -#if defined(PETSC_HAVE_F90_2PTR_ARG) - PetscFortranCallbackId fort_pgiptr; -#endif -} TDyFortFunc; -static TDyFortFunc* createfortfunc(PetscErrorCode (*func)(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode*),void *ctx,PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptr)) -static PetscErrorCode tdyfunc_wrapper(TDy tdy,PetscReal *x,PetscReal *f,void *ctx) -{ - TDyFortFunc* ffunc = (TDyFortFunc*)ctx; -#if defined(PETSC_HAVE_F90_2PTR_ARG) - void* ptr; - PetscObjectGetFortranCallback((PetscObject)tdy,PETSC_FORTRAN_CALLBACK_CLASS,ffunc->fort_pgiptr,NULL,&ptr); -#endif - PetscObjectUseFortranCallback(tdy,ffunc->fort_func,(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode* PETSC_F90_2PTR_PROTO_NOVAR),(&tdy,x,f,ffunc->fort_ctx,&ierr PETSC_F90_2PTR_PARAM(ptr))); -} - #if defined(__cplusplus) extern "C" { #endif @@ -199,16 +169,6 @@ PETSC_EXTERN void tdyfinalize_(int *__ierr){ } #endif -#if defined(__cplusplus) -extern "C" { -#endif -PETSC_EXTERN void tdyregisterpressurefn(char* name, PetscErrorCode (*func)(TDy*,PetscReal*,PetscReal*,void*,PetscErrorCode*),void *ctx,PetscErrorCode *ierr PETSC_F90_2PTR_PROTO(ptr)) -*__ierr = TDyRegisterPressureFn((TDy)PetscToPointer((tdy))); -} -#if defined(__cplusplus) -} -#endif - #if defined(__cplusplus) extern "C" { #endif diff --git a/src/tdyconditions.c b/src/tdyconditions.c index 8193c2f6..0e6defb0 100644 --- a/src/tdyconditions.c +++ b/src/tdyconditions.c @@ -5,46 +5,27 @@ Boundary and source-sink conditions are set by PETSc operations */ -// Here's a registry of pressure, temperature, and velocity functions. -typedef PetscErrorCode(*BCFunction)(TDy, PetscReal*, PetscReal*, void*); -KHASH_MAP_INIT_STR(TDY_BCFUNC_MAP, BCFunction) -static khash_t(TDY_BCFUNC_MAP)* pressure_funcs_ = NULL; -static khash_t(TDY_BCFUNC_MAP)* temperature_funcs_ = NULL; -static khash_t(TDY_BCFUNC_MAP)* velocity_funcs_ = NULL; - -PetscErrorCode TDyRegisterPressureFn(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { - PetscFunctionBegin; - if (pressure_funcs_ == NULL) { - pressure_funcs_ = kh_init(TDY_BCFUNC_MAP); - } - - int retval; - khiter_t iter = kh_put(TDY_BCFUNC_MAP, pressure_funcs_, name, &retval); - kh_val(pressure_funcs_, iter) = f; - PetscFunctionReturn(0); -} - -PetscErrorCode TDyRegisterTemperatureFn(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { - PetscFunctionBegin; - if (temperature_funcs_ == NULL) { - temperature_funcs_ = kh_init(TDY_BCFUNC_MAP); - } - - int retval; - khiter_t iter = kh_put(TDY_BCFUNC_MAP, temperature_funcs_, name, &retval); - kh_val(temperature_funcs_, iter) = f; - PetscFunctionReturn(0); +// Here's a registry of functions that can be used for boundary conditions and +// forcing terms. +typedef PetscErrorCode(*Function)(TDy, PetscReal*, PetscReal*, void*); +KHASH_MAP_INIT_STR(TDY_FUNC_MAP, Function) +static khash_t(TDY_FUNC_MAP)* funcs_ = NULL; + +// This function is called on finalization to destroy the function registry. +static void DestroyFunctionRegistry() { + kh_destroy(TDY_FUNC_MAP, funcs_); } -PetscErrorCode TDyRegisterVelocityFn(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { +PetscErrorCode TDyRegisterFunction(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { PetscFunctionBegin; - if (velocity_funcs_ == NULL) { - velocity_funcs_ = kh_init(TDY_BCFUNC_MAP); + if (funcs_ == NULL) { + funcs_ = kh_init(TDY_FUNC_MAP); + TDyOnFinalize(DestroyFunctionRegistry); } int retval; - khiter_t iter = kh_put(TDY_BCFUNC_MAP, velocity_funcs_, name, &retval); - kh_val(velocity_funcs_, iter) = f; + khiter_t iter = kh_put(TDY_FUNC_MAP, funcs_, name, &retval); + kh_val(funcs_, iter) = f; PetscFunctionReturn(0); } @@ -66,9 +47,9 @@ PetscErrorCode TDySelectBoundaryPressureFn(TDy tdy, const char* name, void* ctx) PetscFunctionBegin; int ierr; - khiter_t iter = kh_get(TDY_BCFUNC_MAP, pressure_funcs_, name); - if (iter != kh_end(pressure_funcs_)) { // found it! - BCFunction f = kh_val(pressure_funcs_, iter); + khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name); + if (iter != kh_end(funcs_)) { // found it! + Function f = kh_val(funcs_, iter); ierr = TDySetBoundaryPressureFn(tdy, f, ctx); CHKERRQ(ierr); } else { printf("Uh oh!\n"); // TODO: handle error condition @@ -80,9 +61,9 @@ PetscErrorCode TDySelectBoundaryTemperatureFn(TDy tdy, const char* name, void* c PetscFunctionBegin; int ierr; - khiter_t iter = kh_get(TDY_BCFUNC_MAP, temperature_funcs_, name); - if (iter != kh_end(temperature_funcs_)) { // found it! - BCFunction f = kh_val(temperature_funcs_, iter); + khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name); + if (iter != kh_end(funcs_)) { // found it! + Function f = kh_val(funcs_, iter); ierr = TDySetBoundaryTemperatureFn(tdy, f, ctx); CHKERRQ(ierr); } else { printf("Uh oh!\n"); // TODO: handle error condition @@ -94,9 +75,9 @@ PetscErrorCode TDySelectBoundaryVelocityFn(TDy tdy, const char* name, void* ctx) PetscFunctionBegin; int ierr; - khiter_t iter = kh_get(TDY_BCFUNC_MAP, velocity_funcs_, name); - if (iter != kh_end(velocity_funcs_)) { // found it! - BCFunction f = kh_val(velocity_funcs_, iter); + khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name); + if (iter != kh_end(funcs_)) { // found it! + Function f = kh_val(funcs_, iter); ierr = TDySetBoundaryVelocityFn(tdy, f, ctx); CHKERRQ(ierr); } else { printf("Uh oh!\n"); // TODO: handle error condition From 97e5ea10c457bce88dec111e1f5de8513e2a9807 Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 5 Aug 2021 09:21:35 -0700 Subject: [PATCH 6/8] Fixed remaining pointer-chasing problems. Works! --- src/f90-mod/tdycoremod.F90 | 28 +++++++++++++++++- src/interface/ftn/tdycoref.c | 6 ---- src/tdyconditions.c | 57 ++++++++++++++++++++---------------- 3 files changed, 59 insertions(+), 32 deletions(-) diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index af49edce..4c4d8fd8 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -504,7 +504,33 @@ function GetFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) ierr = GetFn(FtoCString(name), c_func) call c_f_procpointer(c_func, f_func) - call TDySetBoundaryPressureFn(tdy, f_func, c_null_ptr) + call TDySetBoundaryPressureFn(tdy, f_func, 0, ierr) + end subroutine + + subroutine TDySelectBoundaryVelocityFn(tdy, name, ierr) + use, intrinsic :: iso_c_binding + use tdycoredef + implicit none + TDy :: tdy + character(len=*), intent(in) :: name + PetscErrorCode :: ierr + + type(c_funptr) :: c_func + procedure(TDyFunction), pointer :: f_func + + interface + function GetFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: name + type(c_funptr) :: c_func + integer(c_int) :: ierr + end function + end interface + + ierr = GetFn(FtoCString(name), c_func) + call c_f_procpointer(c_func, f_func) + call TDySetBoundaryVelocityFn(tdy, f_func, 0, ierr) end subroutine ! Here's a function that converts a Fortran string to a C string and diff --git a/src/interface/ftn/tdycoref.c b/src/interface/ftn/tdycoref.c index 38efc9cd..79254367 100644 --- a/src/interface/ftn/tdycoref.c +++ b/src/interface/ftn/tdycoref.c @@ -39,9 +39,6 @@ #define tdysetpermeabilityfunction_ TDYSETPERMEABILITYFUNCTION #define tdysetresidualsaturationfunction_ TDYSETRESIDUALSATURATIONFUNCTION #define tdysetforcingfunction_ TDYSETFORCINGFUNCTION -#define tdyselectboundarypressurefn_ TDYSELECTBOUNDARYPRESSUREFN -#define tdyselectboundarytemperaturefn_ TDYSELECTBOUNDARYTEMPERATUREFN -#define tdyselectboundaryvelocityfn_ TDYSELECTBOUNDARYVELOCITYFN #define tdysetboundarypressurefn_ TDYSETBOUNDARYPRESSUREFN #define tdysetboundaryvelocityfn_ TDYSETBOUNDARYVELOCITYFN #define tdysetporosityvalueslocal0_ TDYSETPOROSITYVALUESLOCAL0 @@ -104,9 +101,6 @@ #define tdysetpermeabilityfunction_ tdysetpermeabilityfunction #define tdysetresidualsaturationfunction_ tdysetresidualsaturationfunction #define tdysetforcingfunction_ tdysetforcingfunction -#define tdyselectboundarypressurefn_ tdyselectboundarypressurefn -#define tdyselectboundarytemperaturefn_ tdyselectboundarytemperaturefn -#define tdyselectboundaryvelocityfn_ tdyselectboundaryvelocityfn #define tdysetboundarypressurefn_ tdysetboundarypressurefn #define tdysetboundaryvelocityfn_ tdysetboundaryvelocityfn #define tdysetporosityvalueslocal0_ tdysetporosityvalueslocal0 diff --git a/src/tdyconditions.c b/src/tdyconditions.c index 0e6defb0..95c43009 100644 --- a/src/tdyconditions.c +++ b/src/tdyconditions.c @@ -16,7 +16,7 @@ static void DestroyFunctionRegistry() { kh_destroy(TDY_FUNC_MAP, funcs_); } -PetscErrorCode TDyRegisterFunction(const char* name, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*)) { +PetscErrorCode TDyRegisterFunction(const char* name, Function f) { PetscFunctionBegin; if (funcs_ == NULL) { funcs_ = kh_init(TDY_FUNC_MAP); @@ -29,6 +29,28 @@ PetscErrorCode TDyRegisterFunction(const char* name, PetscErrorCode(*f)(TDy,Pets PetscFunctionReturn(0); } +PetscErrorCode TDyGetFunction(const char* name, Function* f) { + PetscFunctionBegin; + int ierr; + + if (funcs_ != NULL) { + khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name); + if (iter != kh_end(funcs_)) { // found it! + *f = kh_val(funcs_, iter); + } else { + ierr = -1; + SETERRQ(MPI_COMM_WORLD, ierr, "Function not found!"); + return ierr; + } + } else { + ierr = -1; + SETERRQ(MPI_COMM_WORLD, ierr, "No functions have been registered!"); + return ierr; + } + PetscFunctionReturn(0); +} + + PetscErrorCode TDySetForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,PetscReal*,PetscReal*,void*),void *ctx) { PetscFunctionBegin; if (f) tdy->ops->computeforcing = f; @@ -46,42 +68,27 @@ PetscErrorCode TDySetEnergyForcingFunction(TDy tdy, PetscErrorCode(*f)(TDy,Petsc PetscErrorCode TDySelectBoundaryPressureFn(TDy tdy, const char* name, void* ctx) { PetscFunctionBegin; int ierr; - - khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name); - if (iter != kh_end(funcs_)) { // found it! - Function f = kh_val(funcs_, iter); - ierr = TDySetBoundaryPressureFn(tdy, f, ctx); CHKERRQ(ierr); - } else { - printf("Uh oh!\n"); // TODO: handle error condition - } + Function f; + ierr = TDyGetFunction(name, &f); CHKERRQ(ierr); + ierr = TDySetBoundaryPressureFn(tdy, f, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode TDySelectBoundaryTemperatureFn(TDy tdy, const char* name, void* ctx) { PetscFunctionBegin; int ierr; - - khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name); - if (iter != kh_end(funcs_)) { // found it! - Function f = kh_val(funcs_, iter); - ierr = TDySetBoundaryTemperatureFn(tdy, f, ctx); CHKERRQ(ierr); - } else { - printf("Uh oh!\n"); // TODO: handle error condition - } + Function f; + ierr = TDyGetFunction(name, &f); CHKERRQ(ierr); + ierr = TDySetBoundaryTemperatureFn(tdy, f, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } PetscErrorCode TDySelectBoundaryVelocityFn(TDy tdy, const char* name, void* ctx) { PetscFunctionBegin; int ierr; - - khiter_t iter = kh_get(TDY_FUNC_MAP, funcs_, name); - if (iter != kh_end(funcs_)) { // found it! - Function f = kh_val(funcs_, iter); - ierr = TDySetBoundaryVelocityFn(tdy, f, ctx); CHKERRQ(ierr); - } else { - printf("Uh oh!\n"); // TODO: handle error condition - } + Function f; + ierr = TDyGetFunction(name, &f); CHKERRQ(ierr); + ierr = TDySetBoundaryVelocityFn(tdy, f, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } From 6e702a025bd7dfca0d373f104164311649df96ca Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 5 Aug 2021 09:31:42 -0700 Subject: [PATCH 7/8] Reversed the order in which shutdown functions are called. --- src/tdycore.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tdycore.c b/src/tdycore.c index 3620109e..f8926a3f 100644 --- a/src/tdycore.c +++ b/src/tdycore.c @@ -165,9 +165,9 @@ PetscErrorCode TDyOnFinalize(void (*shutdown_func)(void)) { PetscErrorCode TDyFinalize() { PetscFunctionBegin; - // Call shutdown functions and destroy the list. + // Call shutdown functions in reverse order, and destroy the list. if (shutdown_funcs_ != NULL) { - for (int i = 0; i < num_shutdown_funcs_; ++i) { + for (int i = num_shutdown_funcs_-1; i >= 0; --i) { shutdown_funcs_[i](); } free(shutdown_funcs_); From 130a002f8ec4664d2711373075124c343c9221fb Mon Sep 17 00:00:00 2001 From: "Jeffrey N. Johnson" Date: Thu, 5 Aug 2021 10:45:23 -0700 Subject: [PATCH 8/8] Added support for porosity and forcing functions to registry. --- demo/transient/transient_snes_mpfaof90.F90 | 20 ++-- src/f90-mod/tdycoremod.F90 | 101 ++++++++++++++++----- 2 files changed, 91 insertions(+), 30 deletions(-) diff --git a/demo/transient/transient_snes_mpfaof90.F90 b/demo/transient/transient_snes_mpfaof90.F90 index a61d7b60..8ae7fc31 100644 --- a/demo/transient/transient_snes_mpfaof90.F90 +++ b/demo/transient/transient_snes_mpfaof90.F90 @@ -9,8 +9,8 @@ module snes_mpfaof90mod subroutine PorosityFunction(tdy,x,theta,dummy,ierr) implicit none TDy :: tdy - PetscReal, intent(in) :: x - PetscReal, intent(out):: theta + PetscReal, intent(in) :: x(3) + PetscReal, intent(out) :: theta integer :: dummy(*) PetscErrorCode :: ierr @@ -49,8 +49,8 @@ end subroutine ResidualSaturation subroutine PorosityFunctionPFLOTRAN(tdy,x,theta,dummy,ierr) implicit none TDy :: tdy - PetscReal, intent(in) :: x - PetscReal, intent(out):: theta + PetscReal, intent(in) :: x(3) + PetscReal, intent(out) :: theta integer :: dummy(*) PetscErrorCode :: ierr @@ -100,8 +100,8 @@ end subroutine ResidualSat_PFLOTRAN subroutine PressureFunction(tdy,x,pressure,dummy,ierr) implicit none TDy :: tdy - PetscReal, intent(in) :: x(3) - PetscReal, intent(out):: pressure + PetscReal, intent(in) :: x(3) + PetscReal, intent(out) :: pressure integer :: dummy(*) PetscErrorCode :: ierr @@ -163,6 +163,8 @@ program main ! Register some functions. call TDyRegisterFunction("p0", PressureFunction, ierr) + call TDyRegisterFunction("porosity", PorosityFunction, ierr) + call TDyRegisterFunction("porosity_pflotran", PorosityFunctionPFLOTRAN, ierr) CHKERRA(ierr); call TDyCreate(tdy, ierr); @@ -310,7 +312,8 @@ program main CHKERRA(ierr); if (pflotran_consistent) then - call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr); +! call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,0,ierr); + call TDySelectPorosityFunction(tdy,"porosity_pflotran",ierr); CHKERRA(ierr); do c = 1,ncell @@ -333,7 +336,8 @@ program main else - call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr); +! call TDySetPorosityFunction(tdy,PorosityFunction,0,ierr); + call TDySelectPorosityFunction(tdy,"porosity",ierr); CHKERRA(ierr); call TDySetBlockPermeabilityValuesLocal(tdy,ncell,index,blockPerm,ierr); diff --git a/src/f90-mod/tdycoremod.F90 b/src/f90-mod/tdycoremod.F90 index 4c4d8fd8..48bd6d3c 100644 --- a/src/f90-mod/tdycoremod.F90 +++ b/src/f90-mod/tdycoremod.F90 @@ -448,6 +448,17 @@ subroutine TDyFunction(tdy, x, f, dummy, ierr) end subroutine end interface + ! We use GetRegFn to retrieve function pointers from the C registry. + interface + function GetRegFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) + use, intrinsic :: iso_c_binding + implicit none + type(c_ptr), value :: name + type(c_funptr) :: c_func + integer(c_int) :: ierr + end function + end interface + contains subroutine TDyInit(ierr) @@ -481,7 +492,7 @@ function RegisterFn(name, func) bind (c, name="TDyRegisterFunction") result(ierr ierr = RegisterFn(FtoCString(name), c_funloc(func)) end subroutine - subroutine TDySelectBoundaryPressureFn(tdy, name, ierr) + subroutine TDySelectPorosityFunction(tdy, name, ierr) use, intrinsic :: iso_c_binding use tdycoredef implicit none @@ -492,17 +503,73 @@ subroutine TDySelectBoundaryPressureFn(tdy, name, ierr) type(c_funptr) :: c_func procedure(TDyFunction), pointer :: f_func - interface - function GetFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) - use, intrinsic :: iso_c_binding - implicit none - type(c_ptr), value :: name - type(c_funptr) :: c_func - integer(c_int) :: ierr - end function - end interface + ierr = GetRegFn(FtoCString(name), c_func) + call c_f_procpointer(c_func, f_func) + call TDySetPorosityFunction(tdy, f_func, 0, ierr) + end subroutine + + subroutine TDySelectForcingFunction(tdy, name, ierr) + use, intrinsic :: iso_c_binding + use tdycoredef + implicit none + TDy :: tdy + character(len=*), intent(in) :: name + PetscErrorCode :: ierr + + type(c_funptr) :: c_func + procedure(TDyFunction), pointer :: f_func - ierr = GetFn(FtoCString(name), c_func) + ierr = GetRegFn(FtoCString(name), c_func) + call c_f_procpointer(c_func, f_func) + call TDySetForcingFunction(tdy, f_func, 0, ierr) + end subroutine + +! Uncomment this when we are ready to set energy forcing fns in Fortran. +! subroutine TDySelectEnergyForcingFunction(tdy, name, ierr) +! use, intrinsic :: iso_c_binding +! use tdycoredef +! implicit none +! TDy :: tdy +! character(len=*), intent(in) :: name +! PetscErrorCode :: ierr +! +! type(c_funptr) :: c_func +! procedure(TDyFunction), pointer :: f_func +! +! ierr = GetRegFn(FtoCString(name), c_func) +! call c_f_procpointer(c_func, f_func) +! call TDySetEnergyForcingFunction(tdy, f_func, 0, ierr) +! end subroutine + +! Uncomment this when we are ready to set permeability functions programmatically. +! subroutine TDySelectPermeabilityFunction(tdy, name, ierr) +! use, intrinsic :: iso_c_binding +! use tdycoredef +! implicit none +! TDy :: tdy +! character(len=*), intent(in) :: name +! PetscErrorCode :: ierr +! +! type(c_funptr) :: c_func +! procedure(TDyFunction), pointer :: f_func +! +! ierr = GetRegFn(FtoCString(name), c_func) +! call c_f_procpointer(c_func, f_func) +! call TDySetPermeabilityFunction(tdy, f_func, 0, ierr) +! end subroutine + + subroutine TDySelectBoundaryPressureFn(tdy, name, ierr) + use, intrinsic :: iso_c_binding + use tdycoredef + implicit none + TDy :: tdy + character(len=*), intent(in) :: name + PetscErrorCode :: ierr + + type(c_funptr) :: c_func + procedure(TDyFunction), pointer :: f_func + + ierr = GetRegFn(FtoCString(name), c_func) call c_f_procpointer(c_func, f_func) call TDySetBoundaryPressureFn(tdy, f_func, 0, ierr) end subroutine @@ -518,17 +585,7 @@ subroutine TDySelectBoundaryVelocityFn(tdy, name, ierr) type(c_funptr) :: c_func procedure(TDyFunction), pointer :: f_func - interface - function GetFn(name, c_func) bind (c, name="TDyGetFunction") result(ierr) - use, intrinsic :: iso_c_binding - implicit none - type(c_ptr), value :: name - type(c_funptr) :: c_func - integer(c_int) :: ierr - end function - end interface - - ierr = GetFn(FtoCString(name), c_func) + ierr = GetRegFn(FtoCString(name), c_func) call c_f_procpointer(c_func, f_func) call TDySetBoundaryVelocityFn(tdy, f_func, 0, ierr) end subroutine