@@ -5,12 +5,13 @@ module musica_micm
55#define ASSERT( expr ) call assert( expr, __FILE__, __LINE__ )
66 use iso_c_binding, only: c_ptr, c_char, c_int, c_int64_t, c_bool, c_double, c_null_char, &
77 c_size_t, c_f_pointer, c_funptr, c_null_ptr, c_associated
8- use iso_fortran_env, only: int64
8+ use iso_fortran_env, only: int64, real64
99 use musica_state, only: state_t
1010 use musica_util, only: assert, mappings_t, string_t, string_t_c, error_t_c
1111 implicit none
1212
1313 public :: micm_t, solver_stats_t, get_micm_version, is_cuda_available
14+ public :: rosenbrock_solver_parameters_t, backward_euler_solver_parameters_t
1415 public :: UndefinedSolver, Rosenbrock, RosenbrockStandardOrder, BackwardEuler, BackwardEulerStandardOrder, CudaRosenbrock
1516 private
1617
@@ -37,6 +38,46 @@ module musica_micm
3738 enumerator :: CudaRosenbrock = 5
3839 end enum
3940
41+ ! > C-binding type for Rosenbrock solver parameters
42+ type, bind(c) :: rosenbrock_solver_parameters_t_c
43+ real (c_double) :: relative_tolerance
44+ type (c_ptr) :: absolute_tolerances = c_null_ptr
45+ integer (c_size_t) :: num_absolute_tolerances = 0
46+ real (c_double) :: h_min
47+ real (c_double) :: h_max
48+ real (c_double) :: h_start
49+ integer (c_size_t) :: max_number_of_steps
50+ end type rosenbrock_solver_parameters_t_c
51+
52+ ! > C-binding type for Backward Euler solver parameters
53+ type, bind(c) :: backward_euler_solver_parameters_t_c
54+ real (c_double) :: relative_tolerance
55+ type (c_ptr) :: absolute_tolerances = c_null_ptr
56+ integer (c_size_t) :: num_absolute_tolerances = 0
57+ integer (c_size_t) :: max_number_of_steps
58+ type (c_ptr) :: time_step_reductions = c_null_ptr
59+ integer (c_size_t) :: num_time_step_reductions = 0
60+ end type backward_euler_solver_parameters_t_c
61+
62+ ! > Fortran-native type for Rosenbrock solver parameters
63+ type :: rosenbrock_solver_parameters_t
64+ real (real64) :: relative_tolerance = 1.0e-6_real64
65+ real (real64), allocatable :: absolute_tolerances(:)
66+ real (real64) :: h_min = 0.0_real64
67+ real (real64) :: h_max = 0.0_real64
68+ real (real64) :: h_start = 0.0_real64
69+ integer :: max_number_of_steps = 1000
70+ end type rosenbrock_solver_parameters_t
71+
72+ ! > Fortran-native type for Backward Euler solver parameters
73+ type :: backward_euler_solver_parameters_t
74+ real (real64) :: relative_tolerance = 1.0e-6_real64
75+ real (real64), allocatable :: absolute_tolerances(:)
76+ integer :: max_number_of_steps = 11
77+ real (real64) :: time_step_reductions(5 ) = & ! must have exactly 5 elements
78+ (/ 0.5_real64 , 0.5_real64 , 0.5_real64 , 0.5_real64 , 0.1_real64 / )
79+ end type backward_euler_solver_parameters_t
80+
4081 interface
4182 function create_micm_c (config_path , solver_type , error ) &
4283 bind(C, name= " CreateMicm" )
@@ -125,6 +166,42 @@ integer(c_size_t) function get_maximum_number_of_grid_cells_c(micm) &
125166 import c_ptr, c_size_t
126167 type (c_ptr), value, intent (in ) :: micm
127168 end function get_maximum_number_of_grid_cells_c
169+
170+ subroutine set_rosenbrock_solver_parameters_c (micm , params , error ) &
171+ bind(C, name= " SetRosenbrockSolverParameters" )
172+ use musica_util, only: error_t_c
173+ import c_ptr, rosenbrock_solver_parameters_t_c
174+ type (c_ptr), value, intent (in ) :: micm
175+ type (rosenbrock_solver_parameters_t_c), intent (in ) :: params
176+ type (error_t_c), intent (inout ) :: error
177+ end subroutine set_rosenbrock_solver_parameters_c
178+
179+ subroutine set_backward_euler_solver_parameters_c (micm , params , error ) &
180+ bind(C, name= " SetBackwardEulerSolverParameters" )
181+ use musica_util, only: error_t_c
182+ import c_ptr, backward_euler_solver_parameters_t_c
183+ type (c_ptr), value, intent (in ) :: micm
184+ type (backward_euler_solver_parameters_t_c), intent (in ) :: params
185+ type (error_t_c), intent (inout ) :: error
186+ end subroutine set_backward_euler_solver_parameters_c
187+
188+ subroutine get_rosenbrock_solver_parameters_c (micm , params , error ) &
189+ bind(C, name= " GetRosenbrockSolverParameters" )
190+ use musica_util, only: error_t_c
191+ import c_ptr, rosenbrock_solver_parameters_t_c
192+ type (c_ptr), value, intent (in ) :: micm
193+ type (rosenbrock_solver_parameters_t_c), intent (inout ) :: params
194+ type (error_t_c), intent (inout ) :: error
195+ end subroutine get_rosenbrock_solver_parameters_c
196+
197+ subroutine get_backward_euler_solver_parameters_c (micm , params , error ) &
198+ bind(C, name= " GetBackwardEulerSolverParameters" )
199+ use musica_util, only: error_t_c
200+ import c_ptr, backward_euler_solver_parameters_t_c
201+ type (c_ptr), value, intent (in ) :: micm
202+ type (backward_euler_solver_parameters_t_c), intent (inout ) :: params
203+ type (error_t_c), intent (inout ) :: error
204+ end subroutine get_backward_euler_solver_parameters_c
128205 end interface
129206
130207 type :: micm_t
@@ -140,6 +217,11 @@ end function get_maximum_number_of_grid_cells_c
140217 procedure :: get_species_property_double
141218 procedure :: get_species_property_int
142219 procedure :: get_species_property_bool
220+ ! Solver parameters
221+ procedure :: set_rosenbrock_solver_parameters
222+ procedure :: set_backward_euler_solver_parameters
223+ procedure :: get_rosenbrock_solver_parameters
224+ procedure :: get_backward_euler_solver_parameters
143225 ! Deallocate the micm instance
144226 final :: finalize
145227 end type micm_t
@@ -434,5 +516,123 @@ function is_cuda_available(error) result(value)
434516 value = is_cuda_available_c(error_c)
435517 error = error_t(error_c)
436518 end function is_cuda_available
437-
519+
520+ subroutine set_rosenbrock_solver_parameters (this , params , error )
521+ use musica_util, only: error_t_c, error_t
522+ use iso_c_binding, only: c_loc
523+ class(micm_t), intent (in ) :: this
524+ type (rosenbrock_solver_parameters_t), intent (in ), target :: params
525+ type (error_t), intent (inout ) :: error
526+
527+ type (rosenbrock_solver_parameters_t_c) :: params_c
528+ type (error_t_c) :: error_c
529+ real (c_double), allocatable , target :: abs_tol_c(:)
530+
531+ params_c% relative_tolerance = real (params% relative_tolerance, c_double)
532+ params_c% h_min = real (params% h_min, c_double)
533+ params_c% h_max = real (params% h_max, c_double)
534+ params_c% h_start = real (params% h_start, c_double)
535+ params_c% max_number_of_steps = int (params% max_number_of_steps, c_size_t)
536+ if (allocated (params% absolute_tolerances)) then
537+ allocate (abs_tol_c(size (params% absolute_tolerances)))
538+ abs_tol_c = real (params% absolute_tolerances, c_double)
539+ params_c% absolute_tolerances = c_loc(abs_tol_c(1 ))
540+ params_c% num_absolute_tolerances = int (size (abs_tol_c), c_size_t)
541+ else
542+ params_c% absolute_tolerances = c_null_ptr
543+ params_c% num_absolute_tolerances = 0
544+ end if
545+ call set_rosenbrock_solver_parameters_c(this% ptr, params_c, error_c)
546+ error = error_t(error_c)
547+ end subroutine set_rosenbrock_solver_parameters
548+
549+ subroutine set_backward_euler_solver_parameters (this , params , error )
550+ use musica_util, only: error_t_c, error_t
551+ use iso_c_binding, only: c_loc
552+ class(micm_t), intent (in ) :: this
553+ type (backward_euler_solver_parameters_t), intent (in ), target :: params
554+ type (error_t), intent (inout ) :: error
555+
556+ type (backward_euler_solver_parameters_t_c) :: params_c
557+ type (error_t_c) :: error_c
558+ real (c_double), allocatable , target :: abs_tol_c(:)
559+ real (c_double), allocatable , target :: tsr_c(:)
560+
561+ params_c% relative_tolerance = real (params% relative_tolerance, c_double)
562+ params_c% max_number_of_steps = int (params% max_number_of_steps, c_size_t)
563+ if (allocated (params% absolute_tolerances)) then
564+ allocate (abs_tol_c(size (params% absolute_tolerances)))
565+ abs_tol_c = real (params% absolute_tolerances, c_double)
566+ params_c% absolute_tolerances = c_loc(abs_tol_c(1 ))
567+ params_c% num_absolute_tolerances = int (size (abs_tol_c), c_size_t)
568+ else
569+ params_c% absolute_tolerances = c_null_ptr
570+ params_c% num_absolute_tolerances = 0
571+ end if
572+ allocate (tsr_c(size (params% time_step_reductions)))
573+ tsr_c = real (params% time_step_reductions, c_double)
574+ params_c% time_step_reductions = c_loc(tsr_c(1 ))
575+ params_c% num_time_step_reductions = int (size (tsr_c), c_size_t)
576+ call set_backward_euler_solver_parameters_c(this% ptr, params_c, error_c)
577+ error = error_t(error_c)
578+ end subroutine set_backward_euler_solver_parameters
579+
580+ function get_rosenbrock_solver_parameters (this , error ) result(params)
581+ use musica_util, only: error_t_c, error_t
582+ class(micm_t), intent (in ) :: this
583+ type (error_t), intent (inout ) :: error
584+ type (rosenbrock_solver_parameters_t) :: params
585+
586+ type (rosenbrock_solver_parameters_t_c) :: params_c
587+ type (error_t_c) :: error_c
588+ real (c_double), pointer :: abs_tol_ptr(:)
589+
590+ params_c% absolute_tolerances = c_null_ptr
591+ params_c% num_absolute_tolerances = 0
592+ call get_rosenbrock_solver_parameters_c(this% ptr, params_c, error_c)
593+ error = error_t(error_c)
594+ if (.not. error% is_success()) return
595+ params% relative_tolerance = real (params_c% relative_tolerance, real64)
596+ params% h_min = real (params_c% h_min, real64)
597+ params% h_max = real (params_c% h_max, real64)
598+ params% h_start = real (params_c% h_start, real64)
599+ params% max_number_of_steps = int (params_c% max_number_of_steps)
600+ if (params_c% num_absolute_tolerances > 0 .and. c_associated(params_c% absolute_tolerances)) then
601+ call c_f_pointer(params_c% absolute_tolerances, abs_tol_ptr, [int (params_c% num_absolute_tolerances)])
602+ allocate (params% absolute_tolerances(int (params_c% num_absolute_tolerances)))
603+ params% absolute_tolerances = real (abs_tol_ptr, real64)
604+ end if
605+ end function get_rosenbrock_solver_parameters
606+
607+ function get_backward_euler_solver_parameters (this , error ) result(params)
608+ use musica_util, only: error_t_c, error_t
609+ class(micm_t), intent (in ) :: this
610+ type (error_t), intent (inout ) :: error
611+ type (backward_euler_solver_parameters_t) :: params
612+
613+ type (backward_euler_solver_parameters_t_c) :: params_c
614+ type (error_t_c) :: error_c
615+ real (c_double), pointer :: abs_tol_ptr(:)
616+ real (c_double), pointer :: tsr_ptr(:)
617+
618+ params_c% absolute_tolerances = c_null_ptr
619+ params_c% num_absolute_tolerances = 0
620+ params_c% time_step_reductions = c_null_ptr
621+ params_c% num_time_step_reductions = 0
622+ call get_backward_euler_solver_parameters_c(this% ptr, params_c, error_c)
623+ error = error_t(error_c)
624+ if (.not. error% is_success()) return
625+ params% relative_tolerance = real (params_c% relative_tolerance, real64)
626+ params% max_number_of_steps = int (params_c% max_number_of_steps)
627+ if (params_c% num_absolute_tolerances > 0 .and. c_associated(params_c% absolute_tolerances)) then
628+ call c_f_pointer(params_c% absolute_tolerances, abs_tol_ptr, [int (params_c% num_absolute_tolerances)])
629+ allocate (params% absolute_tolerances(int (params_c% num_absolute_tolerances)))
630+ params% absolute_tolerances = real (abs_tol_ptr, real64)
631+ end if
632+ if (params_c% num_time_step_reductions > 0 .and. c_associated(params_c% time_step_reductions)) then
633+ call c_f_pointer(params_c% time_step_reductions, tsr_ptr, [int (params_c% num_time_step_reductions)])
634+ params% time_step_reductions = real (tsr_ptr, real64)
635+ end if
636+ end function get_backward_euler_solver_parameters
637+
438638end module musica_micm
0 commit comments