OpenCMISS-Iron Internal API Documentation
solver_routines.f90
Go to the documentation of this file.
1 
43 
46 
47  USE base_routines
49 #ifdef WITH_CELLML
50  USE cellml_model_definition
51 #endif
52  USE cmiss_cellml
53  USE cmisspetsc
54  USE cmisspetsctypes
56  USE constants
59  USE field_routines
60  USE kinds
61  USE input_output
65  USE maths
69  USE strings
70  USE timer
71  USE types
72 
73 #include "macros.h"
74 
75  IMPLICIT NONE
76 
77  PRIVATE
78 
79 #include "petscversion.h"
80 
81  !Module parameters
82 
87  INTEGER(INTG), PARAMETER :: solver_number_of_solver_types=9
88  INTEGER(INTG), PARAMETER :: solver_linear_type=1
89  INTEGER(INTG), PARAMETER :: solver_nonlinear_type=2
90  INTEGER(INTG), PARAMETER :: solver_dynamic_type=3
91  INTEGER(INTG), PARAMETER :: solver_dae_type=4
92  INTEGER(INTG), PARAMETER :: solver_eigenproblem_type=5
93  INTEGER(INTG), PARAMETER :: solver_optimiser_type=6
94  INTEGER(INTG), PARAMETER :: solver_cellml_evaluator_type=7
95  INTEGER(INTG), PARAMETER :: solver_state_iteration_type=8
96  INTEGER(INTG), PARAMETER :: solver_geometric_transformation_type=9
98 
103  INTEGER(INTG), PARAMETER :: solver_cmiss_library=library_cmiss_type
104  INTEGER(INTG), PARAMETER :: solver_petsc_library=library_petsc_type
105  INTEGER(INTG), PARAMETER :: solver_mumps_library=library_mumps_type
106  INTEGER(INTG), PARAMETER :: solver_superlu_library=library_superlu_type
107  INTEGER(INTG), PARAMETER :: solver_spooles_library=library_spooles_type
108  INTEGER(INTG), PARAMETER :: solver_umfpack_library=library_umfpack_type
109  INTEGER(INTG), PARAMETER :: solver_lusol_library=library_lusol_type
110  INTEGER(INTG), PARAMETER :: solver_essl_library=library_essl_type
111  INTEGER(INTG), PARAMETER :: solver_lapack_library=library_lapack_type
112  INTEGER(INTG), PARAMETER :: solver_tao_library=library_tao_type
113  INTEGER(INTG), PARAMETER :: solver_hypre_library=library_hypre_type
114  INTEGER(INTG), PARAMETER :: solver_pastix_library=library_pastix_type
116 
121  INTEGER(INTG), PARAMETER :: solver_linear_direct_solve_type=1
122  INTEGER(INTG), PARAMETER :: solver_linear_iterative_solve_type=2
124 
129  INTEGER(INTG), PARAMETER :: solver_direct_lu=1
130  INTEGER(INTG), PARAMETER :: solver_direct_cholesky=2
131  INTEGER(INTG), PARAMETER :: solver_direct_svd=3
133 
138  INTEGER(INTG), PARAMETER :: solver_iterative_richardson=1
139  INTEGER(INTG), PARAMETER :: solver_iterative_chebyshev=2
140  INTEGER(INTG), PARAMETER :: solver_iterative_conjugate_gradient=3
141  INTEGER(INTG), PARAMETER :: solver_iterative_biconjugate_gradient=4
142  INTEGER(INTG), PARAMETER :: solver_iterative_gmres=5
143  INTEGER(INTG), PARAMETER :: solver_iterative_bicgstab=6
144  INTEGER(INTG), PARAMETER :: solver_iterative_conjgrad_squared=7
146 
151  INTEGER(INTG), PARAMETER :: solver_iterative_no_preconditioner=0
152  INTEGER(INTG), PARAMETER :: solver_iterative_jacobi_preconditioner=1
153  INTEGER(INTG), PARAMETER :: solver_iterative_block_jacobi_preconditioner=2
154  INTEGER(INTG), PARAMETER :: solver_iterative_sor_preconditioner=3
156  INTEGER(INTG), PARAMETER :: solver_iterative_incomplete_lu_preconditioner=5
159 
164  INTEGER(INTG), PARAMETER :: solver_nonlinear_newton=1
165  INTEGER(INTG), PARAMETER :: solver_nonlinear_bfgs_inverse=2
166  INTEGER(INTG), PARAMETER :: solver_nonlinear_sqp=3
167  INTEGER(INTG), PARAMETER :: solver_nonlinear_quasi_newton=4
169 
174  INTEGER(INTG), PARAMETER :: solver_quasi_newton_linesearch=1
175  INTEGER(INTG), PARAMETER :: solver_quasi_newton_trustregion=2
177 
182  INTEGER(INTG), PARAMETER :: solver_quasi_newton_lbfgs=1
183  INTEGER(INTG), PARAMETER :: solver_quasi_newton_goodbroyden=2
184  INTEGER(INTG), PARAMETER :: solver_quasi_newton_badbroyden=3
186 
191  INTEGER(INTG), PARAMETER :: solver_quasi_newton_linesearch_basic=1
192  INTEGER(INTG), PARAMETER :: solver_quasi_newton_linesearch_l2=2
193  INTEGER(INTG), PARAMETER :: solver_quasi_newton_linesearch_cp=3
195 
200  INTEGER(INTG), PARAMETER :: solver_quasi_newton_restart_none=1
201  INTEGER(INTG), PARAMETER :: solver_quasi_newton_restart_powell=2
202  INTEGER(INTG), PARAMETER :: solver_quasi_newton_restart_periodic=3
204 
209  INTEGER(INTG), PARAMETER :: solver_quasi_newton_scale_none=1
210  INTEGER(INTG), PARAMETER :: solver_quasi_newton_scale_shanno=2
211  INTEGER(INTG), PARAMETER :: solver_quasi_newton_scale_linesearch=3
212  INTEGER(INTG), PARAMETER :: solver_quasi_newton_scale_jacobian=4
214 
219  INTEGER(INTG), PARAMETER :: solver_newton_linesearch=1
220  INTEGER(INTG), PARAMETER :: solver_newton_trustregion=2
222 
227  INTEGER(INTG), PARAMETER :: solver_newton_linesearch_nonorms=1
228  INTEGER(INTG), PARAMETER :: solver_newton_linesearch_linear=2
229  INTEGER(INTG), PARAMETER :: solver_newton_linesearch_quadratic=3
230  INTEGER(INTG), PARAMETER :: solver_newton_linesearch_cubic=4
232 
237  INTEGER(INTG), PARAMETER :: solver_newton_jacobian_not_calculated=1
238  INTEGER(INTG), PARAMETER :: solver_newton_jacobian_equations_calculated=2
239  INTEGER(INTG), PARAMETER :: solver_newton_jacobian_fd_calculated=3
241 
246  INTEGER(INTG), PARAMETER :: solver_newton_convergence_petsc_default=1
247  INTEGER(INTG), PARAMETER :: solver_newton_convergence_energy_norm=2
250 
255  INTEGER(INTG), PARAMETER :: solver_dynamic_first_order=1
256  INTEGER(INTG), PARAMETER :: solver_dynamic_second_order=2
258 
263  INTEGER(INTG), PARAMETER :: solver_dynamic_linear=1
264  INTEGER(INTG), PARAMETER :: solver_dynamic_nonlinear=2
266 
271  INTEGER(INTG), PARAMETER :: solver_dynamic_first_degree=1
272  INTEGER(INTG), PARAMETER :: solver_dynamic_second_degree=2
273  INTEGER(INTG), PARAMETER :: solver_dynamic_third_degree=3
275 
280  INTEGER(INTG), PARAMETER :: solver_dynamic_euler_scheme=1
281  INTEGER(INTG), PARAMETER :: solver_dynamic_backward_euler_scheme=2
282  INTEGER(INTG), PARAMETER :: solver_dynamic_crank_nicolson_scheme=3
283  INTEGER(INTG), PARAMETER :: solver_dynamic_galerkin_scheme=4
284  INTEGER(INTG), PARAMETER :: solver_dynamic_zlamal_scheme=5
285  INTEGER(INTG), PARAMETER :: solver_dynamic_second_degree_gear_scheme=6
286  INTEGER(INTG), PARAMETER :: solver_dynamic_second_degree_liniger1_scheme=7
287  INTEGER(INTG), PARAMETER :: solver_dynamic_second_degree_liniger2_scheme=8
288  INTEGER(INTG), PARAMETER :: solver_dynamic_newmark1_scheme=9
289  INTEGER(INTG), PARAMETER :: solver_dynamic_newmark2_scheme=10
290  INTEGER(INTG), PARAMETER :: solver_dynamic_newmark3_scheme=11
291  INTEGER(INTG), PARAMETER :: solver_dynamic_third_degree_gear_scheme=12
292  INTEGER(INTG), PARAMETER :: solver_dynamic_third_degree_liniger1_scheme=13
293  INTEGER(INTG), PARAMETER :: solver_dynamic_third_degree_liniger2_scheme=14
294  INTEGER(INTG), PARAMETER :: solver_dynamic_houbolt_scheme=15
295  INTEGER(INTG), PARAMETER :: solver_dynamic_wilson_scheme=16
296  INTEGER(INTG), PARAMETER :: solver_dynamic_bossak_newmark1_scheme=17
297  INTEGER(INTG), PARAMETER :: solver_dynamic_bossak_newmark2_scheme=18
298  INTEGER(INTG), PARAMETER :: solver_dynamic_hilbert_hughes_taylor1_scheme=19
299  INTEGER(INTG), PARAMETER :: solver_dynamic_hilbert_hughes_taylor2_scheme=20
300  INTEGER(INTG), PARAMETER :: solver_dynamic_user_defined_scheme=21
302 
307  INTEGER(INTG), PARAMETER :: solver_dae_differential_only=0
308  INTEGER(INTG), PARAMETER :: solver_dae_index_1=1
309  INTEGER(INTG), PARAMETER :: solver_dae_index_2=2
310  INTEGER(INTG), PARAMETER :: solver_dae_index_3=3
312 
317  INTEGER(INTG), PARAMETER :: solver_dae_euler=1
318  INTEGER(INTG), PARAMETER :: solver_dae_crank_nicolson=2
319  INTEGER(INTG), PARAMETER :: solver_dae_runge_kutta=3
320  INTEGER(INTG), PARAMETER :: solver_dae_adams_moulton=4
321  INTEGER(INTG), PARAMETER :: solver_dae_bdf=5
322  INTEGER(INTG), PARAMETER :: solver_dae_rush_larson=6
323  INTEGER(INTG), PARAMETER :: solver_dae_external=7
324 
326 
331  INTEGER(INTG), PARAMETER :: solver_dae_euler_forward=1
332  INTEGER(INTG), PARAMETER :: solver_dae_euler_backward=2
333  INTEGER(INTG), PARAMETER :: solver_dae_euler_improved=3
335 
340  INTEGER(INTG), PARAMETER :: solver_solution_initialise_zero=0
341  INTEGER(INTG), PARAMETER :: solver_solution_initialise_current_field=1
342  INTEGER(INTG), PARAMETER :: solver_solution_initialise_no_change=2
344 
349  INTEGER(INTG), PARAMETER :: solver_no_output=0
350  INTEGER(INTG), PARAMETER :: solver_progress_output=1
351  INTEGER(INTG), PARAMETER :: solver_timing_output=2
352  INTEGER(INTG), PARAMETER :: solver_solver_output=3
353  INTEGER(INTG), PARAMETER :: solver_matrix_output=4
355 
360  INTEGER(INTG), PARAMETER :: solver_sparse_matrices=1
361  INTEGER(INTG), PARAMETER :: solver_full_matrices=2
363  !Module types
364 
365  !Module variables
366 
367  !Interfaces
368 
369  INTERFACE
370 
371  SUBROUTINE solver_dae_external_integrate(NUMBER_OF_DOFS,START_TIME,END_TIME,INITIAL_STEP, &
372  & only_one_model_index,models_data,number_of_state,state_data,number_of_parameters, &
373  & parameters_data,number_of_intermediate,intermediate_data,err) bind(c, name="SolverDAEExternalIntegrate")
374 
375  USE iso_c_binding
376 
377  INTEGER(C_INT), VALUE, INTENT(IN) :: NUMBER_OF_DOFS
378  REAL(C_DOUBLE), VALUE, INTENT(IN) :: START_TIME
379  REAL(C_DOUBLE), VALUE, INTENT(IN) :: END_TIME
380  REAL(C_DOUBLE), INTENT(INOUT) :: INITIAL_STEP
381  INTEGER(C_INT), VALUE, INTENT(IN) :: ONLY_ONE_MODEL_INDEX
382  INTEGER(C_INT), INTENT(IN) :: MODELS_DATA(*)
383  INTEGER(C_INT), VALUE, INTENT(IN) :: NUMBER_OF_STATE
384  REAL(C_DOUBLE), INTENT(INOUT) :: STATE_DATA(*)
385  INTEGER(C_INT), VALUE, INTENT(IN) :: NUMBER_OF_PARAMETERS
386  REAL(C_DOUBLE), INTENT(IN) :: PARAMETERS_DATA(*)
387  INTEGER(C_INT), VALUE, INTENT(IN) :: NUMBER_OF_INTERMEDIATE
388  REAL(C_DOUBLE), INTENT(OUT) :: INTERMEDIATE_DATA(*)
389  INTEGER(C_INT), INTENT(OUT) :: ERR
390 
391  END SUBROUTINE solver_dae_external_integrate
392 
393  END INTERFACE
394 
396  MODULE PROCEDURE solver_dynamic_theta_set_dp1
397  MODULE PROCEDURE solver_dynamic_theta_set_dp
398  END INTERFACE solver_dynamic_theta_set
399 
401  MODULE PROCEDURE solver_label_get_c
402  MODULE PROCEDURE solver_label_get_vs
403  END INTERFACE solver_label_get
404 
406  MODULE PROCEDURE solver_label_set_c
407  MODULE PROCEDURE solver_label_set_vs
408  END INTERFACE solver_label_set
409 
411 
414 
418 
420 
422 
425 
429 
431 
433 
436 
438 
441 
444 
447 
450 
453 
456 
458 
460 
462 
472 
474 
477 
479 
481 
483 
485 
487 
490 
492 
494 
496 
498 
500 
502 
504 
506 
507  PUBLIC solver_destroy
508 
510 
512 
514 
516 
518 
520 
522 
524 
526 
528 
530 
532 
534 
536 
538 
540 
542 
544 
546 
548 
550 
552 
554 
556 
558 
560 
562 
564 
566 
568 
570 
572 
574 
576 
578 
580 
582 
584 
586 
588 
590 
592 
594 
596 
598 
600 
602 
604 
606 
608 
610 
612 
614 
616 
618 
620 
622 
624 
626 
628 
630 
632 
634 
636 
638 
640 
642 
644 
646 
648 
650 
652 
654 
656 
658 
660 
662 
664 
666 
668 
670 
672 
674 
676 
678 
680 
682 
683  PUBLIC solver_solve
684 
686 
688 
689  PUBLIC solver_type_set
690 
692 
694 
696 
698 
700 
701  PUBLIC solvers_destroy
702 
703  PUBLIC solvers_number_set
704 
705  PUBLIC solvers_solver_get
706 
708 
710 
712 
713 CONTAINS
714 
715  !
716  !================================================================================================================================
717  !
718 
720  SUBROUTINE cellml_equations_cellml_add(CELLML_EQUATIONS,CELLML,CELLML_INDEX,ERR,ERROR,*)
722  !Argument variables
723  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
724  TYPE(cellml_type), POINTER :: CELLML
725  INTEGER(INTG), INTENT(OUT) :: CELLML_INDEX
726  INTEGER(INTG), INTENT(OUT) :: ERR
727  TYPE(varying_string), INTENT(OUT) :: ERROR
728  !Local Variables
729  INTEGER(INTG) :: cellml_idx
730  TYPE(cellml_ptr_type), ALLOCATABLE :: NEW_CELLML_ENVIRONMENTS(:)
731  TYPE(solver_type), POINTER :: SOLVER
732 
733  enters("CELLML_EQUATIONS_CELLML_ADD",err,error,*999)
734 
735  IF(ASSOCIATED(cellml_equations)) THEN
736  IF(cellml_equations%CELLML_EQUATIONS_FINISHED) THEN
737  CALL flagerror("CellML equations has already been finished.",err,error,*999)
738  ELSE
739  solver=>cellml_equations%SOLVER
740  IF(ASSOCIATED(solver)) THEN
741  IF(ASSOCIATED(cellml)) THEN
742  IF(cellml%CELLML_FINISHED) THEN
743  ALLOCATE(new_cellml_environments(cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS+1),stat=err)
744  IF(err/=0) CALL flagerror("Could not allocate new CellML environments.",err,error,*999)
745  DO cellml_idx=1,cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS
746  new_cellml_environments(cellml_idx)%PTR=>cellml_equations%CELLML_ENVIRONMENTS(cellml_idx)%PTR
747  ENDDO !cellml_idx
748  new_cellml_environments(cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS+1)%PTR=>cellml
749  CALL move_alloc(new_cellml_environments,cellml_equations%CELLML_ENVIRONMENTS)
750  cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS=cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS+1
751  cellml_index=cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS
752  ELSE
753  CALL flagerror("CellML environment has not been finished.",err,error,*999)
754  ENDIF
755  ELSE
756  CALL flagerror("CellML environment is not associated.",err,error,*999)
757  ENDIF
758  ELSE
759  CALL flagerror("CellML equations solver is not associated.",err,error,*999)
760  ENDIF
761  ENDIF
762  ELSE
763  CALL flagerror("CellML equations is not associated.",err,error,*999)
764  ENDIF
765 
766  exits("CELLML_EQUATIONS_CELLML_ADD")
767  RETURN
768 999 IF(ALLOCATED(new_cellml_environments)) DEALLOCATE(new_cellml_environments)
769  errorsexits("CELLML_EQUATIONS_CELLML_ADD",err,error)
770  RETURN 1
771 
772  END SUBROUTINE cellml_equations_cellml_add
773 
774  !
775  !================================================================================================================================
776  !
777 
779  SUBROUTINE cellml_equations_create_finish(CELLML_EQUATIONS,ERR,ERROR,*)
781  !Argument variables
782  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
783  INTEGER(INTG), INTENT(OUT) :: ERR
784  TYPE(varying_string), INTENT(OUT) :: ERROR
785  !Local Variables
786  TYPE(solver_type), POINTER :: SOLVER
787 
788  enters("CELLML_EQUATIONS_CREATE_FINISH",err,error,*999)
789 
790  IF(ASSOCIATED(cellml_equations)) THEN
791  IF(cellml_equations%CELLML_EQUATIONS_FINISHED) THEN
792  CALL flagerror("CellML equations has already been finished.",err,error,*999)
793  ELSE
794  solver=>cellml_equations%SOLVER
795  IF(ASSOCIATED(solver)) THEN
796  cellml_equations%CELLML_EQUATIONS_FINISHED=.true.
797  ELSE
798  CALL flagerror("CellML equations solver is not associated.",err,error,*999)
799  ENDIF
800  ENDIF
801  ELSE
802  CALL flagerror("CellML equations is not associated.",err,error,*999)
803  ENDIF
804 
805  exits("CELLML_EQUATIONS_CREATE_FINISH")
806  RETURN
807 999 errorsexits("CELLML_EQUATIONS_CREATE_FINISH",err,error)
808  RETURN 1
809 
810  END SUBROUTINE cellml_equations_create_finish
811 
812  !
813  !================================================================================================================================
814  !
815 
817  SUBROUTINE cellml_equations_create_start(SOLVER,CELLML_EQUATIONS,ERR,ERROR,*)
819  !Argument variables
820  TYPE(solver_type), POINTER :: SOLVER
821  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
822  INTEGER(INTG), INTENT(OUT) :: ERR
823  TYPE(varying_string), INTENT(OUT) :: ERROR
824  !Local Variables
825 
826  enters("CELLML_EQUATIONS_CREATE_START",err,error,*999)
827 
828  IF(ASSOCIATED(solver)) THEN
829  IF(ASSOCIATED(cellml_equations)) THEN
830  CALL flagerror("CellML equations is already associated.",err,error,*999)
831  ELSE
832  NULLIFY(cellml_equations)
833  CALL cellml_equations_initialise(solver,err,error,*999)
834  cellml_equations=>solver%CELLML_EQUATIONS
835  ENDIF
836  ELSE
837  CALL flagerror("Solver is not associated.",err,error,*999)
838  ENDIF
839 
840  exits("CELLML_EQUATIONS_CREATE_START")
841  RETURN
842 999 errorsexits("CELLML_EQUATIONS_CREATE_START",err,error)
843  RETURN 1
844 
845  END SUBROUTINE cellml_equations_create_start
846 
847  !
848  !================================================================================================================================
849  !
850 
852  SUBROUTINE cellml_equations_destroy(CELLML_EQUATIONS,ERR,ERROR,*)
854  !Argument variables
855  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
856  INTEGER(INTG), INTENT(OUT) :: ERR
857  TYPE(varying_string), INTENT(OUT) :: ERROR
858  !Local Variables
859 
860  enters("CELLML_EQUATIONS_DESTROY",err,error,*999)
861 
862  IF(ASSOCIATED(cellml_equations)) THEN
863  CALL cellml_equations_finalise(cellml_equations,err,error,*999)
864  ELSE
865  CALL flagerror("CellML equations is not associated.",err,error,*999)
866  ENDIF
867 
868  exits("CELLML_EQUATIONS_DESTROY")
869  RETURN
870 999 errorsexits("CELLML_EQUATIONS_DESTROY",err,error)
871  RETURN 1
872 
873  END SUBROUTINE cellml_equations_destroy
874 
875  !
876  !================================================================================================================================
877  !
878 
880  SUBROUTINE cellml_equations_finalise(CELLML_EQUATIONS,ERR,ERROR,*)
882  !Argument variables
883  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
884  INTEGER(INTG), INTENT(OUT) :: ERR
885  TYPE(varying_string), INTENT(OUT) :: ERROR
886  !Local Variables
887 
888  enters("CELLML_EQUATIONS_FINALISE",err,error,*999)
889 
890  IF(ASSOCIATED(cellml_equations)) THEN
891  IF(ALLOCATED(cellml_equations%CELLML_ENVIRONMENTS)) DEALLOCATE(cellml_equations%CELLML_ENVIRONMENTS)
892  DEALLOCATE(cellml_equations)
893  ENDIF
894 
895  exits("CELLML_EQUATIONS_FINALISE")
896  RETURN
897 999 errorsexits("CELLML_EQUATIONS_FINALISE",err,error)
898  RETURN 1
899 
900  END SUBROUTINE cellml_equations_finalise
901 
902  !
903  !================================================================================================================================
904  !
905 
907  SUBROUTINE cellml_equations_initialise(SOLVER,ERR,ERROR,*)
909  !Argument variables
910  TYPE(solver_type), POINTER :: SOLVER
911  INTEGER(INTG), INTENT(OUT) :: ERR
912  TYPE(varying_string), INTENT(OUT) :: ERROR
913  !Local Variables
914  INTEGER(INTG) :: DUMMY_ERR
915  TYPE(varying_string) :: DUMMY_ERROR
916 
917  enters("CELLML_EQUATIONS_INITIALISE",err,error,*998)
918 
919  IF(ASSOCIATED(solver)) THEN
920  IF(ASSOCIATED(solver%CELLML_EQUATIONS)) THEN
921  CALL flagerror("CellML equations is already associated for this solver.",err,error,*998)
922  ELSE
923  ALLOCATE(solver%CELLML_EQUATIONS,stat=err)
924  IF(err/=0) CALL flagerror("Could not allocate CellML equations.",err,error,*999)
925  solver%CELLML_EQUATIONS%SOLVER=>solver
926  solver%CELLML_EQUATIONS%CELLML_EQUATIONS_FINISHED=.false.
927  solver%CELLML_EQUATIONS%NUMBER_OF_CELLML_ENVIRONMENTS=0
928  ENDIF
929  ELSE
930  CALL flagerror("Solver is not associated.",err,error,*998)
931  ENDIF
932 
933  exits("CELLML_EQUATIONS_INITIALISE")
934  RETURN
935 999 CALL cellml_equations_finalise(solver%CELLML_EQUATIONS,dummy_err,dummy_error,*998)
936 998 errorsexits("CELLML_EQUATIONS_INITIALISE",err,error)
937  RETURN 1
938 
939  END SUBROUTINE cellml_equations_initialise
940 
941  !
942  !================================================================================================================================
943  !
944 
946  SUBROUTINE solver_cellml_equations_get(SOLVER,CELLML_EQUATIONS,ERR,ERROR,*)
948  !Argument variables
949  TYPE(solver_type), POINTER :: SOLVER
950  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
951  INTEGER(INTG), INTENT(OUT) :: ERR
952  TYPE(varying_string), INTENT(OUT) :: ERROR
953  !Local Variables
954 
955  enters("SOLVER_CELLML_EQUATIONS_GET",err,error,*998)
956 
957  IF(ASSOCIATED(solver)) THEN
958  IF(ASSOCIATED(cellml_equations)) THEN
959  CALL flagerror("CellML equations is already associated.",err,error,*998)
960  ELSE
961  cellml_equations=>solver%CELLML_EQUATIONS
962  IF(.NOT.ASSOCIATED(cellml_equations)) CALL flagerror("CellML equations is not associated.",err,error,*999)
963  ENDIF
964  !ELSE
965  !CALL FlagError("Solver has not been finished.",ERR,ERROR,*998)
966  !ENDIF
967  ELSE
968  CALL flagerror("Solver is not associated.",err,error,*998)
969  ENDIF
970 
971  exits("SOLVER_CELLML_EQUATIONS_GET")
972  RETURN
973 999 NULLIFY(cellml_equations)
974 998 errorsexits("SOLVER_CELLML_EQUATIONS_GET",err,error)
975  RETURN 1
976 
977  END SUBROUTINE solver_cellml_equations_get
978 
979  !
980  !================================================================================================================================
981  !
982 
984  SUBROUTINE solver_cellml_evaluator_create_finish(CELLML_EVALUATOR_SOLVER,ERR,ERROR,*)
986  !Argument variables
987  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
988  INTEGER(INTG), INTENT(OUT) :: ERR
989  TYPE(varying_string), INTENT(OUT) :: ERROR
990  !Local Variables
991 
992  enters("SOLVER_CELLML_EVALUATOR_CREATE_FINISH",err,error,*999)
993 
994  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
995  CALL flagerror("Not implemented.",err,error,*999)
996  ELSE
997  CALL flagerror("CellML evaluastor solver is not associated.",err,error,*999)
998  ENDIF
999 
1000  exits("SOLVER_CELLML_EVALUATOR_CREATE_FINISH")
1001  RETURN
1002 999 errorsexits("SOLVER_CELLML_EVALUATOR_CREATE_FINISH",err,error)
1003  RETURN 1
1004 
1006 
1007  !
1008  !================================================================================================================================
1009  !
1010 
1012  SUBROUTINE solver_cellml_evaluator_finalise(CELLML_EVALUATOR_SOLVER,ERR,ERROR,*)
1014  !Argument variables
1015  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
1016  INTEGER(INTG), INTENT(OUT) :: ERR
1017  TYPE(varying_string), INTENT(OUT) :: ERROR
1018  !Local Variables
1019 
1020  enters("SOLVER_CELLML_EVALUATOR_FINALISE",err,error,*999)
1021 
1022  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
1023  DEALLOCATE(cellml_evaluator_solver)
1024  ENDIF
1025 
1026  exits("SOLVER_CELLML_EVALUATOR_FINALISE")
1027  RETURN
1028 999 errorsexits("SOLVER_CELLML_EVALUATOR_FINALISE",err,error)
1029  RETURN 1
1030 
1031  END SUBROUTINE solver_cellml_evaluator_finalise
1032 
1033  !
1034  !================================================================================================================================
1035  !
1036 
1038  SUBROUTINE solver_cellml_evaluator_initialise(SOLVER,ERR,ERROR,*)
1040  !Argument variables
1041  TYPE(solver_type), POINTER :: SOLVER
1042  INTEGER(INTG), INTENT(OUT) :: ERR
1043  TYPE(varying_string), INTENT(OUT) :: ERROR
1044  !Local Variables
1045  INTEGER(INTG) :: DUMMY_ERR
1046  TYPE(varying_string) :: DUMMY_ERROR
1047 
1048  enters("SOLVER_CELLML_EVALUATOR_INITIALISE",err,error,*998)
1049 
1050  IF(ASSOCIATED(solver)) THEN
1051  IF(ASSOCIATED(solver%CELLML_EVALUATOR_SOLVER)) THEN
1052  CALL flagerror("CellML evaluator solver is already associated for this solver.",err,error,*998)
1053  ELSE
1054  ALLOCATE(solver%CELLML_EVALUATOR_SOLVER,stat=err)
1055  IF(err/=0) CALL flagerror("Could not allocate solver CellML evaluator solver.",err,error,*999)
1056  solver%CELLML_EVALUATOR_SOLVER%SOLVER=>solver
1057  solver%CELLML_EVALUATOR_SOLVER%SOLVER_LIBRARY=solver_cmiss_library
1058  solver%CELLML_EVALUATOR_SOLVER%CURRENT_TIME=0.0_dp
1059  ENDIF
1060  ELSE
1061  CALL flagerror("Solver is not associated.",err,error,*998)
1062  ENDIF
1063 
1064  exits("SOLVER_CELLML_EVALUATOR_INITIALISE")
1065  RETURN
1066 999 CALL solver_cellml_evaluator_finalise(solver%CELLML_EVALUATOR_SOLVER,dummy_err,dummy_error,*998)
1067 998 errorsexits("SOLVER_CELLML_EVALUATOR_INITIALISE",err,error)
1068  RETURN 1
1069 
1070  END SUBROUTINE solver_cellml_evaluator_initialise
1071 
1072  !
1073  !================================================================================================================================
1074  !
1075 
1077  SUBROUTINE solver_cellml_evaluator_library_type_get(CELLML_EVALUATOR_SOLVER,SOLVER_LIBRARY_TYPE,ERR,ERROR,*)
1079  !Argument variables
1080  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
1081  INTEGER(INTG), INTENT(OUT) :: SOLVER_LIBRARY_TYPE
1082  INTEGER(INTG), INTENT(OUT) :: ERR
1083  TYPE(varying_string), INTENT(OUT) :: ERROR
1084  !Local Variables
1085 
1086  enters("SOLVER_CELLML_EVALUATOR_LIBRARY_TYPE_GET",err,error,*999)
1087 
1088  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
1089  solver_library_type=cellml_evaluator_solver%SOLVER_LIBRARY
1090  ELSE
1091  CALL flagerror("CellML evaluator solver is not associated.",err,error,*999)
1092  ENDIF
1093 
1094  exits("SOLVER_CELLML_EVALUATOR_LIBRARY_TYPE_GET")
1095  RETURN
1096 999 errorsexits("SOLVER_CELLML_EVALUATOR_LIBRARY_TYPE_GET",err,error)
1097  RETURN 1
1098 
1100 
1101  !
1102  !================================================================================================================================
1103  !
1104 
1106  SUBROUTINE solver_cellml_evaluator_library_type_set(CELLML_EVALUATOR_SOLVER,SOLVER_LIBRARY_TYPE,ERR,ERROR,*)
1108  !Argument variables
1109  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
1110  INTEGER(INTG), INTENT(IN) :: SOLVER_LIBRARY_TYPE
1111  INTEGER(INTG), INTENT(OUT) :: ERR
1112  TYPE(varying_string), INTENT(OUT) :: ERROR
1113  !Local Variables
1114  TYPE(varying_string) :: LOCAL_ERROR
1115 
1116  enters("SOLVER_CELLML_EVALUATOR_LIBRARY_TYPE_SET",err,error,*999)
1117 
1118  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
1119  SELECT CASE(solver_library_type)
1120  CASE(solver_cmiss_library)
1121  CALL flagerror("Not implemented.",err,error,*999)
1122  CASE DEFAULT
1123  local_error="The specified solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
1124  & " is invalid for a CellML evaluator solver."
1125  CALL flagerror(local_error,err,error,*999)
1126  END SELECT
1127  ELSE
1128  CALL flagerror("CellML evaluator solver is not associated.",err,error,*999)
1129  ENDIF
1130 
1131  exits("SOLVER_CELLML_EVALUATOR_LIBRARY_TYPE_SET")
1132  RETURN
1133 999 errorsexits("SOLVER_CELLML_EVALUATOR_LIBRARY_TYPE_SET",err,error)
1134  RETURN 1
1135 
1137 
1138  !
1139  !================================================================================================================================
1140  !
1141 
1143  SUBROUTINE solver_cellml_evaluator_time_get(CELLML_EVALUATOR_SOLVER,TIME,ERR,ERROR,*)
1145  !Argument variables
1146  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
1147  REAL(DP), INTENT(OUT) :: TIME
1148  INTEGER(INTG), INTENT(OUT) :: ERR
1149  TYPE(varying_string), INTENT(OUT) :: ERROR
1150  !Local Variables
1151 
1152  enters("SOLVER_CELLML_EVALUATOR_TIME_GET",err,error,*999)
1153 
1154  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
1155  time=cellml_evaluator_solver%CURRENT_TIME
1156  ELSE
1157  CALL flagerror("CellML evaluator solver is not associated.",err,error,*999)
1158  ENDIF
1159 
1160  exits("SOLVER_CELLML_EVALUATOR_TIME_GET")
1161  RETURN
1162 999 errorsexits("SOLVER_CELLML_EVALUATOR_TIME_GET",err,error)
1163  RETURN 1
1164 
1165  END SUBROUTINE solver_cellml_evaluator_time_get
1166 
1167  !
1168  !================================================================================================================================
1169  !
1170 
1172  SUBROUTINE solver_cellml_evaluator_time_set(CELLML_EVALUATOR_SOLVER,TIME,ERR,ERROR,*)
1174  !Argument variables
1175  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
1176  REAL(DP), INTENT(IN) :: TIME
1177  INTEGER(INTG), INTENT(OUT) :: ERR
1178  TYPE(varying_string), INTENT(OUT) :: ERROR
1179 
1180  enters("SOLVER_CELLML_EVALUATOR_TIME_SET",err,error,*999)
1181 
1182  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
1183  cellml_evaluator_solver%CURRENT_TIME=time
1184  ELSE
1185  CALL flagerror("CellML evaluator solver is not associated.",err,error,*999)
1186  ENDIF
1187 
1188  exits("SOLVER_CELLML_EVALUATOR_TIME_SET")
1189  RETURN
1190 999 errorsexits("SOLVER_CELLML_EVALUATOR_TIME_SET",err,error)
1191  RETURN 1
1192 
1193  END SUBROUTINE solver_cellml_evaluator_time_set
1194 
1195  !
1196  !================================================================================================================================
1197  !
1198 
1200  SUBROUTINE solver_cellml_evaluator_solve(CELLML_EVALUATOR_SOLVER,ERR,ERROR,*)
1202  !Argument variables
1203  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
1204  INTEGER(INTG), INTENT(OUT) :: ERR
1205  TYPE(varying_string), INTENT(OUT) :: ERROR
1206  !Local Variables
1207  INTEGER(INTG) :: cellml_idx
1208  INTEGER(INTG), POINTER :: MODELS_DATA(:)
1209  REAL(DP), POINTER :: INTERMEDIATE_DATA(:),PARAMETERS_DATA(:),STATE_DATA(:)
1210  TYPE(cellml_type), POINTER :: CELLML_ENVIRONMENT
1211  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
1212  TYPE(cellml_models_field_type), POINTER :: CELLML_MODELS_FIELD
1213  TYPE(field_variable_type), POINTER :: MODELS_VARIABLE
1214  TYPE(field_type), POINTER :: MODELS_FIELD,STATE_FIELD,PARAMETERS_FIELD,INTERMEDIATE_FIELD
1215  TYPE(solver_type), POINTER :: SOLVER
1216  TYPE(varying_string) :: LOCAL_ERROR
1217 
1218  enters("SOLVER_CELLML_EVALUATOR_SOLVE",err,error,*999)
1219 
1220  NULLIFY(models_data)
1221  NULLIFY(intermediate_data)
1222  NULLIFY(parameters_data)
1223  NULLIFY(state_data)
1224 
1225  NULLIFY(models_variable)
1226  NULLIFY(state_field)
1227  NULLIFY(parameters_field)
1228  NULLIFY(intermediate_field)
1229 
1230  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
1231  solver=>cellml_evaluator_solver%SOLVER
1232  IF(ASSOCIATED(solver)) THEN
1233  cellml_equations=>solver%CELLML_EQUATIONS
1234  IF(ASSOCIATED(cellml_equations)) THEN
1235  DO cellml_idx=1,cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS
1236  cellml_environment=>cellml_equations%CELLML_ENVIRONMENTS(cellml_idx)%PTR
1237  IF(ASSOCIATED(cellml_environment)) THEN
1238  cellml_models_field=>cellml_environment%MODELS_FIELD
1239  IF(ASSOCIATED(cellml_models_field)) THEN
1240  models_field=>cellml_models_field%MODELS_FIELD
1241  IF(ASSOCIATED(models_field)) THEN
1242 
1243 !!TODO: Maybe move this getting of fields earlier up the DAE solver chain? For now keep here.
1244 
1245  !Make sure CellML fields have been updated to the current value of any mapped fields
1246  CALL cellml_field_to_cellml_update(cellml_environment,err,error,*999)
1247 
1248  CALL field_variable_get(models_field,field_u_variable_type,models_variable,err,error,*999)
1249  CALL field_parameter_set_data_get(models_field,field_u_variable_type,field_values_set_type, &
1250  & models_data,err,error,*999)
1251 
1252  !Get the state information if this environment has any.
1253  IF(ASSOCIATED(cellml_environment%STATE_FIELD)) THEN
1254  state_field=>cellml_environment%STATE_FIELD%STATE_FIELD
1255  IF(ASSOCIATED(state_field)) THEN
1256  CALL field_parameter_set_data_get(state_field,field_u_variable_type,field_values_set_type, &
1257  & state_data,err,error,*999)
1258  ENDIF
1259  ENDIF
1260 
1261  !Get the parameters information if this environment has any.
1262  IF(ASSOCIATED(cellml_environment%PARAMETERS_FIELD)) THEN
1263  parameters_field=>cellml_environment%PARAMETERS_FIELD%PARAMETERS_FIELD
1264  IF(ASSOCIATED(parameters_field)) THEN
1265  CALL field_parameter_set_data_get(parameters_field,field_u_variable_type,field_values_set_type, &
1266  & parameters_data,err,error,*999)
1267  ENDIF
1268  ENDIF
1269 
1270  !Get the intermediate information if this environment has any.
1271  IF(ASSOCIATED(cellml_environment%INTERMEDIATE_FIELD)) THEN
1272  intermediate_field=>cellml_environment%INTERMEDIATE_FIELD%INTERMEDIATE_FIELD
1273  IF(ASSOCIATED(intermediate_field)) THEN
1274  CALL field_parameter_set_data_get(intermediate_field,field_u_variable_type,field_values_set_type, &
1275  & intermediate_data,err,error,*999)
1276  ENDIF
1277  ENDIF
1278 
1279  !Solve these CellML equations
1280  SELECT CASE(cellml_evaluator_solver%SOLVER_LIBRARY)
1281  CASE(solver_cmiss_library)
1282  CALL solver_cellml_evaluate(cellml_evaluator_solver,cellml_environment,models_variable%TOTAL_NUMBER_OF_DOFS, &
1283  & cellml_environment%MODELS_FIELD%ONLY_ONE_MODEL_INDEX,models_data,cellml_environment% &
1284  & maximum_number_of_state,state_data,cellml_environment%MAXIMUM_NUMBER_OF_PARAMETERS, &
1285  & parameters_data,cellml_environment%MAXIMUM_NUMBER_OF_INTERMEDIATE,intermediate_data,err,error,*999)
1286  CASE DEFAULT
1287  CALL flagerror("Solver library not implemented.",err,error,*999)
1288  END SELECT
1289 
1290  !Restore field data
1291  CALL field_parameter_set_data_restore(models_field,field_u_variable_type,field_values_set_type, &
1292  & models_data,err,error,*999)
1293  IF(ASSOCIATED(state_field)) CALL field_parameter_set_data_restore(state_field,field_u_variable_type, &
1294  & field_values_set_type,state_data,err,error,*999)
1295  IF(ASSOCIATED(parameters_field)) CALL field_parameter_set_data_restore(parameters_field, &
1296  & field_u_variable_type,field_values_set_type,parameters_data,err,error,*999)
1297  IF(ASSOCIATED(intermediate_field)) CALL field_parameter_set_data_restore(intermediate_field, &
1298  & field_u_variable_type,field_values_set_type,intermediate_data,err,error,*999)
1299 
1300  !Make sure fields have been updated to the current value of any mapped CellML fields
1301  CALL cellml_cellml_to_field_update(cellml_environment,err,error,*999)
1302 
1303  ELSE
1304  local_error="The CellML models field is not associated for CellML index "// &
1305  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
1306  CALL flagerror(local_error,err,error,*999)
1307  ENDIF
1308  ELSE
1309  local_error="The CellML models field is not associated for CellML index "// &
1310  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
1311  CALL flagerror(local_error,err,error,*999)
1312  ENDIF
1313  ELSE
1314  local_error="The CellML enviroment is not associated for for CellML index "// &
1315  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
1316  CALL flagerror(local_error,err,error,*999)
1317  ENDIF
1318  ENDDO !cellml_idx
1319  ELSE
1320  CALL flagerror("Solver solver equations is not associated.",err,error,*999)
1321  ENDIF
1322  ELSE
1323  CALL flagerror("Solver is not associated.",err,error,*999)
1324  ENDIF
1325  ELSE
1326  CALL flagerror("CellML evaluator solver is not associated.",err,error,*999)
1327  ENDIF
1328 
1329  exits("SOLVER_CELLML_EVALUATOR_SOLVE")
1330  RETURN
1331 999 errorsexits("SOLVER_CELLML_EVALUATOR_SOLVE",err,error)
1332  RETURN 1
1333 
1334  END SUBROUTINE solver_cellml_evaluator_solve
1335 
1336  !
1337  !================================================================================================================================
1338  !
1339 
1341  SUBROUTINE solver_cellml_evaluate(CELLML_EVALUATOR_SOLVER,CELLML,N, ONLY_ONE_MODEL_INDEX,MODELS_DATA,MAX_NUMBER_STATES, &
1342  & state_data,max_number_parameters,parameters_data,max_number_intermediates,intermediate_data,err,error,*)
1344  !Argument variables
1345  TYPE(cellml_evaluator_solver_type), POINTER :: CELLML_EVALUATOR_SOLVER
1346  TYPE(cellml_type), POINTER :: CELLML
1347  INTEGER(INTG), INTENT(IN) :: N
1348  INTEGER(INTG), INTENT(IN) :: ONLY_ONE_MODEL_INDEX
1349  INTEGER(INTG), POINTER :: MODELS_DATA(:)
1350  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_STATES
1351  REAL(DP), POINTER :: STATE_DATA(:)
1352  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_PARAMETERS
1353  REAL(DP), POINTER :: PARAMETERS_DATA(:)
1354  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_INTERMEDIATES
1355  REAL(DP), POINTER :: INTERMEDIATE_DATA(:)
1356  INTEGER(INTG), INTENT(OUT) :: ERR
1357  TYPE(varying_string), INTENT(OUT) :: ERROR
1358  !Local Variables
1359  INTEGER(INTG) :: dof_idx,DOF_ORDER_TYPE,INTERMEDIATE_END_DOF,intermediate_idx,INTERMEDIATE_START_DOF,model_idx, &
1360  & NUMBER_INTERMEDIATES,NUMBER_PARAMETERS,NUMBER_STATES,PARAMETER_END_DOF,parameter_idx,PARAMETER_START_DOF, &
1361  & STATE_END_DOF,state_idx,STATE_START_DOF
1362  REAL(DP) :: INTERMEDIATES(max(1,max_number_intermediates)),PARAMETERS(max(1,max_number_parameters)), &
1363  & RATES(MAX(1,MAX_NUMBER_STATES)),STATES(MAX(1,MAX_NUMBER_STATES))
1364  TYPE(cellml_model_type), POINTER :: MODEL
1365  TYPE(varying_string) :: LOCAL_ERROR
1366 
1367  enters("SOLVER_CELLML_EVALUATE",err,error,*999)
1368 
1369  IF(ASSOCIATED(cellml_evaluator_solver)) THEN
1370  IF(ASSOCIATED(cellml)) THEN
1371  IF(ASSOCIATED(cellml%MODELS_FIELD)) THEN
1372  CALL field_dof_order_type_get(cellml%MODELS_FIELD%MODELS_FIELD,field_u_variable_type,dof_order_type,err,error,*999)
1373  IF(dof_order_type==field_separated_component_dof_order) THEN
1374  !Dof components are separated. Will need to copy data to temporary arrays.
1375  IF(only_one_model_index==cellml_models_field_not_constant) THEN
1376  !Mulitple models
1377  DO dof_idx=1,n
1378  model_idx=models_data(dof_idx)
1379  IF(model_idx.GT.0) THEN
1380  model=>cellml%MODELS(model_idx)%PTR
1381  IF(ASSOCIATED(model)) THEN
1382  number_states=model%NUMBER_OF_STATE
1383  number_intermediates=model%NUMBER_OF_INTERMEDIATE
1384  number_parameters=model%NUMBER_OF_PARAMETERS
1385 
1386  !Copy CellML data to temporary arrays
1387  DO state_idx=1,number_states
1388  states(state_idx)=state_data((dof_idx-1)*n+state_idx)
1389  ENDDO !state_idx
1390  DO parameter_idx=1,number_parameters
1391  parameters(parameter_idx)=parameters_data((dof_idx-1)*n+parameter_idx)
1392  ENDDO !parameter_idx
1393 
1394 #ifdef WITH_CELLML
1395  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp,states,rates,intermediates, &
1396  & parameters)
1397 #else
1398  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
1399 #endif
1400 
1401  !Copy temporary data back to CellML arrays
1402  DO intermediate_idx=1,number_intermediates
1403  intermediate_data((dof_idx-1)*n+intermediate_idx)=intermediates(intermediate_idx)
1404  ENDDO !intermediate_idx
1405  DO state_idx=1,number_states
1406  state_data((dof_idx-1)*n+state_idx)=states(state_idx)
1407  ENDDO !state_idx
1408 
1409  ELSE
1410  local_error="CellML environment model is not associated for model index "// &
1411  & trim(numbertovstring(only_one_model_index,"*",err,error))//" belonging to dof index "// &
1412  & trim(numbertovstring(dof_idx,"*",err,error))//"."
1413  CALL flagerror(local_error,err,error,*999)
1414  ENDIF
1415  ENDIF !model_idx
1416  ENDDO !dof_idx
1417  ELSE
1418  !One one model is used.
1419  model=>cellml%MODELS(only_one_model_index)%PTR
1420  IF(ASSOCIATED(model)) THEN
1421  number_states=model%NUMBER_OF_STATE
1422  number_intermediates=model%NUMBER_OF_INTERMEDIATE
1423  number_parameters=model%NUMBER_OF_PARAMETERS
1424  DO dof_idx=1,n
1425  model_idx=models_data(dof_idx)
1426  IF(model_idx.GT.0) THEN
1427 
1428  !Copy CellML data to temporary arrays
1429  DO state_idx=1,number_states
1430  states(state_idx)=state_data((dof_idx-1)*n+state_idx)
1431  ENDDO !state_idx
1432  DO parameter_idx=1,number_parameters
1433  parameters(parameter_idx)=parameters_data((dof_idx-1)*n+parameter_idx)
1434  ENDDO !parameter_idx
1435 
1436 #ifdef WITH_CELLML
1437  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp,states,rates,intermediates, &
1438  & parameters)
1439 #else
1440  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
1441 #endif
1442 
1443  !Copy temporary data back to CellML arrays
1444  DO intermediate_idx=1,number_intermediates
1445  intermediate_data((dof_idx-1)*n+intermediate_idx)=intermediates(intermediate_idx)
1446  ENDDO !intermediate_idx
1447  DO state_idx=1,number_states
1448  state_data((dof_idx-1)*n+state_idx)=states(state_idx)
1449  ENDDO !state_idx
1450  ENDIF !model_idx
1451  ENDDO !dof_idx
1452  ELSE
1453  local_error="CellML environment model is not associated for model index "// &
1454  & trim(numbertovstring(only_one_model_index,"*",err,error))//"."
1455  CALL flagerror(local_error,err,error,*999)
1456  ENDIF
1457  ENDIF
1458  ELSE
1459  !Dof components are continguous. Can pass data directly.
1460  IF(only_one_model_index==cellml_models_field_not_constant) THEN
1461  !Mulitple models
1462 
1463 #ifdef WITH_CELLML
1464 
1465  DO dof_idx=1,n
1466  model_idx=models_data(dof_idx)
1467  IF(model_idx.GT.0) THEN
1468  model=>cellml%MODELS(model_idx)%PTR
1469  IF(ASSOCIATED(model)) THEN
1470  number_states=model%NUMBER_OF_STATE
1471  number_intermediates=model%NUMBER_OF_INTERMEDIATE
1472  number_parameters=model%NUMBER_OF_PARAMETERS
1473  !Call RHS. Note some models might not have state, rates, intermediate or parameter data so call accordingly
1474  !to avoid indexing in to null pointers
1475  IF(number_states>0) THEN
1476  IF(number_intermediates>0) THEN
1477  IF(number_parameters>0) THEN
1478  !We have state, intermediate and parameters in the model
1479  state_start_dof=(dof_idx-1)*max_number_states+1
1480  state_end_dof=state_start_dof+number_states-1
1481  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1482  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1483  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
1484  parameter_end_dof=parameter_start_dof+number_parameters-1
1485 
1486  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp, &
1487  & state_data(state_start_dof:state_end_dof), &
1488  & rates,intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters_data( &
1489  & parameter_start_dof:parameter_end_dof))
1490 
1491  ELSE
1492  !We do not have parameters in the model
1493  state_start_dof=(dof_idx-1)*max_number_states+1
1494  state_end_dof=state_start_dof+number_states-1
1495  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1496  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1497 
1498  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp, &
1499  & state_data(state_start_dof:state_end_dof), &
1500  & rates,intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters)
1501 
1502  ENDIF
1503  ELSE
1504  IF(number_parameters>0) THEN
1505  !We do not have intermediates in the model
1506  state_start_dof=(dof_idx-1)*max_number_states+1
1507  state_end_dof=state_start_dof+number_states-1
1508  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
1509  parameter_end_dof=parameter_start_dof+number_parameters-1
1510 
1511  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp, &
1512  & state_data(state_start_dof:state_end_dof), &
1513  & rates,intermediates,parameters_data(parameter_start_dof:parameter_end_dof))
1514 
1515  ELSE
1516  !We do not have intermediates or parameters in the model
1517  state_start_dof=(dof_idx-1)*max_number_states+1
1518  state_end_dof=state_start_dof+number_states-1
1519 
1520  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp, &
1521  & state_data(state_start_dof:state_end_dof), &
1522  & rates,intermediates,parameters)
1523 
1524  ENDIF
1525  ENDIF
1526  ELSE
1527  IF(number_intermediates>0) THEN
1528  IF(number_parameters>0) THEN
1529  !We do not have any states in the model
1530  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1531  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1532  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
1533  parameter_end_dof=parameter_start_dof+number_parameters-1
1534 
1535  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp,states,rates, &
1536  & intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters_data( &
1537  & parameter_start_dof:parameter_end_dof))
1538  ELSE
1539  !We do not have any states or parameters in the model
1540  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1541  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1542 
1543  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp,states,rates, &
1544  & intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters)
1545 
1546  ENDIF
1547  ELSE
1548  CALL flagerror("Invalid CellML model - there are no states or intermediates.",err,error,*999)
1549  ENDIF
1550  ENDIF
1551 
1552 
1553  ELSE
1554  local_error="CellML environment model is not associated for model index "// &
1555  & trim(numbertovstring(only_one_model_index,"*",err,error))//" belonging to dof index "// &
1556  & trim(numbertovstring(dof_idx,"*",err,error))//"."
1557  CALL flagerror(local_error,err,error,*999)
1558  ENDIF
1559  ENDIF !model_idx
1560  ENDDO !dof_idx
1561 #else
1562  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
1563 #endif
1564 
1565  ELSE
1566  !One model is used.
1567  model=>cellml%MODELS(only_one_model_index)%PTR
1568  IF(ASSOCIATED(model)) THEN
1569  number_states=model%NUMBER_OF_STATE
1570  number_intermediates=model%NUMBER_OF_INTERMEDIATE
1571  number_parameters=model%NUMBER_OF_PARAMETERS
1572 #ifdef WITH_CELLML
1573  !Call RHS. Note some models might not have state, rates, intermediate or parameter data so call accordingly
1574  !to avoid referencing null pointers
1575  IF(number_states>0) THEN
1576  IF(number_intermediates>0) THEN
1577  IF(number_parameters>0) THEN
1578  !We have states, intermediate and parameters for the model
1579  DO dof_idx=1,n
1580  model_idx=models_data(dof_idx)
1581  IF(model_idx.GT.0) THEN
1582  state_start_dof=(dof_idx-1)*max_number_states+1
1583  state_end_dof=state_start_dof+number_states-1
1584  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1585  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1586  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
1587  parameter_end_dof=parameter_start_dof+number_parameters-1
1588 
1589  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp, &
1590  & state_data(state_start_dof:state_end_dof), &
1591  & rates,intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters_data( &
1592  & parameter_start_dof:parameter_end_dof))
1593  ENDIF !model_idx
1594  ENDDO !dof_idx
1595  ELSE
1596  !We do not have parameters in the model
1597  DO dof_idx=1,n
1598  model_idx=models_data(dof_idx)
1599  IF(model_idx.GT.0) THEN
1600  state_start_dof=(dof_idx-1)*max_number_states+1
1601  state_end_dof=state_start_dof+number_states-1
1602  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1603  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1604 
1605  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp, &
1606  & state_data(state_start_dof:state_end_dof), &
1607  & rates,intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters)
1608  ENDIF !model_idx
1609  ENDDO !dof_idx
1610 
1611  ENDIF
1612  ELSE
1613  IF(number_parameters>0) THEN
1614  !We do not have any intermediates in the model
1615  DO dof_idx=1,n
1616  model_idx=models_data(dof_idx)
1617  IF(model_idx.GT.0) THEN
1618 
1619  state_start_dof=(dof_idx-1)*max_number_states+1
1620  state_end_dof=state_start_dof+number_states-1
1621  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
1622  parameter_end_dof=parameter_start_dof+number_parameters-1
1623 
1624  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp, &
1625  & state_data(state_start_dof:state_end_dof), &
1626  & rates,intermediates,parameters_data(parameter_start_dof:parameter_end_dof))
1627  ENDIF !model_idx
1628  ENDDO !dof_idx
1629  ELSE
1630  !We do not have any intermediates or parameters in the model
1631  DO dof_idx=1,n
1632  model_idx=models_data(dof_idx)
1633  IF(model_idx.GT.0) THEN
1634 
1635  state_start_dof=(dof_idx-1)*max_number_states+1
1636  state_end_dof=state_start_dof+number_states-1
1637 
1638  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp,&
1639  & state_data(state_start_dof:state_end_dof), &
1640  & rates,intermediates,parameters)
1641  ENDIF !model_idx
1642  ENDDO !dof_idx
1643  ENDIF
1644  ENDIF
1645  ELSE
1646  IF(number_intermediates>0) THEN
1647  IF(number_parameters>0) THEN
1648  !We do not have any states in the model
1649  DO dof_idx=1,n
1650  model_idx=models_data(dof_idx)
1651  IF(model_idx.GT.0) THEN
1652 
1653  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1654  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1655  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
1656  parameter_end_dof=parameter_start_dof+number_parameters-1
1657 
1658  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp,states,rates, &
1659  & intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters_data( &
1660  & parameter_start_dof:parameter_end_dof))
1661  ENDIF !model_idx
1662  ENDDO !dof_idx
1663  ELSE
1664  !We do not have any states or parameters the model
1665  DO dof_idx=1,n
1666  model_idx=models_data(dof_idx)
1667  IF(model_idx.GT.0) THEN
1668 
1669  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
1670  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
1671 
1672  CALL cellml_model_definition_call_rhs_routine(model%PTR,0.0_dp,states,rates, &
1673  & intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters)
1674  ENDIF !model_idx
1675  ENDDO !dof_idx
1676  ENDIF
1677  ELSE
1678  CALL flagerror("Invalid CellML model - there are no states or intermediates.",err,error,*999)
1679  ENDIF
1680  ENDIF
1681 #else
1682  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
1683 #endif
1684  ELSE
1685  local_error="CellML environment model is not associated for model index "// &
1686  & trim(numbertovstring(only_one_model_index,"*",err,error))//"."
1687  CALL flagerror(local_error,err,error,*999)
1688  ENDIF
1689  ENDIF
1690  ENDIF
1691  ELSE
1692  CALL flagerror("CellML environment models field is not associated.",err,error,*999)
1693  ENDIF
1694  ELSE
1695  CALL flagerror("CellML environment is not associated.",err,error,*999)
1696  ENDIF
1697  ELSE
1698  CALL flagerror("CellML evaluator solver is not associated.",err,error,*999)
1699  ENDIF
1700 
1701  exits("SOLVER_CELLML_EVALUATE")
1702  RETURN
1703 999 errorsexits("SOLVER_CELLML_EVALUATE",err,error)
1704  RETURN 1
1705 
1706  END SUBROUTINE solver_cellml_evaluate
1707 
1708  !
1709  !================================================================================================================================
1710  !
1711 
1713  SUBROUTINE solver_create_finish(SOLVER,ERR,ERROR,*)
1715  !Argument variables
1716  TYPE(solver_type), POINTER :: SOLVER
1717  INTEGER(INTG), INTENT(OUT) :: ERR
1718  TYPE(varying_string), INTENT(OUT) :: ERROR
1719  !Local Variables
1720  INTEGER(INTG) :: solver_idx
1721 
1722  enters("SOLVER_CREATE_FINISH",err,error,*999)
1723 
1724  IF(ASSOCIATED(solver)) THEN
1725  IF(solver%SOLVER_FINISHED) THEN
1726  CALL flagerror("Solver has already been finished.",err,error,*999)
1727  ELSE
1728  !Set the finished flag. The final solver finish will be done once the solver equations have been finished.
1729  DO solver_idx=1,solver%NUMBER_OF_LINKED_SOLVERS
1730  solver%LINKED_SOLVERS(solver_idx)%PTR%SOLVER_FINISHED=.true.
1731  ENDDO !solver_idx
1732 
1733  solver%SOLVER_FINISHED=.true.
1734  ENDIF
1735  ELSE
1736  CALL flagerror("Solver is not associated.",err,error,*999)
1737  ENDIF
1738 
1739  exits("SOLVER_CREATE_FINISH")
1740  RETURN
1741 999 errorsexits("SOLVER_CREATE_FINISH",err,error)
1742  RETURN 1
1743 
1744  END SUBROUTINE solver_create_finish
1745 
1746  !
1747  !================================================================================================================================
1748  !
1749 
1751  SUBROUTINE solver_dae_adams_moulton_finalise(ADAMS_MOULTON_SOLVER,ERR,ERROR,*)
1753  !Argument variables
1754  TYPE(adams_moulton_dae_solver_type), POINTER :: ADAMS_MOULTON_SOLVER
1755  INTEGER(INTG), INTENT(OUT) :: ERR
1756  TYPE(varying_string), INTENT(OUT) :: ERROR
1757  !Local Variables
1758 
1759  enters("SOLVER_DAE_ADAMS_MOULTON_FINALISE",err,error,*999)
1760 
1761  IF(ASSOCIATED(adams_moulton_solver)) THEN
1762  DEALLOCATE(adams_moulton_solver)
1763  ENDIF
1764 
1765  exits("SOLVER_DAE_ADAMS_MOULTON_FINALISE")
1766  RETURN
1767 999 errorsexits("SOLVER_DAE_ADAMS_MOULTON_FINALISE",err,error)
1768  RETURN 1
1769 
1770  END SUBROUTINE solver_dae_adams_moulton_finalise
1771 
1772  !
1773  !================================================================================================================================
1774  !
1775 
1777  SUBROUTINE solver_dae_adams_moulton_initialise(DAE_SOLVER,ERR,ERROR,*)
1779  !Argument variables
1780  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
1781  INTEGER(INTG), INTENT(OUT) :: ERR
1782  TYPE(varying_string), INTENT(OUT) :: ERROR
1783  !Local Variables
1784  INTEGER(INTG) :: DUMMY_ERR
1785  TYPE(varying_string) :: DUMMY_ERROR
1786 
1787  enters("SOLVER_DAE_ADAMS_MOULTON_INITIALISE",err,error,*998)
1788 
1789  IF(ASSOCIATED(dae_solver)) THEN
1790  IF(ASSOCIATED(dae_solver%ADAMS_MOULTON_SOLVER)) THEN
1791  CALL flagerror("Adams-Moulton solver is already associated for this differential-algebraic equation solver.", &
1792  & err,error,*998)
1793  ELSE
1794  !Allocate the Adams-Moulton solver
1795  ALLOCATE(dae_solver%ADAMS_MOULTON_SOLVER,stat=err)
1796  IF(err/=0) CALL flagerror("Could not allocate Adams-Moulton solver.",err,error,*999)
1797  !Initialise
1798  dae_solver%ADAMS_MOULTON_SOLVER%DAE_SOLVER=>dae_solver
1799  dae_solver%ADAMS_MOULTON_SOLVER%SOLVER_LIBRARY=0
1800  !Defaults
1801  ENDIF
1802  ELSE
1803  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*998)
1804  ENDIF
1805 
1806  exits("SOLVER_DAE_ADAMS_MOULTON_INITIALISE")
1807  RETURN
1808 999 CALL solver_dae_adams_moulton_finalise(dae_solver%ADAMS_MOULTON_SOLVER,dummy_err,dummy_error,*998)
1809 998 errorsexits("SOLVER_DAE_ADAMS_MOULTON_INITIALISE",err,error)
1810  RETURN 1
1811 
1813 
1814  !
1815  !================================================================================================================================
1816  !
1817 
1819  SUBROUTINE solver_dae_adams_moulton_solve(ADAMS_MOULTON_SOLVER,ERR,ERROR,*)
1821  !Argument variables
1822  TYPE(adams_moulton_dae_solver_type), POINTER :: ADAMS_MOULTON_SOLVER
1823  INTEGER(INTG), INTENT(OUT) :: ERR
1824  TYPE(varying_string), INTENT(OUT) :: ERROR
1825  !Local Variables
1826 
1827  enters("SOLVER_DAE_ADAMS_MOULTON_SOLVE",err,error,*999)
1828 
1829  IF(ASSOCIATED(adams_moulton_solver)) THEN
1830  CALL flagerror("Not implemented.",err,error,*999)
1831  ELSE
1832  CALL flagerror("Adams-Moulton differential-algebraic equation solver is not associated.",err,error,*999)
1833  ENDIF
1834 
1835  exits("SOLVER_DAE_ADAMS_MOULTON_SOLVE")
1836  RETURN
1837 999 errorsexits("SOLVER_DAE_ADAMS_MOULTON_SOLVE",err,error)
1838  RETURN 1
1839 
1840  END SUBROUTINE solver_dae_adams_moulton_solve
1841 
1842  !
1843  !================================================================================================================================
1844  !
1845 
1847  SUBROUTINE solver_dae_create_finish(DAE_SOLVER,ERR,ERROR,*)
1849  !Argument variables
1850  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
1851  INTEGER(INTG), INTENT(OUT) :: ERR
1852  TYPE(varying_string), INTENT(OUT) :: ERROR
1853  !Local Variables
1854 
1855  enters("SOLVER_DAE_CREATE_FINISH",err,error,*999)
1856 
1857  IF(ASSOCIATED(dae_solver)) THEN
1858  CALL flagerror("Not implemented.",err,error,*999)
1859  ELSE
1860  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*999)
1861  ENDIF
1862 
1863  exits("SOLVER_DAE_CREATE_FINISH")
1864  RETURN
1865 999 errorsexits("SOLVER_DAE_CREATE_FINISH",err,error)
1866  RETURN 1
1867 
1868  END SUBROUTINE solver_dae_create_finish
1869 
1870  !
1871  !================================================================================================================================
1872  !
1873 
1875  SUBROUTINE solver_dae_euler_backward_finalise(BACKWARD_EULER_SOLVER,ERR,ERROR,*)
1877  !Argument variables
1878  TYPE(backward_euler_dae_solver_type), POINTER :: BACKWARD_EULER_SOLVER
1879  INTEGER(INTG), INTENT(OUT) :: ERR
1880  TYPE(varying_string), INTENT(OUT) :: ERROR
1881  !Local Variables
1882 
1883  enters("SOLVER_DAE_EULER_BACKWARD_FINALISE",err,error,*999)
1884 
1885  IF(ASSOCIATED(backward_euler_solver)) THEN
1886  DEALLOCATE(backward_euler_solver)
1887  ENDIF
1888 
1889  exits("SOLVER_DAE_EULER_BACKWARD_FINALISE")
1890  RETURN
1891 999 errorsexits("SOLVER_DAE_EULER_BACKWARD_FINALISE",err,error)
1892  RETURN 1
1893 
1894  END SUBROUTINE solver_dae_euler_backward_finalise
1895 
1896  !
1897  !================================================================================================================================
1898  !
1899 
1901  SUBROUTINE solver_dae_euler_backward_initialise(EULER_DAE_SOLVER,ERR,ERROR,*)
1903  !Argument variables
1904  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
1905  INTEGER(INTG), INTENT(OUT) :: ERR
1906  TYPE(varying_string), INTENT(OUT) :: ERROR
1907  !Local Variables
1908  INTEGER(INTG) :: DUMMY_ERR
1909  TYPE(varying_string) :: DUMMY_ERROR
1910 
1911  enters("SOLVER_DAE_EULER_BACKWARD_INITIALISE",err,error,*998)
1912 
1913  IF(ASSOCIATED(euler_dae_solver)) THEN
1914  IF(ASSOCIATED(euler_dae_solver%BACKWARD_EULER_SOLVER)) THEN
1915  CALL flagerror("Backward Euler solver is already associated for this Euler differential-algebraic equation solver.", &
1916  & err,error,*998)
1917  ELSE
1918  !Allocate the backward Euler solver
1919  ALLOCATE(euler_dae_solver%BACKWARD_EULER_SOLVER,stat=err)
1920  IF(err/=0) CALL flagerror("Could not allocate backward Euler solver.",err,error,*999)
1921  !Initialise
1922  euler_dae_solver%BACKWARD_EULER_SOLVER%EULER_DAE_SOLVER=>euler_dae_solver
1923  euler_dae_solver%BACKWARD_EULER_SOLVER%SOLVER_LIBRARY=0
1924  !Defaults
1925  ENDIF
1926  ELSE
1927  CALL flagerror("Euler differential-algebraic equation solver is not associated.",err,error,*998)
1928  ENDIF
1929 
1930  exits("SOLVER_DAE_EULER_BACKWARD_INITIALISE")
1931  RETURN
1932 999 CALL solver_dae_euler_backward_finalise(euler_dae_solver%BACKWARD_EULER_SOLVER,dummy_err,dummy_error,*998)
1933 998 errorsexits("SOLVER_DAE_EULER_BACKWARD_INITIALISE",err,error)
1934  RETURN 1
1935 
1937 
1938  !
1939  !================================================================================================================================
1940  !
1941 
1943  SUBROUTINE solver_dae_euler_backward_solve(BACKWARD_EULER_SOLVER,ERR,ERROR,*)
1945  !Argument variables
1946  TYPE(backward_euler_dae_solver_type), POINTER :: BACKWARD_EULER_SOLVER
1947  INTEGER(INTG), INTENT(OUT) :: ERR
1948  TYPE(varying_string), INTENT(OUT) :: ERROR
1949  !Local Variables
1950 
1951  enters("SOLVER_DAE_EULER_BACKWARD_SOLVE",err,error,*999)
1952 
1953  IF(ASSOCIATED(backward_euler_solver)) THEN
1954  CALL flagerror("Not implemented.",err,error,*999)
1955  ELSE
1956  CALL flagerror("Backward Euler differential-algebraic equation solver is not associated.",err,error,*999)
1957  ENDIF
1958 
1959  exits("SOLVER_DAE_EULER_BACKWARD_SOLVE")
1960  RETURN
1961 999 errorsexits("SOLVER_DAE_EULER_BACKWARD_SOLVE",err,error)
1962  RETURN 1
1963 
1964  END SUBROUTINE solver_dae_euler_backward_solve
1965 
1966  !
1967  !================================================================================================================================
1968  !
1969 
1971  SUBROUTINE solver_dae_euler_finalise(EULER_SOLVER,ERR,ERROR,*)
1973  !Argument variables
1974  TYPE(euler_dae_solver_type), POINTER :: EULER_SOLVER
1975  INTEGER(INTG), INTENT(OUT) :: ERR
1976  TYPE(varying_string), INTENT(OUT) :: ERROR
1977  !Local Variables
1978 
1979  enters("SOLVER_DAE_EULER_FINALISE",err,error,*999)
1980 
1981  IF(ASSOCIATED(euler_solver)) THEN
1982  CALL solver_dae_euler_forward_finalise(euler_solver%FORWARD_EULER_SOLVER,err,error,*999)
1983  CALL solver_dae_euler_backward_finalise(euler_solver%BACKWARD_EULER_SOLVER,err,error,*999)
1984  CALL solver_dae_euler_improved_finalise(euler_solver%IMPROVED_EULER_SOLVER,err,error,*999)
1985  DEALLOCATE(euler_solver)
1986  ENDIF
1987 
1988  exits("SOLVER_DAE_EULER_FINALISE")
1989  RETURN
1990 999 errorsexits("SOLVER_DAE_EULER_FINALISE",err,error)
1991  RETURN 1
1992 
1993  END SUBROUTINE solver_dae_euler_finalise
1994 
1995  !
1996  !================================================================================================================================
1997  !
1998 
2000  SUBROUTINE solver_dae_euler_forward_finalise(FORWARD_EULER_SOLVER,ERR,ERROR,*)
2002  !Argument variables
2003  TYPE(forward_euler_dae_solver_type), POINTER :: FORWARD_EULER_SOLVER
2004  INTEGER(INTG), INTENT(OUT) :: ERR
2005  TYPE(varying_string), INTENT(OUT) :: ERROR
2006  !Local Variables
2007 
2008  enters("SOLVER_DAE_EULER_FORWARD_FINALISE",err,error,*999)
2009 
2010  IF(ASSOCIATED(forward_euler_solver)) THEN
2011  DEALLOCATE(forward_euler_solver)
2012  ENDIF
2013 
2014  exits("SOLVER_DAE_EULER_FORWARD_FINALISE")
2015  RETURN
2016 999 errorsexits("SOLVER_DAE_EULER_FORWARD_FINALISE",err,error)
2017  RETURN 1
2018 
2019  END SUBROUTINE solver_dae_euler_forward_finalise
2020 
2021  !
2022  !================================================================================================================================
2023  !
2024 
2026  SUBROUTINE solver_dae_euler_forward_initialise(EULER_DAE_SOLVER,ERR,ERROR,*)
2028  !Argument variables
2029  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
2030  INTEGER(INTG), INTENT(OUT) :: ERR
2031  TYPE(varying_string), INTENT(OUT) :: ERROR
2032  !Local Variables
2033  INTEGER(INTG) :: DUMMY_ERR
2034  TYPE(varying_string) :: DUMMY_ERROR
2035 
2036  enters("SOLVER_DAE_EULER_FORWARD_INITIALISE",err,error,*998)
2037 
2038  IF(ASSOCIATED(euler_dae_solver)) THEN
2039  IF(ASSOCIATED(euler_dae_solver%FORWARD_EULER_SOLVER)) THEN
2040  CALL flagerror("Forward Euler solver is already associated for this Euler differential-algebraic equation solver.", &
2041  & err,error,*998)
2042  ELSE
2043  !Allocate the forward Euler solver
2044  ALLOCATE(euler_dae_solver%FORWARD_EULER_SOLVER,stat=err)
2045  IF(err/=0) CALL flagerror("Could not allocate forward Euler solver.",err,error,*999)
2046  !Initialise
2047  euler_dae_solver%FORWARD_EULER_SOLVER%EULER_DAE_SOLVER=>euler_dae_solver
2048  euler_dae_solver%FORWARD_EULER_SOLVER%SOLVER_LIBRARY=solver_cmiss_library
2049  !Defaults
2050  ENDIF
2051  ELSE
2052  CALL flagerror("Euler differential-algebraic equation solver is not associated.",err,error,*998)
2053  ENDIF
2054 
2055  exits("SOLVER_DAE_EULER_FORWARD_INITIALISE")
2056  RETURN
2057 999 CALL solver_dae_euler_forward_finalise(euler_dae_solver%FORWARD_EULER_SOLVER,dummy_err,dummy_error,*998)
2058 998 errorsexits("SOLVER_DAE_EULER_FORWARD_INITIALISE",err,error)
2059  RETURN 1
2060 
2062 
2063  !
2064  !================================================================================================================================
2065  !
2066 
2068  SUBROUTINE solver_dae_euler_forward_integrate(FORWARD_EULER_SOLVER,CELLML,N,START_TIME,END_TIME,TIME_INCREMENT, &
2069  & only_one_model_index,models_data,max_number_states,state_data,max_number_parameters,parameters_data, &
2070  & max_number_intermediates,intermediate_data,err,error,*)
2072  !Argument variables
2073  TYPE(forward_euler_dae_solver_type), POINTER :: FORWARD_EULER_SOLVER
2074  TYPE(cellml_type), POINTER :: CELLML
2075  INTEGER(INTG), INTENT(IN) :: N
2076  REAL(DP), INTENT(IN) :: START_TIME
2077  REAL(DP), INTENT(IN) :: END_TIME
2078  REAL(DP), INTENT(INOUT) :: TIME_INCREMENT
2079  INTEGER(INTG), INTENT(IN) :: ONLY_ONE_MODEL_INDEX
2080  INTEGER(INTG), POINTER :: MODELS_DATA(:)
2081  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_STATES
2082  REAL(DP), POINTER :: STATE_DATA(:)
2083  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_PARAMETERS
2084  REAL(DP), POINTER :: PARAMETERS_DATA(:)
2085  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_INTERMEDIATES
2086  REAL(DP), POINTER :: INTERMEDIATE_DATA(:)
2087  INTEGER(INTG), INTENT(OUT) :: ERR
2088  TYPE(varying_string), INTENT(OUT) :: ERROR
2089  !Local Variables
2090  INTEGER(INTG) :: dof_idx,DOF_ORDER_TYPE,INTERMEDIATE_END_DOF,intermediate_idx,INTERMEDIATE_START_DOF,model_idx, &
2091  & NUMBER_INTERMEDIATES,NUMBER_PARAMETERS,NUMBER_STATES,PARAMETER_END_DOF,parameter_idx,PARAMETER_START_DOF, &
2092  & STATE_END_DOF,state_idx,STATE_START_DOF
2093  REAL(DP) :: INTERMEDIATES(max(1,max_number_intermediates)),PARAMETERS(max(1,max_number_parameters)), &
2094  & RATES(MAX(1,MAX_NUMBER_STATES)),STATES(MAX(1,MAX_NUMBER_STATES)),TIME
2095  TYPE(cellml_model_type), POINTER :: MODEL
2096  TYPE(varying_string) :: LOCAL_ERROR
2097 
2098  enters("SOLVER_DAE_EULER_FORWARD_INTEGRATE",err,error,*999)
2099 
2100  IF(ASSOCIATED(forward_euler_solver)) THEN
2101  IF(ASSOCIATED(cellml)) THEN
2102  IF(ASSOCIATED(cellml%MODELS_FIELD)) THEN
2103  CALL field_dof_order_type_get(cellml%MODELS_FIELD%MODELS_FIELD,field_u_variable_type,dof_order_type,err,error,*999)
2104  IF(dof_order_type==field_separated_component_dof_order) THEN
2105  !Dof components are separated. Will need to copy data to temporary arrays.
2106  IF(only_one_model_index==cellml_models_field_not_constant) THEN
2107  !Mulitple models
2108  DO WHILE(time<=end_time)
2109  DO dof_idx=1,n
2110  model_idx=models_data(dof_idx)
2111  IF(model_idx.GT.0) THEN
2112  model=>cellml%MODELS(model_idx)%PTR
2113  IF(ASSOCIATED(model)) THEN
2114  number_states=model%NUMBER_OF_STATE
2115  number_intermediates=model%NUMBER_OF_INTERMEDIATE
2116  number_parameters=model%NUMBER_OF_PARAMETERS
2117 
2118  !Copy CellML data to temporary arrays
2119  DO state_idx=1,number_states
2120  states(state_idx)=state_data((dof_idx-1)*n+state_idx)
2121  ENDDO !state_idx
2122  DO parameter_idx=1,number_parameters
2123  parameters(parameter_idx)=parameters_data((dof_idx-1)*n+parameter_idx)
2124  ENDDO !parameter_idx
2125 
2126 #ifdef WITH_CELLML
2127  CALL cellml_model_definition_call_rhs_routine(model%PTR,time,states,rates,intermediates, &
2128  & parameters)
2129 #else
2130  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
2131 #endif
2132 
2133  !Copy temporary data back to CellML arrays
2134  DO intermediate_idx=1,number_intermediates
2135  intermediate_data((dof_idx-1)*n+intermediate_idx)=intermediates(intermediate_idx)
2136  ENDDO !intermediate_idx
2137  DO state_idx=1,number_states
2138  state_data((dof_idx-1)*n+state_idx)=states(state_idx)+time_increment*rates(state_idx)
2139  ENDDO !state_idx
2140 
2141  ELSE
2142  local_error="CellML environment model is not associated for model index "// &
2143  & trim(numbertovstring(only_one_model_index,"*",err,error))//" belonging to dof index "// &
2144  & trim(numbertovstring(dof_idx,"*",err,error))//"."
2145  CALL flagerror(local_error,err,error,*999)
2146  ENDIF
2147  ENDIF !model_idx
2148  ENDDO !dof_idx
2149  time=time+time_increment
2150  ENDDO !time
2151  ELSE
2152  !One one model is used.
2153  model=>cellml%MODELS(only_one_model_index)%PTR
2154  IF(ASSOCIATED(model)) THEN
2155  number_states=model%NUMBER_OF_STATE
2156  number_intermediates=model%NUMBER_OF_INTERMEDIATE
2157  number_parameters=model%NUMBER_OF_PARAMETERS
2158  time=start_time
2159  DO WHILE(time<=end_time)
2160  DO dof_idx=1,n
2161 
2162  model_idx=models_data(dof_idx)
2163  IF(model_idx.GT.0) THEN
2164  !Copy CellML data to temporary arrays
2165  DO state_idx=1,number_states
2166  states(state_idx)=state_data((dof_idx-1)*n+state_idx)
2167  ENDDO !state_idx
2168  DO parameter_idx=1,number_parameters
2169  parameters(parameter_idx)=parameters_data((dof_idx-1)*n+parameter_idx)
2170  ENDDO !parameter_idx
2171 
2172 #ifdef WITH_CELLML
2173  CALL cellml_model_definition_call_rhs_routine(model%PTR,time,states,rates,intermediates, &
2174  & parameters)
2175 #else
2176  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
2177 #endif
2178 
2179  !Copy temporary data back to CellML arrays
2180  DO intermediate_idx=1,number_intermediates
2181  intermediate_data((dof_idx-1)*n+intermediate_idx)=intermediates(intermediate_idx)
2182  ENDDO !intermediate_idx
2183  DO state_idx=1,number_states
2184  state_data((dof_idx-1)*n+state_idx)=states(state_idx)+time_increment*rates(state_idx)
2185  ENDDO !state_idx
2186  ENDIF !model_idx
2187  ENDDO !dof_idx
2188  time=time+time_increment
2189  ENDDO !time
2190  ELSE
2191  local_error="CellML environment model is not associated for model index "// &
2192  & trim(numbertovstring(only_one_model_index,"*",err,error))//"."
2193  CALL flagerror(local_error,err,error,*999)
2194  ENDIF
2195  ENDIF
2196  ELSE
2197  !Dof components are continguous. Can pass data directly.
2198  IF(only_one_model_index==cellml_models_field_not_constant) THEN
2199  !Mulitple models
2200  time=start_time
2201  DO WHILE(time<=end_time)
2202  DO dof_idx=1,n
2203  model_idx=models_data(dof_idx)
2204  IF(model_idx==0) THEN
2205  ! Do nothing- empty model index specified
2206  ELSE IF(model_idx > 0 .AND. model_idx <= cellml%NUMBER_OF_MODELS) THEN
2207  model=>cellml%MODELS(model_idx)%PTR
2208  IF(ASSOCIATED(model)) THEN
2209  number_states=model%NUMBER_OF_STATE
2210  number_intermediates=model%NUMBER_OF_INTERMEDIATE
2211  number_parameters=model%NUMBER_OF_PARAMETERS
2212 
2213 #ifdef WITH_CELLML
2214  !Call RHS. Note some models might not have state, rates, intermediate or parameter data so call accordingly
2215  !to avoid referencing null pointers
2216  IF(number_states>0) THEN
2217  IF(number_intermediates>0) THEN
2218  IF(number_parameters>0) THEN
2219  !We have states, intermediate and parameters for the model
2220 
2221  state_start_dof=(dof_idx-1)*max_number_states+1
2222  state_end_dof=state_start_dof+number_states-1
2223  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
2224  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
2225  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
2226  parameter_end_dof=parameter_start_dof+number_parameters-1
2227 
2228  CALL cellml_model_definition_call_rhs_routine(model%PTR,time,state_data(state_start_dof: &
2229  & state_end_dof),rates,intermediate_data(intermediate_start_dof:intermediate_end_dof), &
2230  & parameters_data(parameter_start_dof:parameter_end_dof))
2231 
2232  ELSE
2233  !We do not have parameters in the model
2234 
2235  state_start_dof=(dof_idx-1)*max_number_states+1
2236  state_end_dof=state_start_dof+number_states-1
2237  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
2238  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
2239 
2240  CALL cellml_model_definition_call_rhs_routine(model%PTR,time,state_data(state_start_dof: &
2241  & state_end_dof),rates,intermediate_data(intermediate_start_dof:intermediate_end_dof), &
2242  & parameters)
2243 
2244  ENDIF
2245  ELSE
2246  IF(number_parameters>0) THEN
2247  !We do not have intermediates in the model
2248  state_start_dof=(dof_idx-1)*max_number_states+1
2249  state_end_dof=state_start_dof+number_states-1
2250  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
2251  parameter_end_dof=parameter_start_dof+number_parameters-1
2252 
2253  CALL cellml_model_definition_call_rhs_routine(model%PTR,time,state_data(state_start_dof: &
2254  & state_end_dof),rates,intermediates,parameters_data(parameter_start_dof:parameter_end_dof))
2255 
2256  ELSE
2257  !We do not have intermediates or parameters in the model
2258  state_start_dof=(dof_idx-1)*max_number_states+1
2259  state_end_dof=state_start_dof+number_states-1
2260 
2261  CALL cellml_model_definition_call_rhs_routine(model%PTR,time,state_data(state_start_dof: &
2262  & state_end_dof),rates,intermediates,parameters)
2263 
2264  ENDIF
2265  ENDIF
2266  ELSE
2267  CALL flagerror("Invalid CellML model for integration - there are no states.",err,error,*999)
2268  ENDIF
2269 
2270 #else
2271  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
2272 #endif
2273  state_data(state_start_dof:state_end_dof)=state_data(state_start_dof:state_end_dof)+ &
2274  & time_increment*rates(1:number_states)
2275  ELSE
2276  local_error="CellML environment model is not associated for model index "// &
2277  & trim(numbertovstring(only_one_model_index,"*",err,error))//" belonging to dof index "// &
2278  & trim(numbertovstring(dof_idx,"*",err,error))//"."
2279  CALL flagerror(local_error,err,error,*999)
2280  ENDIF
2281  ELSE
2282  local_error="Invalid CellML model index: "// &
2283  & trim(numbertovstring(model_idx,"*",err,error))//". The specified index should be between 1 and "// &
2284  & trim(numbertovstring(cellml%NUMBER_OF_MODELS,"*",err,error))//"."
2285  CALL flagerror(local_error,err,error,*999)
2286  ENDIF
2287  ENDDO !dof_idx
2288  time=time+time_increment
2289  ENDDO !time
2290  ELSE
2291  !One one model is used.
2292  model=>cellml%MODELS(only_one_model_index)%PTR
2293  IF(ASSOCIATED(model)) THEN
2294  number_states=model%NUMBER_OF_STATE
2295  number_intermediates=model%NUMBER_OF_INTERMEDIATE
2296  number_parameters=model%NUMBER_OF_PARAMETERS
2297 #ifdef WITH_CELLML
2298 
2299  !Call RHS. Note some models might not have state, rates, intermediate or parameter data so call accordingly
2300  !to avoid referencing null pointers
2301  IF(number_states>0) THEN
2302  IF(number_intermediates>0) THEN
2303  IF(number_parameters>0) THEN
2304  !We have states, intermediate and parameters for the model
2305 
2306  time=start_time
2307  DO WHILE(time<=end_time)
2308  DO dof_idx=1,n
2309  model_idx=models_data(dof_idx)
2310  IF(model_idx.GT.0) THEN
2311 
2312  state_start_dof=(dof_idx-1)*max_number_states+1
2313  state_end_dof=state_start_dof+number_states-1
2314  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
2315  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
2316  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
2317  parameter_end_dof=parameter_start_dof+number_parameters-1
2318 
2319  CALL cellml_model_definition_call_rhs_routine(model%PTR,time, &
2320  & state_data(state_start_dof:state_end_dof), &
2321  & rates,intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters_data( &
2322  & parameter_start_dof:parameter_end_dof))
2323 
2324  state_data(state_start_dof:state_end_dof)=state_data(state_start_dof:state_end_dof)+ &
2325  & time_increment*rates(1:number_states)
2326  ENDIF !model_idx
2327  ENDDO !dof_idx
2328  time=time+time_increment
2329  ENDDO !time
2330  ELSE
2331  !We do not have parameters in the model
2332  time=start_time
2333  DO WHILE(time<=end_time)
2334  DO dof_idx=1,n
2335  model_idx=models_data(dof_idx)
2336  IF(model_idx.GT.0) THEN
2337 
2338  state_start_dof=(dof_idx-1)*max_number_states+1
2339  state_end_dof=state_start_dof+number_states-1
2340  intermediate_start_dof=(dof_idx-1)*max_number_intermediates+1
2341  intermediate_end_dof=intermediate_start_dof+number_intermediates-1
2342 
2343  CALL cellml_model_definition_call_rhs_routine(model%PTR,time, &
2344  & state_data(state_start_dof:state_end_dof), &
2345  & rates,intermediate_data(intermediate_start_dof:intermediate_end_dof),parameters)
2346 
2347  state_data(state_start_dof:state_end_dof)=state_data(state_start_dof:state_end_dof)+ &
2348  & time_increment*rates(1:number_states)
2349  ENDIF !model_idx
2350  ENDDO !dof_idx
2351  time=time+time_increment
2352  ENDDO !time
2353  ENDIF
2354  ELSE
2355  IF(number_parameters>0) THEN
2356  !We do not have intermediates in the model
2357 
2358  time=start_time
2359  DO WHILE(time<=end_time)
2360  DO dof_idx=1,n
2361  model_idx=models_data(dof_idx)
2362  IF(model_idx.GT.0) THEN
2363 
2364  state_start_dof=(dof_idx-1)*max_number_states+1
2365  state_end_dof=state_start_dof+number_states-1
2366  parameter_start_dof=(dof_idx-1)*max_number_parameters+1
2367  parameter_end_dof=parameter_start_dof+number_parameters-1
2368 
2369  CALL cellml_model_definition_call_rhs_routine(model%PTR,time, &
2370  & state_data(state_start_dof:state_end_dof), &
2371  & rates,intermediates,parameters_data(parameter_start_dof:parameter_end_dof))
2372 
2373  state_data(state_start_dof:state_end_dof)=state_data(state_start_dof:state_end_dof)+ &
2374  & time_increment*rates(1:number_states)
2375  ENDIF !model_idx
2376  ENDDO !dof_idx
2377  time=time+time_increment
2378  ENDDO !time
2379  ELSE
2380  !We do not have intermediates or parameters in the model
2381  time=start_time
2382  DO WHILE(time<=end_time)
2383  DO dof_idx=1,n
2384  model_idx=models_data(dof_idx)
2385  IF(model_idx.GT.0) THEN
2386 
2387  state_start_dof=(dof_idx-1)*max_number_states+1
2388  state_end_dof=state_start_dof+number_states-1
2389 
2390  CALL cellml_model_definition_call_rhs_routine(model%PTR,time, &
2391  & state_data(state_start_dof:state_end_dof), &
2392  & rates,intermediates,parameters)
2393 
2394  state_data(state_start_dof:state_end_dof)=state_data(state_start_dof:state_end_dof)+ &
2395  & time_increment*rates(1:number_states)
2396  ENDIF !model_idx
2397  ENDDO !dof_idx
2398  time=time+time_increment
2399  ENDDO !time
2400  ENDIF
2401  ENDIF
2402  ELSE
2403  CALL flagerror("Invalid CellML model for integration - there are no states.",err,error,*999)
2404  ENDIF
2405 
2406 #else
2407  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
2408 #endif
2409 
2410  ELSE
2411  local_error="CellML environment model is not associated for model index "// &
2412  & trim(numbertovstring(only_one_model_index,"*",err,error))//"."
2413  CALL flagerror(local_error,err,error,*999)
2414  ENDIF
2415  ENDIF
2416  ENDIF
2417  ELSE
2418  CALL flagerror("CellML environment models field is not associated.",err,error,*999)
2419  ENDIF
2420  ELSE
2421  CALL flagerror("CellML environment is not associated.",err,error,*999)
2422  ENDIF
2423  ELSE
2424  CALL flagerror("Forward Euler solver is not associated.",err,error,*999)
2425  ENDIF
2426 
2427  exits("SOLVER_DAE_EULER_FORWARD_INTEGRATE")
2428  RETURN
2429 999 errorsexits("SOLVER_DAE_EULER_FORWARD_INTEGRATE",err,error)
2430  RETURN 1
2431 
2432  END SUBROUTINE solver_dae_euler_forward_integrate
2433 
2434  !
2435  !================================================================================================================================
2436  !
2437 
2439  SUBROUTINE solver_dae_euler_forward_solve(FORWARD_EULER_SOLVER,ERR,ERROR,*)
2441  !Argument variables
2442  TYPE(forward_euler_dae_solver_type), POINTER :: FORWARD_EULER_SOLVER
2443  INTEGER(INTG), INTENT(OUT) :: ERR
2444  TYPE(varying_string), INTENT(OUT) :: ERROR
2445  !Local Variables
2446  INTEGER(INTG) :: cellml_idx
2447  INTEGER(INTG), POINTER :: MODELS_DATA(:)
2448  REAL(DP), POINTER :: INTERMEDIATE_DATA(:),PARAMETERS_DATA(:),STATE_DATA(:)
2449  TYPE(cellml_type), POINTER :: CELLML_ENVIRONMENT
2450  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
2451  TYPE(cellml_models_field_type), POINTER :: CELLML_MODELS_FIELD
2452  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
2453  TYPE(euler_dae_solver_type), POINTER :: EULER_SOLVER
2454  TYPE(field_variable_type), POINTER :: MODELS_VARIABLE
2455  TYPE(field_type), POINTER :: MODELS_FIELD,STATE_FIELD,PARAMETERS_FIELD,INTERMEDIATE_FIELD
2456  TYPE(solver_type), POINTER :: SOLVER
2457  TYPE(varying_string) :: LOCAL_ERROR
2458 
2459  enters("SOLVER_DAE_EULER_FORWARD_SOLVE",err,error,*999)
2460 
2461  NULLIFY(models_data)
2462  NULLIFY(intermediate_data)
2463  NULLIFY(parameters_data)
2464  NULLIFY(state_data)
2465  NULLIFY(models_variable)
2466  NULLIFY(models_field)
2467  NULLIFY(state_field)
2468  NULLIFY(parameters_field)
2469  NULLIFY(intermediate_field)
2470 
2471  IF(ASSOCIATED(forward_euler_solver)) THEN
2472  euler_solver=>forward_euler_solver%EULER_DAE_SOLVER
2473  IF(ASSOCIATED(euler_solver)) THEN
2474  dae_solver=>euler_solver%DAE_SOLVER
2475  IF(ASSOCIATED(dae_solver)) THEN
2476  solver=>dae_solver%SOLVER
2477  IF(ASSOCIATED(solver)) THEN
2478  cellml_equations=>solver%CELLML_EQUATIONS
2479  IF(ASSOCIATED(cellml_equations)) THEN
2480  DO cellml_idx=1,cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS
2481  cellml_environment=>cellml_equations%CELLML_ENVIRONMENTS(cellml_idx)%PTR
2482  IF(ASSOCIATED(cellml_environment)) THEN
2483  cellml_models_field=>cellml_environment%MODELS_FIELD
2484  IF(ASSOCIATED(cellml_models_field)) THEN
2485  models_field=>cellml_models_field%MODELS_FIELD
2486  IF(ASSOCIATED(models_field)) THEN
2487 
2488 !!TODO: Maybe move this getting of fields earlier up the DAE solver chain? For now keep here.
2489 
2490  !Make sure CellML fields have been updated to the current value of any mapped fields
2491  CALL cellml_field_to_cellml_update(cellml_environment,err,error,*999)
2492 
2493  CALL field_variable_get(models_field,field_u_variable_type,models_variable,err,error,*999)
2494  CALL field_parameter_set_data_get(models_field,field_u_variable_type,field_values_set_type, &
2495  & models_data,err,error,*999)
2496 
2497  !Get the state information if this environment has any.
2498  IF(ASSOCIATED(cellml_environment%STATE_FIELD)) THEN
2499  state_field=>cellml_environment%STATE_FIELD%STATE_FIELD
2500  IF(ASSOCIATED(state_field)) THEN
2501  CALL field_parameter_set_data_get(state_field,field_u_variable_type,field_values_set_type, &
2502  & state_data,err,error,*999)
2503  ENDIF
2504  ENDIF
2505 
2506  !Get the parameters information if this environment has any.
2507  IF(ASSOCIATED(cellml_environment%PARAMETERS_FIELD)) THEN
2508  parameters_field=>cellml_environment%PARAMETERS_FIELD%PARAMETERS_FIELD
2509  IF(ASSOCIATED(parameters_field)) THEN
2510  CALL field_parameter_set_data_get(parameters_field,field_u_variable_type,field_values_set_type, &
2511  & parameters_data,err,error,*999)
2512  ENDIF
2513  ENDIF
2514 
2515  !Get the intermediate information if this environment has any.
2516  IF(ASSOCIATED(cellml_environment%INTERMEDIATE_FIELD)) THEN
2517  intermediate_field=>cellml_environment%INTERMEDIATE_FIELD%INTERMEDIATE_FIELD
2518  IF(ASSOCIATED(intermediate_field)) THEN
2519  CALL field_parameter_set_data_get(intermediate_field,field_u_variable_type,field_values_set_type, &
2520  & intermediate_data,err,error,*999)
2521  ENDIF
2522  ENDIF
2523 
2524  !Integrate these CellML equations
2525  CALL solver_dae_euler_forward_integrate(forward_euler_solver,cellml_environment,models_variable% &
2526  & total_number_of_dofs,dae_solver%START_TIME,dae_solver%END_TIME,dae_solver%INITIAL_STEP, &
2527  & cellml_environment%MODELS_FIELD%ONLY_ONE_MODEL_INDEX,models_data,cellml_environment% &
2528  & maximum_number_of_state,state_data,cellml_environment%MAXIMUM_NUMBER_OF_PARAMETERS, &
2529  & parameters_data,cellml_environment%MAXIMUM_NUMBER_OF_INTERMEDIATE,intermediate_data,err,error,*999)
2530 
2531  !Restore field data
2532  CALL field_parameter_set_data_restore(models_field,field_u_variable_type,field_values_set_type, &
2533  & models_data,err,error,*999)
2534  IF(ASSOCIATED(state_field)) CALL field_parameter_set_data_restore(state_field,field_u_variable_type, &
2535  & field_values_set_type,state_data,err,error,*999)
2536  IF(ASSOCIATED(parameters_field)) CALL field_parameter_set_data_restore(parameters_field, &
2537  & field_u_variable_type,field_values_set_type,parameters_data,err,error,*999)
2538  IF(ASSOCIATED(intermediate_field)) CALL field_parameter_set_data_restore(intermediate_field, &
2539  & field_u_variable_type,field_values_set_type,intermediate_data,err,error,*999)
2540 
2541  !Make sure fields have been updated to the current value of any mapped CellML fields
2542  CALL cellml_cellml_to_field_update(cellml_environment,err,error,*999)
2543 
2544  ELSE
2545  local_error="The CellML models field is not associated for CellML index "// &
2546  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
2547  CALL flagerror(local_error,err,error,*999)
2548  ENDIF
2549  ELSE
2550  local_error="The CellML models field is not associated for CellML index "// &
2551  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
2552  CALL flagerror(local_error,err,error,*999)
2553  ENDIF
2554  ELSE
2555  local_error="The CellML enviroment is not associated for for CellML index "// &
2556  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
2557  CALL flagerror(local_error,err,error,*999)
2558  ENDIF
2559  ENDDO !cellml_idx
2560  ELSE
2561  CALL flagerror("Solver solver equations is not associated.",err,error,*999)
2562  ENDIF
2563  ELSE
2564  CALL flagerror("Solver is not associated.",err,error,*999)
2565  ENDIF
2566  ELSE
2567  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*999)
2568  ENDIF
2569  ELSE
2570  CALL flagerror("Euler differential-algebraic equation solver is not associated.",err,error,*999)
2571  ENDIF
2572  ELSE
2573  CALL flagerror("Forward Euler differential-algebraic equation solver is not associated.",err,error,*999)
2574  ENDIF
2575 
2576  exits("SOLVER_DAE_EULER_FORWARD_SOLVE")
2577  RETURN
2578 999 errorsexits("SOLVER_DAE_EULER_FORWARD_SOLVE",err,error)
2579  RETURN 1
2580 
2581  END SUBROUTINE solver_dae_euler_forward_solve
2582 
2583  !
2584  !================================================================================================================================
2585  !
2586 
2588  SUBROUTINE solver_dae_euler_improved_finalise(IMPROVED_EULER_SOLVER,ERR,ERROR,*)
2590  !Argument variables
2591  TYPE(improved_euler_dae_solver_type), POINTER :: IMPROVED_EULER_SOLVER
2592  INTEGER(INTG), INTENT(OUT) :: ERR
2593  TYPE(varying_string), INTENT(OUT) :: ERROR
2594  !Local Variables
2595 
2596  enters("SOLVER_DAE_EULER_IMPROVED_FINALISE",err,error,*999)
2597 
2598  IF(ASSOCIATED(improved_euler_solver)) THEN
2599  DEALLOCATE(improved_euler_solver)
2600  ENDIF
2601 
2602  exits("SOLVER_DAE_EULER_IMPROVED_FINALISE")
2603  RETURN
2604 999 errorsexits("SOLVER_DAE_EULER_IMPROVED_FINALISE",err,error)
2605  RETURN 1
2606 
2607  END SUBROUTINE solver_dae_euler_improved_finalise
2608 
2609  !
2610  !================================================================================================================================
2611  !
2612 
2614  SUBROUTINE solver_dae_euler_improved_initialise(EULER_DAE_SOLVER,ERR,ERROR,*)
2616  !Argument variables
2617  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
2618  INTEGER(INTG), INTENT(OUT) :: ERR
2619  TYPE(varying_string), INTENT(OUT) :: ERROR
2620  !Local Variables
2621  INTEGER(INTG) :: DUMMY_ERR
2622  TYPE(varying_string) :: DUMMY_ERROR
2623 
2624  enters("SOLVER_DAE_EULER_IMPROVED_INITIALISE",err,error,*998)
2625 
2626  IF(ASSOCIATED(euler_dae_solver)) THEN
2627  IF(ASSOCIATED(euler_dae_solver%IMPROVED_EULER_SOLVER)) THEN
2628  CALL flagerror("Improved Euler solver is already associated for this Euler differential-algebraic equation solver.", &
2629  & err,error,*998)
2630  ELSE
2631  !Allocate the improved Euler solver
2632  ALLOCATE(euler_dae_solver%IMPROVED_EULER_SOLVER,stat=err)
2633  IF(err/=0) CALL flagerror("Could not allocate improved Euler solver.",err,error,*999)
2634  !Initialise
2635  euler_dae_solver%IMPROVED_EULER_SOLVER%EULER_DAE_SOLVER=>euler_dae_solver
2636  euler_dae_solver%IMPROVED_EULER_SOLVER%SOLVER_LIBRARY=0
2637  !Defaults
2638  ENDIF
2639  ELSE
2640  CALL flagerror("Euler differential-algebraic equation solver is not associated.",err,error,*998)
2641  ENDIF
2642 
2643  exits("SOLVER_DAE_EULER_IMPROVED_INITIALISE")
2644  RETURN
2645 999 CALL solver_dae_euler_improved_finalise(euler_dae_solver%IMPROVED_EULER_SOLVER,dummy_err,dummy_error,*998)
2646 998 errorsexits("SOLVER_DAE_EULER_IMPROVED_INITIALISE",err,error)
2647  RETURN 1
2648 
2650 
2651  !
2652  !================================================================================================================================
2653  !
2654 
2656  SUBROUTINE solver_dae_euler_improved_solve(IMPROVED_EULER_SOLVER,ERR,ERROR,*)
2658  !Argument variables
2659  TYPE(improved_euler_dae_solver_type), POINTER :: IMPROVED_EULER_SOLVER
2660  INTEGER(INTG), INTENT(OUT) :: ERR
2661  TYPE(varying_string), INTENT(OUT) :: ERROR
2662  !Local Variables
2663 
2664  enters("SOLVER_DAE_EULER_IMPROVED_SOLVE",err,error,*999)
2665 
2666  IF(ASSOCIATED(improved_euler_solver)) THEN
2667  CALL flagerror("Not implemented.",err,error,*999)
2668  ELSE
2669  CALL flagerror("Improved Euler differential-algebraic equation solver is not associated.",err,error,*999)
2670  ENDIF
2671 
2672  exits("SOLVER_DAE_EULER_IMPROVED_SOLVE")
2673  RETURN
2674 999 errorsexits("SOLVER_DAE_EULER_IMPROVED_SOLVE",err,error)
2675  RETURN 1
2676 
2677  END SUBROUTINE solver_dae_euler_improved_solve
2678 
2679  !
2680  !================================================================================================================================
2681  !
2682 
2684  SUBROUTINE solver_dae_euler_initialise(DAE_SOLVER,ERR,ERROR,*)
2686  !Argument variables
2687  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
2688  INTEGER(INTG), INTENT(OUT) :: ERR
2689  TYPE(varying_string), INTENT(OUT) :: ERROR
2690  !Local Variables
2691  INTEGER(INTG) :: DUMMY_ERR
2692  TYPE(varying_string) :: DUMMY_ERROR
2693 
2694  enters("SOLVER_DAE_EULER_INITIALISE",err,error,*998)
2695 
2696  IF(ASSOCIATED(dae_solver)) THEN
2697  IF(ASSOCIATED(dae_solver%EULER_SOLVER)) THEN
2698  CALL flagerror("Euler solver is already associated for this differential-algebraic equation solver.",err,error,*998)
2699  ELSE
2700  !Allocate the Euler solver
2701  ALLOCATE(dae_solver%EULER_SOLVER,stat=err)
2702  IF(err/=0) CALL flagerror("Could not allocate Euler solver.",err,error,*999)
2703  !Initialise
2704  dae_solver%EULER_SOLVER%DAE_SOLVER=>dae_solver
2705  NULLIFY(dae_solver%EULER_SOLVER%FORWARD_EULER_SOLVER)
2706  NULLIFY(dae_solver%EULER_SOLVER%BACKWARD_EULER_SOLVER)
2707  NULLIFY(dae_solver%EULER_SOLVER%IMPROVED_EULER_SOLVER)
2708  !Default to a forward Euler solver
2709  CALL solver_dae_euler_forward_initialise(dae_solver%EULER_SOLVER,err,error,*999)
2710  dae_solver%EULER_SOLVER%EULER_TYPE=solver_dae_euler_forward
2711  ENDIF
2712  ELSE
2713  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*998)
2714  ENDIF
2715 
2716  exits("SOLVER_DAE_EULER_INITIALISE")
2717  RETURN
2718 999 CALL solver_dae_euler_finalise(dae_solver%EULER_SOLVER,dummy_err,dummy_error,*998)
2719 998 errorsexits("SOLVER_DAE_EULER_INITIALISE",err,error)
2720  RETURN 1
2721 
2722  END SUBROUTINE solver_dae_euler_initialise
2723 
2724  !
2725  !================================================================================================================================
2726  !
2727 
2729  SUBROUTINE solver_dae_euler_library_type_get(EULER_DAE_SOLVER,SOLVER_LIBRARY_TYPE,ERR,ERROR,*)
2731  !Argument variables
2732  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
2733  INTEGER(INTG), INTENT(OUT) :: SOLVER_LIBRARY_TYPE
2734  INTEGER(INTG), INTENT(OUT) :: ERR
2735  TYPE(varying_string), INTENT(OUT) :: ERROR
2736  !Local Variables
2737  TYPE(backward_euler_dae_solver_type), POINTER :: BACKWARD_EULER_DAE_SOLVER
2738  TYPE(forward_euler_dae_solver_type), POINTER :: FORWARD_EULER_DAE_SOLVER
2739  TYPE(improved_euler_dae_solver_type), POINTER :: IMPROVED_EULER_DAE_SOLVER
2740  TYPE(varying_string) :: LOCAL_ERROR
2741 
2742  enters("SOLVER_DAE_EULER_LIBRARY_TYPE_GET",err,error,*999)
2743 
2744  IF(ASSOCIATED(euler_dae_solver)) THEN
2745  SELECT CASE(euler_dae_solver%EULER_TYPE)
2747  forward_euler_dae_solver=>euler_dae_solver%FORWARD_EULER_SOLVER
2748  IF(ASSOCIATED(forward_euler_dae_solver)) THEN
2749  solver_library_type=forward_euler_dae_solver%SOLVER_LIBRARY
2750  ELSE
2751  CALL flagerror("The forward Euler differntial-algebraic equations solver is not associated.",err,error,*999)
2752  ENDIF
2754  backward_euler_dae_solver=>euler_dae_solver%BACKWARD_EULER_SOLVER
2755  IF(ASSOCIATED(backward_euler_dae_solver)) THEN
2756  solver_library_type=backward_euler_dae_solver%SOLVER_LIBRARY
2757  ELSE
2758  CALL flagerror("The backward Euler differntial-algebraic equations solver is not associated.",err,error,*999)
2759  ENDIF
2761  improved_euler_dae_solver=>euler_dae_solver%IMPROVED_EULER_SOLVER
2762  IF(ASSOCIATED(improved_euler_dae_solver)) THEN
2763  solver_library_type=improved_euler_dae_solver%SOLVER_LIBRARY
2764  ELSE
2765  CALL flagerror("The improved Euler differntial-algebraic equations solver is not associated.",err,error,*999)
2766  ENDIF
2767  CASE DEFAULT
2768  local_error="The Euler differential-algebraic equations solver type of "// &
2769  & trim(numbertovstring(euler_dae_solver%EULER_TYPE,"*",err,error))//" is invalid."
2770  CALL flagerror(local_error,err,error,*999)
2771  END SELECT
2772  ELSE
2773  CALL flagerror("Euler DAE solver is not associated.",err,error,*999)
2774  ENDIF
2775 
2776  exits("SOLVER_DAE_EULER_LIBRARY_TYPE_GET")
2777  RETURN
2778 999 errorsexits("SOLVER_DAE_EULER_LIBRARY_TYPE_GET",err,error)
2779  RETURN 1
2780 
2781  END SUBROUTINE solver_dae_euler_library_type_get
2782 
2783  !
2784  !================================================================================================================================
2785  !
2786 
2788  SUBROUTINE solver_dae_euler_library_type_set(EULER_DAE_SOLVER,SOLVER_LIBRARY_TYPE,ERR,ERROR,*)
2790  !Argument variables
2791  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
2792  INTEGER(INTG), INTENT(IN) :: SOLVER_LIBRARY_TYPE
2793  INTEGER(INTG), INTENT(OUT) :: ERR
2794  TYPE(varying_string), INTENT(OUT) :: ERROR
2795  !Local Variables
2796  TYPE(backward_euler_dae_solver_type), POINTER :: BACKWARD_EULER_DAE_SOLVER
2797  TYPE(forward_euler_dae_solver_type), POINTER :: FORWARD_EULER_DAE_SOLVER
2798  TYPE(improved_euler_dae_solver_type), POINTER :: IMPROVED_EULER_DAE_SOLVER
2799  TYPE(varying_string) :: LOCAL_ERROR
2800 
2801  enters("SOLVER_DAE_EULER_LIBRARY_TYPE_SET",err,error,*999)
2802 
2803  IF(ASSOCIATED(euler_dae_solver)) THEN
2804  SELECT CASE(euler_dae_solver%EULER_TYPE)
2806  forward_euler_dae_solver=>euler_dae_solver%FORWARD_EULER_SOLVER
2807  IF(ASSOCIATED(forward_euler_dae_solver)) THEN
2808  SELECT CASE(solver_library_type)
2809  CASE(solver_cmiss_library)
2810  CALL flagerror("Not implemented.",err,error,*999)
2811  CASE(solver_petsc_library)
2812  CALL flagerror("Not implemented.",err,error,*999)
2813  CASE DEFAULT
2814  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
2815  & " is invalid for a forward Euler DAE solver."
2816  CALL flagerror(local_error,err,error,*999)
2817  END SELECT
2818  ELSE
2819  CALL flagerror("The forward Euler differential-algebraic equation solver is not associated.",err,error,*999)
2820  ENDIF
2822  backward_euler_dae_solver=>euler_dae_solver%BACKWARD_EULER_SOLVER
2823  IF(ASSOCIATED(backward_euler_dae_solver)) THEN
2824  SELECT CASE(solver_library_type)
2825  CASE(solver_cmiss_library)
2826  CALL flagerror("Not implemented.",err,error,*999)
2827  CASE(solver_petsc_library)
2828  CALL flagerror("Not implemented.",err,error,*999)
2829  CASE DEFAULT
2830  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
2831  & " is invalid for a backward Euler DAE solver."
2832  CALL flagerror(local_error,err,error,*999)
2833  END SELECT
2834  ELSE
2835  CALL flagerror("The backward Euler differential-algebraic equation solver is not associated.",err,error,*999)
2836  ENDIF
2838  improved_euler_dae_solver=>euler_dae_solver%IMPROVED_EULER_SOLVER
2839  IF(ASSOCIATED(improved_euler_dae_solver)) THEN
2840  SELECT CASE(solver_library_type)
2841  CASE(solver_cmiss_library)
2842  CALL flagerror("Not implemented.",err,error,*999)
2843  CASE(solver_petsc_library)
2844  CALL flagerror("Not implemented.",err,error,*999)
2845  CASE DEFAULT
2846  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
2847  & " is invalid for an improved Euler DAE solver."
2848  CALL flagerror(local_error,err,error,*999)
2849  END SELECT
2850  ELSE
2851  CALL flagerror("The improved Euler differential-algebraic equation solver is not associated.",err,error,*999)
2852  ENDIF
2853  CASE DEFAULT
2854  local_error="The Euler differential-algebraic equations solver type of "// &
2855  & trim(numbertovstring(euler_dae_solver%EULER_TYPE,"*",err,error))//" is invalid."
2856  CALL flagerror(local_error,err,error,*999)
2857  END SELECT
2858  ELSE
2859  CALL flagerror("The Euler differential-algebraic equation solver is not associated.",err,error,*999)
2860  ENDIF
2861 
2862  exits("SOLVER_DAE_EULER_LIBRARY_TYPE_SET")
2863  RETURN
2864 999 errorsexits("SOLVER_DAE_EULER_LIBRARY_TYPE_SET",err,error)
2865  RETURN 1
2866 
2867  END SUBROUTINE solver_dae_euler_library_type_set
2868 
2869  !
2870  !================================================================================================================================
2871  !
2872 
2874  SUBROUTINE solver_dae_euler_solve(EULER_SOLVER,ERR,ERROR,*)
2876  !Argument variables
2877  TYPE(euler_dae_solver_type), POINTER :: EULER_SOLVER
2878  INTEGER(INTG), INTENT(OUT) :: ERR
2879  TYPE(varying_string), INTENT(OUT) :: ERROR
2880  !Local Variables
2881  TYPE(varying_string) :: LOCAL_ERROR
2882 
2883  enters("SOLVER_DAE_EULER_SOLVE",err,error,*999)
2884 
2885  IF(ASSOCIATED(euler_solver)) THEN
2886  SELECT CASE(euler_solver%EULER_TYPE)
2888  CALL solver_dae_euler_forward_solve(euler_solver%FORWARD_EULER_SOLVER,err,error,*999)
2890  CALL solver_dae_euler_backward_solve(euler_solver%BACKWARD_EULER_SOLVER,err,error,*999)
2892  CALL solver_dae_euler_improved_solve(euler_solver%IMPROVED_EULER_SOLVER,err,error,*999)
2893  CASE DEFAULT
2894  local_error="The Euler differential-algebraic equation solver type of "// &
2895  & trim(numbertovstring(euler_solver%EULER_TYPE,"*",err,error))//" is invalid."
2896  CALL flagerror(local_error,err,error,*999)
2897  END SELECT
2898  ELSE
2899  CALL flagerror("Euler differential-algebraic equation solver is not associated.",err,error,*999)
2900  ENDIF
2901 
2902  exits("SOLVER_DAE_EULER_SOLVE")
2903  RETURN
2904 999 errorsexits("SOLVER_DAE_EULER_SOLVE",err,error)
2905  RETURN 1
2906 
2907  END SUBROUTINE solver_dae_euler_solve
2908 
2909  !
2910  !================================================================================================================================
2911  !
2912 
2914  SUBROUTINE solver_dae_euler_solver_type_get(SOLVER,DAE_EULER_TYPE,ERR,ERROR,*)
2916  !Argument variables
2917  TYPE(solver_type), POINTER :: SOLVER
2918  INTEGER(INTG), INTENT(OUT) :: DAE_EULER_TYPE
2919  INTEGER(INTG), INTENT(OUT) :: ERR
2920  TYPE(varying_string), INTENT(OUT) :: ERROR
2921  !Local Variables
2922  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
2923  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
2924 
2925  enters("SOLVER_DAE_EULER_SOLVER_TYPE_GET",err,error,*999)
2926 
2927  IF(ASSOCIATED(solver)) THEN
2928  IF(solver%SOLVER_FINISHED) THEN
2929  IF(solver%SOLVE_TYPE==solver_dae_type) THEN
2930  dae_solver=>solver%DAE_SOLVER
2931  IF(ASSOCIATED(dae_solver)) THEN
2932  IF(dae_solver%DAE_SOLVE_TYPE==solver_dae_euler) THEN
2933  euler_dae_solver=>dae_solver%EULER_SOLVER
2934  IF(ASSOCIATED(euler_dae_solver)) THEN
2935  dae_euler_type=euler_dae_solver%EULER_TYPE
2936  ELSE
2937  CALL flagerror("The differential-algebraic equation solver Euler solver is not associated.",err,error,*999)
2938  ENDIF
2939  ELSE
2940  CALL flagerror("The solver differential-algebraic equation solver is not an Euler differential-algebraic "// &
2941  & "equation solver.",err,error,*999)
2942  ENDIF
2943  ELSE
2944  CALL flagerror("The solver differential-algebraic equation solver is not associated.",err,error,*999)
2945  ENDIF
2946  ELSE
2947  CALL flagerror("The solver is not a differential-algebraic equation solver.",err,error,*999)
2948  ENDIF
2949  ELSE
2950  CALL flagerror("Solver has not been finished.",err,error,*999)
2951  ENDIF
2952  ELSE
2953  CALL flagerror("Solver is not associated.",err,error,*999)
2954  ENDIF
2955 
2956  exits("SOLVER_DAE_EULER_SOLVER_TYPE_GET")
2957  RETURN
2958 999 errorsexits("SOLVER_DAE_EULER_SOLVER_TYPE_GET",err,error)
2959  RETURN 1
2960 
2961  END SUBROUTINE solver_dae_euler_solver_type_get
2962 
2963  !
2964  !================================================================================================================================
2965  !
2966 
2968  SUBROUTINE solver_dae_euler_solver_type_set(SOLVER,DAE_EULER_TYPE,ERR,ERROR,*)
2970  !Argument variables
2971  TYPE(solver_type), POINTER :: SOLVER
2972  INTEGER(INTG), INTENT(IN) :: DAE_EULER_TYPE
2973  INTEGER(INTG), INTENT(OUT) :: ERR
2974  TYPE(varying_string), INTENT(OUT) :: ERROR
2975  !Local Variables
2976  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
2977  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
2978  TYPE(varying_string) :: LOCAL_ERROR
2979 
2980  enters("SOLVER_DAE_EULER_SOLVER_TYPE_SET",err,error,*999)
2981 
2982  IF(ASSOCIATED(solver)) THEN
2983  IF(solver%SOLVER_FINISHED) THEN
2984  CALL flagerror("Solver has already been finished.",err,error,*999)
2985  ELSE
2986  IF(solver%SOLVE_TYPE==solver_dae_type) THEN
2987  dae_solver=>solver%DAE_SOLVER
2988  IF(ASSOCIATED(dae_solver)) THEN
2989  IF(dae_solver%DAE_SOLVE_TYPE==solver_dae_euler) THEN
2990  euler_dae_solver=>dae_solver%EULER_SOLVER
2991  IF(ASSOCIATED(euler_dae_solver)) THEN
2992  IF(dae_euler_type/=euler_dae_solver%EULER_TYPE) THEN
2993  !Intialise the new Euler differential-algebraic equation solver type
2994  SELECT CASE(dae_euler_type)
2996  CALL solver_dae_euler_forward_initialise(euler_dae_solver,err,error,*999)
2998  CALL solver_dae_euler_backward_initialise(euler_dae_solver,err,error,*999)
3000  CALL solver_dae_euler_improved_initialise(euler_dae_solver,err,error,*999)
3001  CASE DEFAULT
3002  local_error="The specified Euler differential-algebraic equation solver type of "// &
3003  & trim(numbertovstring(dae_euler_type,"*",err,error))//" is invalid."
3004  CALL flagerror(local_error,err,error,*999)
3005  END SELECT
3006  !Finalise the old Euler differential-algebraic equation solver type
3007  SELECT CASE(euler_dae_solver%EULER_TYPE)
3009  CALL solver_dae_euler_forward_finalise(euler_dae_solver%FORWARD_EULER_SOLVER,err,error,*999)
3011  CALL solver_dae_euler_backward_finalise(euler_dae_solver%BACKWARD_EULER_SOLVER,err,error,*999)
3013  CALL solver_dae_euler_improved_finalise(euler_dae_solver%IMPROVED_EULER_SOLVER,err,error,*999)
3014  CASE DEFAULT
3015  local_error="The Euler differential-algebraic equation solver type of "// &
3016  & trim(numbertovstring(euler_dae_solver%EULER_TYPE,"*",err,error))//" is invalid."
3017  CALL flagerror(local_error,err,error,*999)
3018  END SELECT
3019  euler_dae_solver%EULER_TYPE=dae_euler_type
3020  ENDIF
3021  ELSE
3022  CALL flagerror("The differential-algebraic equation solver Euler solver is not associated.",err,error,*999)
3023  ENDIF
3024  ELSE
3025  CALL flagerror("The solver differential-algebraic equation solver is not an Euler differential-algebraic "// &
3026  & "equation solver.",err,error,*999)
3027  ENDIF
3028  ELSE
3029  CALL flagerror("The solver differential-algebraic equation solver is not associated.",err,error,*999)
3030  ENDIF
3031  ELSE
3032  CALL flagerror("The solver is not a differential-algebraic equation solver.",err,error,*999)
3033  ENDIF
3034  ENDIF
3035  ELSE
3036  CALL flagerror("Solver is not associated.",err,error,*999)
3037  ENDIF
3038 
3039  exits("SOLVER_DAE_EULER_SOLVER_TYPE_SET")
3040  RETURN
3041 999 errorsexits("SOLVER_DAE_EULER_SOLVER_TYPE_SET",err,error)
3042  RETURN 1
3043 
3044  END SUBROUTINE solver_dae_euler_solver_type_set
3045 
3046  !
3047  !================================================================================================================================
3048  !
3049 
3051  SUBROUTINE solver_dae_finalise(DAE_SOLVER,ERR,ERROR,*)
3053  !Argument variables
3054  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
3055  INTEGER(INTG), INTENT(OUT) :: ERR
3056  TYPE(varying_string), INTENT(OUT) :: ERROR
3057  !Local Variables
3058 
3059  enters("SOLVER_DAE_FINALISE",err,error,*999)
3060 
3061  IF(ASSOCIATED(dae_solver)) THEN
3062  CALL solver_dae_euler_finalise(dae_solver%EULER_SOLVER,err,error,*999)
3063  CALL solver_dae_crank_nicolson_finalise(dae_solver%CRANK_NICOLSON_SOLVER,err,error,*999)
3064  CALL solver_dae_runge_kutta_finalise(dae_solver%RUNGE_KUTTA_SOLVER,err,error,*999)
3065  CALL solver_dae_adams_moulton_finalise(dae_solver%ADAMS_MOULTON_SOLVER,err,error,*999)
3066  CALL solver_dae_bdf_finalise(dae_solver%BDF_SOLVER,err,error,*999)
3067  CALL solver_dae_rush_larson_finalise(dae_solver%RUSH_LARSON_SOLVER,err,error,*999)
3068  CALL solver_dae_external_finalise(dae_solver%EXTERNAL_SOLVER,err,error,*999)
3069  DEALLOCATE(dae_solver)
3070  ENDIF
3071 
3072  exits("SOLVER_DAE_FINALISE")
3073  RETURN
3074 999 errorsexits("SOLVER_DAE_FINALISE",err,error)
3075  RETURN 1
3076 
3077  END SUBROUTINE solver_dae_finalise
3078 
3079  !
3080  !================================================================================================================================
3081  !
3082 
3084  SUBROUTINE solver_dae_initialise(SOLVER,ERR,ERROR,*)
3086  !Argument variables
3087  TYPE(solver_type), POINTER :: SOLVER
3088  INTEGER(INTG), INTENT(OUT) :: ERR
3089  TYPE(varying_string), INTENT(OUT) :: ERROR
3090  !Local Variables
3091  INTEGER(INTG) :: DUMMY_ERR
3092  TYPE(varying_string) :: DUMMY_ERROR
3093 
3094  enters("SOLVER_DAE_INITIALISE",err,error,*998)
3095 
3096  IF(ASSOCIATED(solver)) THEN
3097  IF(ASSOCIATED(solver%DAE_SOLVER)) THEN
3098  CALL flagerror("Differential-algebraic equation solver is already associated for this solver.",err,error,*998)
3099  ELSE
3100  !Allocate the differential-algebraic equation solver
3101  ALLOCATE(solver%DAE_SOLVER,stat=err)
3102  IF(err/=0) CALL flagerror("Could not allocate solver differential-algebraic equation solver.",err,error,*999)
3103  !Initialise
3104  solver%DAE_SOLVER%SOLVER=>solver
3105  solver%DAE_SOLVER%DAE_TYPE=0
3106  solver%DAE_SOLVER%DAE_SOLVE_TYPE=0
3107  solver%DAE_SOLVER%START_TIME=0.0_dp
3108  solver%DAE_SOLVER%END_TIME=0.1_dp
3109  solver%DAE_SOLVER%INITIAL_STEP=0.1_dp
3110  NULLIFY(solver%DAE_SOLVER%EULER_SOLVER)
3111  NULLIFY(solver%DAE_SOLVER%CRANK_NICOLSON_SOLVER)
3112  NULLIFY(solver%DAE_SOLVER%RUNGE_KUTTA_SOLVER)
3113  NULLIFY(solver%DAE_SOLVER%ADAMS_MOULTON_SOLVER)
3114  NULLIFY(solver%DAE_SOLVER%BDF_SOLVER)
3115  NULLIFY(solver%DAE_SOLVER%RUSH_LARSON_SOLVER)
3116  NULLIFY(solver%DAE_SOLVER%EXTERNAL_SOLVER)
3117  !Default to an Euler differential equation solver
3118  CALL solver_dae_euler_initialise(solver%DAE_SOLVER,err,error,*999)
3119  solver%DAE_SOLVER%DAE_SOLVE_TYPE=solver_dae_euler
3120  ENDIF
3121  ELSE
3122  CALL flagerror("Solver is not associated.",err,error,*998)
3123  ENDIF
3124 
3125  exits("SOLVER_DAE_INITIALISE")
3126  RETURN
3127 999 CALL solver_dae_finalise(solver%DAE_SOLVER,dummy_err,dummy_error,*998)
3128 998 errorsexits("SOLVER_DAE_INITIALISE",err,error)
3129  RETURN 1
3130 
3131  END SUBROUTINE solver_dae_initialise
3132 
3133  !
3134  !================================================================================================================================
3135  !
3136 
3138  SUBROUTINE solver_dae_library_type_get(DAE_SOLVER,SOLVER_LIBRARY_TYPE,ERR,ERROR,*)
3140  !Argument variables
3141  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
3142  INTEGER(INTG), INTENT(OUT) :: SOLVER_LIBRARY_TYPE
3143  INTEGER(INTG), INTENT(OUT) :: ERR
3144  TYPE(varying_string), INTENT(OUT) :: ERROR
3145  !Local Variables
3146  TYPE(adams_moulton_dae_solver_type), POINTER :: ADAMS_MOULTON_DAE_SOLVER
3147  TYPE(bdf_dae_solver_type), POINTER :: BDF_DAE_SOLVER
3148  TYPE(crank_nicolson_dae_solver_type), POINTER :: CRANK_NICOLSON_DAE_SOLVER
3149  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
3150  TYPE(runge_kutta_dae_solver_type), POINTER :: RUNGE_KUTTA_DAE_SOLVER
3151  TYPE(rush_larson_dae_solver_type), POINTER :: RUSH_LARSON_DAE_SOLVER
3152  TYPE(varying_string) :: LOCAL_ERROR
3153 
3154  enters("SOLVER_DAE_LIBRARY_TYPE_GET",err,error,*999)
3155 
3156  IF(ASSOCIATED(dae_solver)) THEN
3157  SELECT CASE(dae_solver%DAE_SOLVE_TYPE)
3158  CASE(solver_dae_euler)
3159  euler_dae_solver=>dae_solver%EULER_SOLVER
3160  IF(ASSOCIATED(euler_dae_solver)) THEN
3161  CALL solver_dae_euler_library_type_get(euler_dae_solver,solver_library_type,err,error,*999)
3162  ELSE
3163  CALL flagerror("Euler differential-algebraic solver is not associated.",err,error,*999)
3164  ENDIF
3166  crank_nicolson_dae_solver=>dae_solver%CRANK_NICOLSON_SOLVER
3167  IF(ASSOCIATED(crank_nicolson_dae_solver)) THEN
3168  solver_library_type=crank_nicolson_dae_solver%SOLVER_LIBRARY
3169  ELSE
3170  CALL flagerror("The Crank-Nicolson differntial-algebraic equations solver is not associated.",err,error,*999)
3171  ENDIF
3173  runge_kutta_dae_solver=>dae_solver%RUNGE_KUTTA_SOLVER
3174  IF(ASSOCIATED(runge_kutta_dae_solver)) THEN
3175  solver_library_type=runge_kutta_dae_solver%SOLVER_LIBRARY
3176  ELSE
3177  CALL flagerror("The Runge-Kutta differntial-algebraic equations solver is not associated.",err,error,*999)
3178  ENDIF
3180  adams_moulton_dae_solver=>dae_solver%ADAMS_MOULTON_SOLVER
3181  IF(ASSOCIATED(adams_moulton_dae_solver)) THEN
3182  solver_library_type=adams_moulton_dae_solver%SOLVER_LIBRARY
3183  ELSE
3184  CALL flagerror("The Adams-Moulton differntial-algebraic equations solver is not associated.",err,error,*999)
3185  ENDIF
3186  CASE(solver_dae_bdf)
3187  bdf_dae_solver=>dae_solver%BDF_SOLVER
3188  IF(ASSOCIATED(bdf_dae_solver)) THEN
3189  solver_library_type=bdf_dae_solver%SOLVER_LIBRARY
3190  ELSE
3191  CALL flagerror("The BDF differntial-algebraic equations solver is not associated.",err,error,*999)
3192  ENDIF
3194  rush_larson_dae_solver=>dae_solver%RUSH_LARSON_SOLVER
3195  IF(ASSOCIATED(rush_larson_dae_solver)) THEN
3196  solver_library_type=rush_larson_dae_solver%SOLVER_LIBRARY
3197  ELSE
3198  CALL flagerror("The Rush-Larson differntial-algebraic equations solver is not associated.",err,error,*999)
3199  ENDIF
3200  CASE(solver_dae_external)
3201  CALL flagerror("Can not get the solver library for an external differntial-algebraic equations solver.",err,error,*999)
3202  CASE DEFAULT
3203  local_error="The differential-algebraic equations solver type of "// &
3204  & trim(numbertovstring(dae_solver%DAE_SOLVE_TYPE,"*",err,error))//" is invalid."
3205  CALL flagerror(local_error,err,error,*999)
3206  END SELECT
3207  ELSE
3208  CALL flagerror("DAE solver is not associated.",err,error,*999)
3209  ENDIF
3210 
3211  exits("SOLVER_DAE_LIBRARY_TYPE_GET")
3212  RETURN
3213 999 errorsexits("SOLVER_DAE_LIBRARY_TYPE_GET",err,error)
3214  RETURN 1
3215 
3216  END SUBROUTINE solver_dae_library_type_get
3217 
3218  !
3219  !================================================================================================================================
3220  !
3221 
3223  SUBROUTINE solver_dae_library_type_set(DAE_SOLVER,SOLVER_LIBRARY_TYPE,ERR,ERROR,*)
3225  !Argument variables
3226  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
3227  INTEGER(INTG), INTENT(IN) :: SOLVER_LIBRARY_TYPE
3228  INTEGER(INTG), INTENT(OUT) :: ERR
3229  TYPE(varying_string), INTENT(OUT) :: ERROR
3230  !Local Variables
3231  TYPE(adams_moulton_dae_solver_type), POINTER :: ADAMS_MOULTON_DAE_SOLVER
3232  TYPE(backward_euler_dae_solver_type), POINTER :: BACKWARD_EULER_DAE_SOLVER
3233  TYPE(bdf_dae_solver_type), POINTER :: BDF_DAE_SOLVER
3234  TYPE(crank_nicolson_dae_solver_type), POINTER :: CRANK_NICOLSON_DAE_SOLVER
3235  TYPE(euler_dae_solver_type), POINTER :: EULER_DAE_SOLVER
3236  TYPE(forward_euler_dae_solver_type), POINTER :: FORWARD_EULER_DAE_SOLVER
3237  TYPE(improved_euler_dae_solver_type), POINTER :: IMPROVED_EULER_DAE_SOLVER
3238  TYPE(runge_kutta_dae_solver_type), POINTER :: RUNGE_KUTTA_DAE_SOLVER
3239  TYPE(rush_larson_dae_solver_type), POINTER :: RUSH_LARSON_DAE_SOLVER
3240  TYPE(varying_string) :: LOCAL_ERROR
3241 
3242  enters("SOLVER_DAE_LIBRARY_TYPE_SET",err,error,*999)
3243 
3244  IF(ASSOCIATED(dae_solver)) THEN
3245  SELECT CASE(dae_solver%DAE_SOLVE_TYPE)
3246  CASE(solver_dae_euler)
3247  euler_dae_solver=>dae_solver%EULER_SOLVER
3248  IF(ASSOCIATED(euler_dae_solver)) THEN
3249  SELECT CASE(euler_dae_solver%EULER_TYPE)
3251  forward_euler_dae_solver=>euler_dae_solver%FORWARD_EULER_SOLVER
3252  IF(ASSOCIATED(forward_euler_dae_solver)) THEN
3253  SELECT CASE(solver_library_type)
3254  CASE(solver_cmiss_library)
3255  forward_euler_dae_solver%SOLVER_LIBRARY=solver_cmiss_library
3256  CASE(solver_petsc_library)
3257  CALL flagerror("Not implemented.",err,error,*999)
3258  CASE DEFAULT
3259  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3260  & " is invalid."
3261  CALL flagerror(local_error,err,error,*999)
3262  END SELECT
3263  ELSE
3264  CALL flagerror("The forward Euler differential-algebraic equation solver is not associated.",err,error,*999)
3265  ENDIF
3267  backward_euler_dae_solver=>euler_dae_solver%BACKWARD_EULER_SOLVER
3268  IF(ASSOCIATED(backward_euler_dae_solver)) THEN
3269  SELECT CASE(solver_library_type)
3270  CASE(solver_cmiss_library)
3271  CALL flagerror("Not implemented.",err,error,*999)
3272  CASE(solver_petsc_library)
3273  CALL flagerror("Not implemented.",err,error,*999)
3274  CASE DEFAULT
3275  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3276  & " is invalid."
3277  CALL flagerror(local_error,err,error,*999)
3278  END SELECT
3279  ELSE
3280  CALL flagerror("The backward Euler differential-algebraic equation solver is not associated.",err,error,*999)
3281  ENDIF
3283  improved_euler_dae_solver=>euler_dae_solver%IMPROVED_EULER_SOLVER
3284  IF(ASSOCIATED(improved_euler_dae_solver)) THEN
3285  SELECT CASE(solver_library_type)
3286  CASE(solver_cmiss_library)
3287  CALL flagerror("Not implemented.",err,error,*999)
3288  CASE(solver_petsc_library)
3289  CALL flagerror("Not implemented.",err,error,*999)
3290  CASE DEFAULT
3291  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3292  & " is invalid."
3293  CALL flagerror(local_error,err,error,*999)
3294  END SELECT
3295  ELSE
3296  CALL flagerror("The improved Euler differential-algebraic equation solver is not associated.",err,error,*999)
3297  ENDIF
3298  CASE DEFAULT
3299  local_error="The Euler differential-algebraic equations solver type of "// &
3300  & trim(numbertovstring(euler_dae_solver%EULER_TYPE,"*",err,error))//" is invalid."
3301  CALL flagerror(local_error,err,error,*999)
3302  END SELECT
3303  ELSE
3304  CALL flagerror("The Euler differential-algebraic equation solver is not associated.",err,error,*999)
3305  ENDIF
3307  crank_nicolson_dae_solver=>dae_solver%CRANK_NICOLSON_SOLVER
3308  IF(ASSOCIATED(crank_nicolson_dae_solver)) THEN
3309  SELECT CASE(solver_library_type)
3310  CASE(solver_cmiss_library)
3311  CALL flagerror("Not implemented.",err,error,*999)
3312  CASE(solver_petsc_library)
3313  CALL flagerror("Not implemented.",err,error,*999)
3314  CASE DEFAULT
3315  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3316  & " is invalid."
3317  CALL flagerror(local_error,err,error,*999)
3318  END SELECT
3319  ELSE
3320  CALL flagerror("The Crank-Nicolson differential-algebraic equation solver is not associated.",err,error,*999)
3321  ENDIF
3323  runge_kutta_dae_solver=>dae_solver%RUNGE_KUTTA_SOLVER
3324  IF(ASSOCIATED(runge_kutta_dae_solver)) THEN
3325  SELECT CASE(solver_library_type)
3326  CASE(solver_cmiss_library)
3327  CALL flagerror("Not implemented.",err,error,*999)
3328  CASE(solver_petsc_library)
3329  CALL flagerror("Not implemented.",err,error,*999)
3330  CASE DEFAULT
3331  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3332  & " is invalid."
3333  CALL flagerror(local_error,err,error,*999)
3334  END SELECT
3335  ELSE
3336  CALL flagerror("The Runge-Kutta differential-algebraic equation solver is not associated.",err,error,*999)
3337  ENDIF
3339  adams_moulton_dae_solver=>dae_solver%ADAMS_MOULTON_SOLVER
3340  IF(ASSOCIATED(adams_moulton_dae_solver)) THEN
3341  SELECT CASE(solver_library_type)
3342  CASE(solver_cmiss_library)
3343  CALL flagerror("Not implemented.",err,error,*999)
3344  CASE(solver_petsc_library)
3345  CALL flagerror("Not implemented.",err,error,*999)
3346  CASE DEFAULT
3347  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3348  & " is invalid."
3349  CALL flagerror(local_error,err,error,*999)
3350  END SELECT
3351  ELSE
3352  CALL flagerror("The Adams-Moulton differential-algebraic equation solver is not associated.",err,error,*999)
3353  ENDIF
3354  CASE(solver_dae_bdf)
3355  bdf_dae_solver=>dae_solver%BDF_SOLVER
3356  IF(ASSOCIATED(bdf_dae_solver)) THEN
3357  SELECT CASE(solver_library_type)
3358  CASE(solver_cmiss_library)
3359  CALL flagerror("Not implemented.",err,error,*999)
3360  CASE(solver_petsc_library)
3361  bdf_dae_solver%SOLVER_LIBRARY = solver_petsc_library
3362  CASE DEFAULT
3363  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3364  & " is invalid."
3365  CALL flagerror(local_error,err,error,*999)
3366  END SELECT
3367  ELSE
3368  CALL flagerror("The BDF differential-algebraic equation solver is not associated.",err,error,*999)
3369  ENDIF
3371  rush_larson_dae_solver=>dae_solver%RUSH_LARSON_SOLVER
3372  IF(ASSOCIATED(rush_larson_dae_solver)) THEN
3373  SELECT CASE(solver_library_type)
3374  CASE(solver_cmiss_library)
3375  CALL flagerror("Not implemented.",err,error,*999)
3376  CASE(solver_petsc_library)
3377  CALL flagerror("Not implemented.",err,error,*999)
3378  CASE DEFAULT
3379  local_error="The solver library type of "//trim(numbertovstring(solver_library_type,"*",err,error))// &
3380  & " is invalid."
3381  CALL flagerror(local_error,err,error,*999)
3382  END SELECT
3383  ELSE
3384  CALL flagerror("The Rush-Larson differential-algebraic equation solver is not associated.",err,error,*999)
3385  ENDIF
3386  CASE(solver_dae_external)
3387  CALL flagerror("Can not set the library type for an external differential-algebraic equation solver is not associated.", &
3388  & err,error,*999)
3389  CASE DEFAULT
3390  local_error="The differential-algebraic equations solver type of "// &
3391  & trim(numbertovstring(dae_solver%DAE_SOLVE_TYPE,"*",err,error))//" is invalid."
3392  CALL flagerror(local_error,err,error,*999)
3393  END SELECT
3394  ELSE
3395  CALL flagerror("DAE solver is not associated.",err,error,*999)
3396  ENDIF
3397 
3398  exits("SOLVER_DAE_LIBRARY_TYPE_SET")
3399  RETURN
3400 999 errorsexits("SOLVER_DAE_LIBRARY_TYPE_SET",err,error)
3401  RETURN 1
3402 
3403  END SUBROUTINE solver_dae_library_type_set
3404 
3405  !
3406  !================================================================================================================================
3407  !
3408 
3410  SUBROUTINE solver_dae_bdf_finalise(BDF_SOLVER,ERR,ERROR,*)
3412  !Argument variables
3413  TYPE(bdf_dae_solver_type), POINTER :: BDF_SOLVER
3414  INTEGER(INTG), INTENT(OUT) :: ERR
3415  TYPE(varying_string), INTENT(OUT) :: ERROR
3416  !Local Variables
3417 
3418  enters("SOLVER_DAE_BDF_FINALISE",err,error,*999)
3419 
3420  IF(ASSOCIATED(bdf_solver)) THEN
3421  DEALLOCATE(bdf_solver)
3422  ENDIF
3423 
3424  exits("SOLVER_DAE_BDF_FINALISE")
3425  RETURN
3426 999 errorsexits("SOLVER_DAE_BDF_FINALISE",err,error)
3427  RETURN 1
3428 
3429  END SUBROUTINE solver_dae_bdf_finalise
3430 
3431  !
3432  !================================================================================================================================
3433  !
3434 
3436  SUBROUTINE solver_dae_bdf_initialise(DAE_SOLVER,ERR,ERROR,*)
3438  !Argument variables
3439  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
3440  INTEGER(INTG), INTENT(OUT) :: ERR
3441  TYPE(varying_string), INTENT(OUT) :: ERROR
3442  !Local Variables
3443  INTEGER(INTG) :: DUMMY_ERR
3444  TYPE(varying_string) :: DUMMY_ERROR
3445 
3446  enters("SOLVER_DAE_BDF_INITIALISE",err,error,*998)
3447 
3448  IF(ASSOCIATED(dae_solver)) THEN
3449  IF(ASSOCIATED(dae_solver%BDF_SOLVER)) THEN
3450  CALL flagerror("BDF solver is already associated for this differential-algebraic equation solver.",err,error,*998)
3451  ELSE
3452  !Allocate the BDF solver
3453  ALLOCATE(dae_solver%BDF_SOLVER,stat=err)
3454  IF(err/=0) CALL flagerror("Could not allocate BDF solver.",err,error,*999)
3455  !Initialise
3456  dae_solver%BDF_SOLVER%DAE_SOLVER=>dae_solver
3457  dae_solver%BDF_SOLVER%SOLVER_LIBRARY=solver_petsc_library
3458  !Defaults
3459  ENDIF
3460  ELSE
3461  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*998)
3462  ENDIF
3463 
3464  exits("SOLVER_DAE_BDF_INITIALISE")
3465  RETURN
3466 999 CALL solver_dae_bdf_finalise(dae_solver%BDF_SOLVER,dummy_err,dummy_error,*998)
3467 998 errorsexits("SOLVER_DAE_BDF_INITIALISE",err,error)
3468  RETURN 1
3469 
3470  END SUBROUTINE solver_dae_bdf_initialise
3471  !
3472  !================================================================================================================================
3473  !
3474 
3476  SUBROUTINE solver_daecellmlpetsccontextfinalise(ctx,err,error,*)
3478  !Argument variables
3479  TYPE(cellmlpetsccontexttype), POINTER :: ctx
3480  INTEGER(INTG), INTENT(OUT) :: err
3481  TYPE(varying_string), INTENT(OUT) :: error
3482  !Local Variables
3483 
3484  enters("Solver_DAECellMLPETScContextFinalise",err,error,*999)
3485 
3486  IF(ASSOCIATED(ctx)) THEN
3487  IF(ASSOCIATED(ctx%rates)) DEALLOCATE(ctx%rates)
3488  IF(ALLOCATED(ctx%ratesIndices)) DEALLOCATE(ctx%ratesIndices)
3489  DEALLOCATE(ctx)
3490  ENDIF
3491 
3492  exits("Solver_DAECellMLPETScContextFinalise")
3493  RETURN
3494 999 errorsexits("Solver_DAECellMLPETScContextFinalise",err,error)
3495  RETURN 1
3496 
3498 
3499 
3500  !
3501  !================================================================================================================================
3502  !
3503 
3505  SUBROUTINE solver_daecellmlpetsccontextinitialise(ctx,err,error,*)
3507  !Argument variables
3508  TYPE(cellmlpetsccontexttype), INTENT(OUT), POINTER :: ctx
3509  INTEGER(INTG), INTENT(OUT) :: err
3510  TYPE(varying_string), INTENT(OUT) :: error
3511  !Local Variables
3512  INTEGER(INTG) :: dummyErr
3513  TYPE(varying_string) :: dummyError
3514 
3515  enters("Solver_DAECellMLPETScContextInitialise",err,error,*998)
3516 
3517  IF(ASSOCIATED(ctx)) THEN
3518  CALL flagerror("Context is already associated.",err,error,*998)
3519  ELSE
3520  !Allocate the CTX
3521  ALLOCATE(ctx,stat=err)
3522  IF(err/=0) CALL flagerror("Could not allocate context.",err,error,*999)
3523  !Initialise
3524  NULLIFY(ctx%solver)
3525  NULLIFY(ctx%cellml)
3526  NULLIFY(ctx%rates)
3527  ctx%dofIdx=0
3528  ENDIF
3529 
3530  exits("Solver_DAECellMLPETScContextInitialise")
3531  RETURN
3532 999 CALL solver_daecellmlpetsccontextfinalise(ctx,dummyerr,dummyerror,*998)
3533 998 errorsexits("Solver_DAECellMLPETScContextInitialise",err,error)
3534  RETURN 1
3535 
3537  !
3538  !================================================================================================================================
3539  !
3540 
3542  SUBROUTINE solver_daecellmlpetsccontextset(ctx,solver,cellml,dofIdx,err,error,*)
3544  !Argument variables
3545  TYPE(cellmlpetsccontexttype), INTENT(IN), POINTER :: ctx
3546  TYPE(solver_type), POINTER, INTENT(IN) :: solver
3547  TYPE(cellml_type), POINTER, INTENT(IN) :: cellml
3548  INTEGER(INTG), INTENT(IN) :: dofIdx
3549  INTEGER(INTG), INTENT(OUT) :: err
3550  TYPE(varying_string), INTENT(OUT) :: error
3551  !Local Variables
3552  INTEGER(INTG) :: arrayIdx,dummyErr
3553  TYPE(varying_string) :: dummyError
3554 
3555  enters("Solver_DAECellMLPETScContextSet",err,error,*998)
3556 
3557  IF(ASSOCIATED(ctx)) THEN
3558  IF(ASSOCIATED(solver)) THEN
3559  IF(ASSOCIATED(cellml)) THEN
3560  !Set
3561  ctx%solver=>solver
3562  ctx%cellml=>cellml
3563  ctx%dofIdx=dofidx
3564  ALLOCATE(ctx%rates(cellml%MAXIMUM_NUMBER_OF_STATE),stat=err)
3565  IF(err/=0) CALL flagerror("Could not allocate context rates.",err,error,*999)
3566  ALLOCATE(ctx%ratesIndices(cellml%MAXIMUM_NUMBER_OF_STATE),stat=err)
3567  IF(err/=0) CALL flagerror("Could not allocate context rates.",err,error,*999)
3568  ctx%ratesIndices=[(arrayidx,arrayidx=0,(cellml%MAXIMUM_NUMBER_OF_STATE-1))]
3569  ELSE
3570  CALL flagerror("CellML environment is not associated.",err,error,*999)
3571  ENDIF
3572  ELSE
3573  CALL flagerror("Solver is not associated.",err,error,*998)
3574  ENDIF
3575  ELSE
3576  CALL flagerror("ctx is not associated.",err,error,*998)
3577  ENDIF
3578 
3579  exits("Solver_DAECellMLPETScContextSet")
3580  RETURN
3581 999 CALL solver_daecellmlpetsccontextfinalise(ctx,dummyerr,dummyerror,*998)
3582 998 errorsexits("Solver_DAECellMLPETScContextSet",err,error)
3583  RETURN 1
3584 
3585  END SUBROUTINE solver_daecellmlpetsccontextset
3586 
3587  !
3588  !================================================================================================================================
3589  !
3591  SUBROUTINE solver_dae_bdf_integrate(BDF_SOLVER,CELLML,N,START_TIME,END_TIME,TIME_INCREMENT, &
3592  & only_one_model_index,models_data,max_number_states,state_data,max_number_parameters,parameters_data, &
3593  & max_number_intermediates,intermediate_data,err,error,*)
3595  !Argument variables
3596  TYPE(bdf_dae_solver_type), POINTER :: BDF_SOLVER
3597  TYPE(cellml_type), POINTER :: CELLML
3598  INTEGER(INTG), INTENT(IN) :: N
3599  REAL(DP), INTENT(IN) :: START_TIME
3600  REAL(DP), INTENT(IN) :: END_TIME
3601  REAL(DP), INTENT(INOUT) :: TIME_INCREMENT
3602  INTEGER(INTG), INTENT(IN) :: ONLY_ONE_MODEL_INDEX
3603  INTEGER(INTG), POINTER, INTENT(IN) :: MODELS_DATA(:)
3604  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_STATES
3605  REAL(DP), POINTER, INTENT (INOUT) :: STATE_DATA(:)
3606  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_PARAMETERS
3607  REAL(DP), POINTER, INTENT(INOUT) :: PARAMETERS_DATA(:)
3608  INTEGER(INTG), INTENT(IN) :: MAX_NUMBER_INTERMEDIATES
3609  REAL(DP), POINTER, INTENT(INOUT) :: INTERMEDIATE_DATA(:)
3610  INTEGER(INTG), INTENT(OUT) :: ERR
3611  TYPE(varying_string), INTENT(OUT) :: ERROR
3612  !Local Variables
3613  TYPE(petsctstype) :: ts
3614  REAL(DP) :: FINALSOLVEDTIME,TIMESTEP
3615  TYPE(petscvectype) :: PETSC_CURRENT_STATES
3616  TYPE(cellmlpetsccontexttype), POINTER :: CTX
3617  INTEGER(INTG) :: dof_idx,DOF_ORDER_TYPE,model_idx, NUMBER_STATES,STATE_END_DOF,state_idx,STATE_START_DOF,array_idx
3618  REAL(DP), ALLOCATABLE :: STATES_TEMP(:),RATES_TEMP(:)
3619  INTEGER(INTG), ALLOCATABLE :: ARRAY_INDICES(:)
3620  TYPE(cellml_model_type), POINTER :: MODEL
3621  TYPE(varying_string) :: LOCAL_ERROR
3622  TYPE(petscvectype) :: PETSC_RATES
3623  EXTERNAL :: problem_solverdaecellmlrhspetsc
3624 
3625 
3626  enters("SOLVER_DAE_BFD_INTEGRATE",err,error,*999)
3627 
3628  NULLIFY(ctx)
3629  timestep=end_time-start_time
3630  IF(ASSOCIATED(bdf_solver)) THEN
3631  IF(ASSOCIATED(cellml)) THEN
3632  IF(ASSOCIATED(cellml%MODELS_FIELD)) THEN
3633  SELECT CASE(bdf_solver%SOLVER_LIBRARY)
3634  CASE(solver_petsc_library)
3635  CALL field_dof_order_type_get(cellml%MODELS_FIELD%MODELS_FIELD, &
3636  & field_u_variable_type,dof_order_type,err,error,*999)
3637  IF(dof_order_type==field_separated_component_dof_order) THEN
3638 
3639  ELSE !dof component order is contiguous
3640  IF(only_one_model_index==cellml_models_field_not_constant) THEN
3641 
3642  ELSE !only one model
3643  model=>cellml%MODELS(only_one_model_index)%PTR
3644  IF(ASSOCIATED(model)) THEN
3645  !determine no. of states in model and allocate necessary arrays
3646  number_states = model%NUMBER_OF_STATE
3647  ALLOCATE(states_temp(0:number_states-1),stat=err)
3648  ALLOCATE(rates_temp(0:number_states-1),stat=err)
3649  ALLOCATE(array_indices(0:number_states-1),stat=err)
3650  array_indices = (/(array_idx,array_idx=0,(number_states-1))/)
3651 
3652 
3653  !initialize context for petsc solving.
3654  CALL solver_daecellmlpetsccontextinitialise(ctx,err,error,*999)
3655  DO dof_idx=1,n
3656  model_idx = models_data(dof_idx)
3657  IF(model_idx>0) THEN !if model is assigned to dof
3658  !access the state field data
3659  state_start_dof=(dof_idx-1)*max_number_states+1
3660  state_end_dof=state_start_dof+number_states-1
3661  DO state_idx=1,number_states
3662  states_temp(state_idx-1) = state_data(state_start_dof+state_idx-1)
3663  ENDDO
3664 
3665  !create PETSC states vector to initialize solver
3666  CALL petsc_veccreateseq(petsc_comm_self, &
3667  & number_states,petsc_current_states,err,error,*999)
3668  !CALL Petsc_VecSetSizes(PETSC_CURRENT_STATES, &
3669  ! & PETSC_DECIDE,(NUMBER_STATES),ERR,ERROR,*999)
3670  !CALL Petsc_VecSetFromOptions(PETSC_CURRENT_STATES,ERR,ERROR,*999)
3671 
3672  !create PETSC rates vector to return values from evaluating rhs routine
3673  CALL petsc_veccreateseq(petsc_comm_self, &
3674  & number_states,petsc_rates,err,error,*999)
3675  !CALL Petsc_VecSetSizes(PETSC_RATES, &
3676  ! & PETSC_DECIDE,(NUMBER_STATES),ERR,ERROR,*999)
3677  !CALL Petsc_VecSetFromOptions(PETSC_RATES,ERR,ERROR,*999)
3678 
3679  !Set up PETSC TS context for sundials BDF solver
3680  CALL petsc_tscreate(petsc_comm_self,ts,err,error,*999)
3681  CALL petsc_tssetproblemtype(ts,petsc_ts_nonlinear,err,error,*999)
3682  CALL petsc_tssettype(ts,petsc_ts_sundials,err,error,*999)
3683  CALL petsc_tssundialssettype(ts,petsc_sundials_bdf,err,error,*999)
3684  CALL petsc_tssundialssettolerance(ts,0.0000001_dp, &
3685  & 0.0000001_dp,err,error,*999)
3686  !set the initial solution to the current state
3687  CALL petsc_vecsetvalues(petsc_current_states,(number_states), &
3688  & array_indices,states_temp, &
3689  & petsc_insert_values,err,error,*999)
3690  CALL petsc_vecassemblybegin(petsc_current_states,err,error,*999)
3691  CALL petsc_vecassemblyend(petsc_current_states,err,error,*999)
3692  CALL petsc_tssetsolution(ts,petsc_current_states,err,error,*999)
3693 
3694  !set up the time data
3695  CALL petsc_tssetinitialtimestep(ts,start_time,time_increment,err,error,*999)
3696  CALL petsc_tssetduration(ts,5000,end_time,err,error,*999)
3697  CALL petsc_tssetexactfinaltime(ts,.true.,err,error,*999)
3698 
3699  IF(diagnostics1) THEN
3700  CALL write_string_value(diagnostic_output_type," DAE START TIME = ",start_time,err,error,*999)
3701  CALL write_string_value(diagnostic_output_type," DAE END TIME = ",end_time,err,error,*999)
3702  ENDIF
3703 
3704  !set rhs function and pass through the cellml model context
3705  CALL solver_daecellmlpetsccontextset(ctx,bdf_solver%DAE_SOLVER%SOLVER,cellml,dof_idx,err,error,*999)
3706  CALL petsc_tssetrhsfunction(ts,petsc_rates,problem_solverdaecellmlrhspetsc,ctx,err,error,*999)
3707 
3708  CALL petsc_tssolve(ts,petsc_current_states,finalsolvedtime,err,error,*999)
3709  IF(diagnostics1) THEN
3710  CALL write_string_value(diagnostic_output_type," FINAL SOLVED TIME = ", &
3711  & finalsolvedtime,err,error,*999)
3712  ENDIF
3713 
3714 
3715  !update the states to new integrated values
3716  CALL petsc_vecassemblybegin(petsc_current_states,err,error,*999)
3717  CALL petsc_vecassemblyend(petsc_current_states,err,error,*999)
3718  CALL petsc_vecgetvalues(petsc_current_states, &
3719  & number_states, array_indices, &
3720  & states_temp, &
3721  & err,error,*999)
3722 
3723  DO state_idx=1,number_states
3724  state_data(state_start_dof+state_idx-1)= &
3725  & states_temp(state_idx-1)
3726  ENDDO
3727  CALL petsc_tsfinalise(ts,err,error,*999)
3728  ENDIF !model_idx
3729  CALL petsc_vecdestroy(petsc_current_states,err,error,*999)
3730  CALL petsc_vecdestroy(petsc_rates,err,error,*999)
3731  ENDDO !dof_idx
3732 
3733  ELSE
3734  CALL flagerror("Cellml model is not associated.",err,error,*999)
3735  ENDIF
3736  ENDIF
3737  ENDIF !dof continguous
3738  CASE DEFAULT
3739  local_error="The BDF solver library type of "// &
3740  & trim(numbertovstring(bdf_solver%SOLVER_LIBRARY,"*",err,error))//" is not implemented."
3741  CALL flagerror(local_error,err,error,*999)
3742  END SELECT
3743  ELSE
3744  CALL flagerror("CELLML models field is not associated.",err,error,*999)
3745  ENDIF
3746  ELSE
3747  CALL flagerror("CELLML environment is not associated.",err,error,*999)
3748  ENDIF
3749  ELSE
3750  CALL flagerror("BDF solver is not associated.",err,error,*999)
3751  ENDIF
3752 
3753  exits("SOLVER_DAE_BDF_INTEGRATE")
3754  RETURN
3755 999 errorsexits("SOLVER_DAE_BDF_INTEGRATE",err,error)
3756  RETURN 1
3757 
3758  END SUBROUTINE solver_dae_bdf_integrate
3759  !
3760  !================================================================================================================================
3761  !
3762 
3764  SUBROUTINE solver_dae_bdf_solve(BDF_SOLVER,ERR,ERROR,*)
3766  !Argument variables
3767  TYPE(bdf_dae_solver_type), POINTER :: BDF_SOLVER
3768  INTEGER(INTG), INTENT(OUT) :: ERR
3769  TYPE(varying_string), INTENT(OUT) :: ERROR
3770  !Local Variables
3771  INTEGER(INTG) :: cellml_idx
3772  INTEGER(INTG), POINTER :: MODELS_DATA(:)
3773  REAL(DP), POINTER :: INTERMEDIATE_DATA(:),PARAMETERS_DATA(:),STATE_DATA(:)
3774  TYPE(cellml_type), POINTER :: CELLML_ENVIRONMENT
3775  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
3776  TYPE(cellml_models_field_type), POINTER :: CELLML_MODELS_FIELD
3777  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
3778  TYPE(field_variable_type), POINTER :: MODELS_VARIABLE
3779  TYPE(field_type), POINTER :: MODELS_FIELD,STATE_FIELD,PARAMETERS_FIELD,INTERMEDIATE_FIELD
3780  TYPE(solver_type), POINTER :: SOLVER
3781  TYPE(varying_string) :: LOCAL_ERROR
3782 
3783  enters("SOLVER_DAE_BDF_SOLVE",err,error,*999)
3784 
3785  NULLIFY(models_data)
3786  NULLIFY(intermediate_data)
3787  NULLIFY(parameters_data)
3788  NULLIFY(state_data)
3789  NULLIFY(models_variable)
3790 
3791  IF(ASSOCIATED(bdf_solver)) THEN
3792  dae_solver=>bdf_solver%DAE_SOLVER
3793  IF(ASSOCIATED(dae_solver)) THEN
3794  solver=>dae_solver%SOLVER
3795  IF(ASSOCIATED(solver)) THEN
3796  cellml_equations=>solver%CELLML_EQUATIONS
3797  IF(ASSOCIATED(cellml_equations)) THEN
3798  DO cellml_idx=1,cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS
3799  cellml_environment=>cellml_equations%CELLML_ENVIRONMENTS(cellml_idx)%PTR
3800  IF(ASSOCIATED(cellml_environment)) THEN
3801  cellml_models_field=>cellml_environment%MODELS_FIELD
3802  IF(ASSOCIATED(cellml_models_field)) THEN
3803  models_field=>cellml_models_field%MODELS_FIELD
3804  IF(ASSOCIATED(models_field)) THEN
3805 
3806 !!TODO: Maybe move this getting of fields earlier up the DAE solver chain? For now keep here.
3807 
3808  !Make sure CellML fields have been updated to the current value of any mapped fields
3809  CALL cellml_field_to_cellml_update(cellml_environment,err,error,*999)
3810 
3811  CALL field_variable_get(models_field,field_u_variable_type,models_variable,err,error,*999)
3812  CALL field_parameter_set_data_get(models_field,field_u_variable_type,field_values_set_type, &
3813  & models_data,err,error,*999)
3814 
3815  !Get the state information if this environment has any.
3816  IF(ASSOCIATED(cellml_environment%STATE_FIELD)) THEN
3817  state_field=>cellml_environment%STATE_FIELD%STATE_FIELD
3818  IF(ASSOCIATED(state_field)) THEN
3819  CALL field_parameter_set_data_get(state_field,field_u_variable_type,field_values_set_type, &
3820  & state_data,err,error,*999)
3821  ENDIF
3822  ENDIF
3823 
3824  !Get the parameters information if this environment has any.
3825  IF(ASSOCIATED(cellml_environment%PARAMETERS_FIELD)) THEN
3826  parameters_field=>cellml_environment%PARAMETERS_FIELD%PARAMETERS_FIELD
3827  IF(ASSOCIATED(parameters_field)) THEN
3828  CALL field_parameter_set_data_get(parameters_field,field_u_variable_type,field_values_set_type, &
3829  & parameters_data,err,error,*999)
3830  ENDIF
3831  ENDIF
3832 
3833  !Get the intermediate information if this environment has any.
3834  IF(ASSOCIATED(cellml_environment%INTERMEDIATE_FIELD)) THEN
3835  intermediate_field=>cellml_environment%INTERMEDIATE_FIELD%INTERMEDIATE_FIELD
3836  IF(ASSOCIATED(intermediate_field)) THEN
3837  CALL field_parameter_set_data_get(intermediate_field,field_u_variable_type,field_values_set_type, &
3838  & intermediate_data,err,error,*999)
3839  ENDIF
3840  ENDIF
3841 
3842  !Integrate these CellML equations
3843 
3844  CALL solver_dae_bdf_integrate(bdf_solver,cellml_environment,models_variable% &
3845  & total_number_of_dofs,dae_solver%START_TIME,dae_solver%END_TIME,dae_solver%INITIAL_STEP, &
3846  & cellml_environment%MODELS_FIELD%ONLY_ONE_MODEL_INDEX,models_data,cellml_environment% &
3847  & maximum_number_of_state,state_data,cellml_environment%MAXIMUM_NUMBER_OF_PARAMETERS, &
3848  & parameters_data,cellml_environment%MAXIMUM_NUMBER_OF_INTERMEDIATE,intermediate_data,err,error,*999)
3849 
3850  !Restore field data
3851  CALL field_parameter_set_data_restore(models_field,field_u_variable_type,field_values_set_type, &
3852  & models_data,err,error,*999)
3853  IF(ASSOCIATED(state_field)) CALL field_parameter_set_data_restore(state_field,field_u_variable_type, &
3854  & field_values_set_type,state_data,err,error,*999)
3855  IF(ASSOCIATED(parameters_field)) CALL field_parameter_set_data_restore(parameters_field, &
3856  & field_u_variable_type,field_values_set_type,parameters_data,err,error,*999)
3857  IF(ASSOCIATED(intermediate_field)) CALL field_parameter_set_data_restore(intermediate_field, &
3858  & field_u_variable_type,field_values_set_type,intermediate_data,err,error,*999)
3859 
3860  !Make sure fields have been updated to the current value of any mapped CellML fields
3861  CALL cellml_cellml_to_field_update(cellml_environment,err,error,*999)
3862 
3863  ELSE
3864  local_error="The CellML models field is not associated for CellML index "// &
3865  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
3866  CALL flagerror(local_error,err,error,*999)
3867  ENDIF
3868  ELSE
3869  local_error="The CellML models field is not associated for CellML index "// &
3870  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
3871  CALL flagerror(local_error,err,error,*999)
3872  ENDIF
3873  ELSE
3874  local_error="The CellML enviroment is not associated for for CellML index "// &
3875  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
3876  CALL flagerror(local_error,err,error,*999)
3877  ENDIF
3878  ENDDO !cellml_idx
3879  ELSE
3880  CALL flagerror("Solver solver equations is not associated.",err,error,*999)
3881  ENDIF
3882  ELSE
3883  CALL flagerror("Solver is not associated.",err,error,*999)
3884  ENDIF
3885  ELSE
3886  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*999)
3887  ENDIF
3888  ELSE
3889  CALL flagerror("BDF differential-algebraic equation solver is not associated.",err,error,*999)
3890  ENDIF
3891 
3892  exits("SOLVER_DAE_BDF_SOLVE")
3893  RETURN
3894 999 errorsexits("SOLVER_DAE_BDF_SOLVE",err,error)
3895  RETURN 1
3896 
3897  END SUBROUTINE solver_dae_bdf_solve
3898 
3899  !
3900  !================================================================================================================================
3901  !
3902 
3904  SUBROUTINE solver_dae_crank_nicolson_finalise(CRANK_NICOLSON_SOLVER,ERR,ERROR,*)
3906  !Argument variables
3907  TYPE(crank_nicolson_dae_solver_type), POINTER :: CRANK_NICOLSON_SOLVER
3908  INTEGER(INTG), INTENT(OUT) :: ERR
3909  TYPE(varying_string), INTENT(OUT) :: ERROR
3910  !Local Variables
3911 
3912  enters("SOLVER_DAE_CRANK_NICOLSON_FINALISE",err,error,*999)
3913 
3914  IF(ASSOCIATED(crank_nicolson_solver)) THEN
3915  DEALLOCATE(crank_nicolson_solver)
3916  ENDIF
3917 
3918  exits("SOLVER_DAE_CRANK_NICOLSON_FINALISE")
3919  RETURN
3920 999 errorsexits("SOLVER_DAE_CRANK_NICOLSON_FINALISE",err,error)
3921  RETURN 1
3922 
3923  END SUBROUTINE solver_dae_crank_nicolson_finalise
3924 
3925  !
3926  !================================================================================================================================
3927  !
3928 
3930  SUBROUTINE solver_dae_crank_nicolson_initialise(DAE_SOLVER,ERR,ERROR,*)
3932  !Argument variables
3933  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
3934  INTEGER(INTG), INTENT(OUT) :: ERR
3935  TYPE(varying_string), INTENT(OUT) :: ERROR
3936  !Local Variables
3937  INTEGER(INTG) :: DUMMY_ERR
3938  TYPE(varying_string) :: DUMMY_ERROR
3939 
3940  enters("SOLVER_DAE_CRANK_NICOLSON_INITIALISE",err,error,*998)
3941 
3942  IF(ASSOCIATED(dae_solver)) THEN
3943  IF(ASSOCIATED(dae_solver%CRANK_NICOLSON_SOLVER)) THEN
3944  CALL flagerror("Crank-Nicolson solver is already associated for this differential-algebraic equation solver.", &
3945  & err,error,*998)
3946  ELSE
3947  !Allocate the Crank-Nicholson solver
3948  ALLOCATE(dae_solver%CRANK_NICOLSON_SOLVER,stat=err)
3949  IF(err/=0) CALL flagerror("Could not allocate Crank-Nicolson solver.",err,error,*999)
3950  !Initialise
3951  dae_solver%CRANK_NICOLSON_SOLVER%DAE_SOLVER=>dae_solver
3952  dae_solver%CRANK_NICOLSON_SOLVER%SOLVER_LIBRARY=0
3953  !Defaults
3954  ENDIF
3955  ELSE
3956  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*998)
3957  ENDIF
3958 
3959  exits("SOLVER_DAE_CRANK_NICOLSON_INITIALISE")
3960  RETURN
3961 999 CALL solver_dae_crank_nicolson_finalise(dae_solver%CRANK_NICOLSON_SOLVER,dummy_err,dummy_error,*998)
3962 998 errorsexits("SOLVER_DAE_CRANK_NICOLSON_INITIALISE",err,error)
3963  RETURN 1
3964 
3966 
3967  !
3968  !================================================================================================================================
3969  !
3970 
3972  SUBROUTINE solver_dae_crank_nicolson_solve(CRANK_NICOLSON_SOLVER,ERR,ERROR,*)
3974  !Argument variables
3975  TYPE(crank_nicolson_dae_solver_type), POINTER :: CRANK_NICOLSON_SOLVER
3976  INTEGER(INTG), INTENT(OUT) :: ERR
3977  TYPE(varying_string), INTENT(OUT) :: ERROR
3978  !Local Variables
3979 
3980  enters("SOLVER_DAE_CRANK_NICOLSON_SOLVE",err,error,*999)
3981 
3982  IF(ASSOCIATED(crank_nicolson_solver)) THEN
3983  CALL flagerror("Not implemented.",err,error,*999)
3984  ELSE
3985  CALL flagerror("Crank-Nicolson differential-algebraic equation solver is not associated.",err,error,*999)
3986  ENDIF
3987 
3988  exits("SOLVER_DAE_CRANK_NICOLSON_SOLVE")
3989  RETURN
3990 999 errorsexits("SOLVER_DAE_CRANK_NICOLSON_SOLVE",err,error)
3991  RETURN 1
3992 
3993  END SUBROUTINE solver_dae_crank_nicolson_solve
3994 
3995  !
3996  !================================================================================================================================
3997  !
3998 
4000  SUBROUTINE solver_dae_external_finalise(EXTERNAL_SOLVER,ERR,ERROR,*)
4002  !Argument variables
4003  TYPE(external_dae_solver_type), POINTER :: EXTERNAL_SOLVER
4004  INTEGER(INTG), INTENT(OUT) :: ERR
4005  TYPE(varying_string), INTENT(OUT) :: ERROR
4006  !Local Variables
4007 
4008  enters("SOLVER_DAE_EXTERNAL_FINALISE",err,error,*999)
4009 
4010  IF(ASSOCIATED(external_solver)) THEN
4011  DEALLOCATE(external_solver)
4012  ENDIF
4013 
4014  exits("SOLVER_DAE_CRANK_NICOLSON_FINALISE")
4015  RETURN
4016 999 errorsexits("SOLVER_DAE_CRANK_NICOLSON_FINALISE",err,error)
4017  RETURN 1
4018 
4019  END SUBROUTINE solver_dae_external_finalise
4020 
4021  !
4022  !================================================================================================================================
4023  !
4024 
4026  SUBROUTINE solver_dae_external_initialise(DAE_SOLVER,ERR,ERROR,*)
4028  !Argument variables
4029  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
4030  INTEGER(INTG), INTENT(OUT) :: ERR
4031  TYPE(varying_string), INTENT(OUT) :: ERROR
4032  !Local Variables
4033  INTEGER(INTG) :: DUMMY_ERR
4034  TYPE(varying_string) :: DUMMY_ERROR
4035 
4036  enters("SOLVER_DAE_EXTERNAL_INITIALISE",err,error,*998)
4037 
4038  IF(ASSOCIATED(dae_solver)) THEN
4039  IF(ASSOCIATED(dae_solver%EXTERNAL_SOLVER)) THEN
4040  CALL flagerror("External solver is already associated for this differential-algebraic equation solver.", &
4041  & err,error,*998)
4042  ELSE
4043  !Allocate the external solver
4044  ALLOCATE(dae_solver%EXTERNAL_SOLVER,stat=err)
4045  IF(err/=0) CALL flagerror("Could not allocate external solver.",err,error,*999)
4046  !Initialise
4047  dae_solver%EXTERNAL_SOLVER%DAE_SOLVER=>dae_solver
4048  !Defaults
4049  ENDIF
4050  ELSE
4051  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*998)
4052  ENDIF
4053 
4054  exits("SOLVER_DAE_EXTERNAL_INITIALISE")
4055  RETURN
4056 999 CALL solver_dae_external_finalise(dae_solver%EXTERNAL_SOLVER,dummy_err,dummy_error,*998)
4057 998 errorsexits("SOLVER_DAE_EXTERNAL_INITIALISE",err,error)
4058  RETURN 1
4059 
4060  END SUBROUTINE solver_dae_external_initialise
4061 
4062  !
4063  !================================================================================================================================
4064  !
4065 
4067  SUBROUTINE solver_dae_external_solve(EXTERNAL_SOLVER,ERR,ERROR,*)
4069  !Argument variables
4070  TYPE(external_dae_solver_type), POINTER :: EXTERNAL_SOLVER
4071  INTEGER(INTG), INTENT(OUT) :: ERR
4072  TYPE(varying_string), INTENT(OUT) :: ERROR
4073  !Local Variables
4074  INTEGER(INTG) :: cellml_idx
4075  INTEGER(INTG), POINTER :: MODELS_DATA(:)
4076  REAL(DP), POINTER :: INTERMEDIATE_DATA(:),PARAMETERS_DATA(:),STATE_DATA(:)
4077  TYPE(cellml_type), POINTER :: CELLML_ENVIRONMENT
4078  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
4079  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
4080  TYPE(field_variable_type), POINTER :: MODELS_VARIABLE
4081  TYPE(field_type), POINTER :: MODELS_FIELD,STATE_FIELD,PARAMETERS_FIELD,INTERMEDIATE_FIELD
4082  TYPE(solver_type), POINTER :: SOLVER
4083  TYPE(varying_string) :: LOCAL_ERROR
4084 
4085  enters("SOLVER_DAE_EXTERNAL_SOLVE",err,error,*999)
4086 
4087  NULLIFY(models_data)
4088  NULLIFY(intermediate_data)
4089  NULLIFY(parameters_data)
4090  NULLIFY(state_data)
4091 
4092  IF(ASSOCIATED(external_solver)) THEN
4093  dae_solver=>external_solver%DAE_SOLVER
4094  IF(ASSOCIATED(dae_solver)) THEN
4095  solver=>dae_solver%SOLVER
4096  IF(ASSOCIATED(solver)) THEN
4097  cellml_equations=>solver%CELLML_EQUATIONS
4098  IF(ASSOCIATED(cellml_equations)) THEN
4099  DO cellml_idx=1,cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS
4100  cellml_environment=>cellml_equations%CELLML_ENVIRONMENTS(cellml_idx)%PTR
4101  IF(ASSOCIATED(cellml_environment)) THEN
4102  IF(ASSOCIATED(cellml_environment%MODELS_FIELD)) THEN
4103  models_field=>cellml_environment%MODELS_FIELD%MODELS_FIELD
4104  IF(ASSOCIATED(models_field)) THEN
4105 
4106  !Make sure CellML fields have been updated to the current value of any mapped fields
4107  CALL cellml_field_to_cellml_update(cellml_environment,err,error,*999)
4108 
4109  NULLIFY(models_variable)
4110  CALL field_variable_get(models_field,field_u_variable_type,models_variable,err,error,*999)
4111  CALL field_parameter_set_data_get(models_field,field_u_variable_type,field_values_set_type, &
4112  & models_data,err,error,*999)
4113 
4114  !Get the state information if this environment has any.
4115  IF(ASSOCIATED(cellml_environment%STATE_FIELD)) THEN
4116  state_field=>cellml_environment%STATE_FIELD%STATE_FIELD
4117  IF(ASSOCIATED(state_field)) THEN
4118  CALL field_parameter_set_data_get(state_field,field_u_variable_type,field_values_set_type, &
4119  & state_data,err,error,*999)
4120  ELSE
4121  NULLIFY(state_data)
4122  ENDIF
4123  ELSE
4124  NULLIFY(state_data)
4125  ENDIF
4126 
4127  !Get the parameters information if this environment has any.
4128  IF(ASSOCIATED(cellml_environment%PARAMETERS_FIELD)) THEN
4129  parameters_field=>cellml_environment%PARAMETERS_FIELD%PARAMETERS_FIELD
4130  IF(ASSOCIATED(parameters_field)) THEN
4131  CALL field_parameter_set_data_get(parameters_field,field_u_variable_type,field_values_set_type, &
4132  & parameters_data,err,error,*999)
4133  ELSE
4134  NULLIFY(parameters_data)
4135  ENDIF
4136  ELSE
4137  NULLIFY(parameters_data)
4138  ENDIF
4139 
4140  !Get the intermediate information if this environment has any.
4141  IF(ASSOCIATED(cellml_environment%INTERMEDIATE_FIELD)) THEN
4142  intermediate_field=>cellml_environment%INTERMEDIATE_FIELD%INTERMEDIATE_FIELD
4143  IF(ASSOCIATED(intermediate_field)) THEN
4144  CALL field_parameter_set_data_get(intermediate_field,field_u_variable_type,field_values_set_type, &
4145  & intermediate_data,err,error,*999)
4146  ELSE
4147  NULLIFY(intermediate_data)
4148  ENDIF
4149  ELSE
4150  NULLIFY(intermediate_data)
4151  ENDIF
4152 
4153  !Call the external solver to integrate these CellML equations
4154  CALL solver_dae_external_integrate(models_variable%TOTAL_NUMBER_OF_DOFS,dae_solver%START_TIME, &
4155  & dae_solver%END_TIME,dae_solver%INITIAL_STEP,cellml_environment%MODELS_FIELD% &
4156  & only_one_model_index,models_data,cellml_environment%MAXIMUM_NUMBER_OF_STATE,state_data, &
4157  & cellml_environment%MAXIMUM_NUMBER_OF_PARAMETERS,parameters_data,cellml_environment% &
4158  & maximum_number_of_intermediate,intermediate_data,err)
4159  IF(err/=0) THEN
4160  error="Error from external solver integrate."
4161  GOTO 999
4162  ENDIF
4163 
4164  !Restore field data
4165  CALL field_parameter_set_data_restore(models_field,field_u_variable_type,field_values_set_type, &
4166  & models_data,err,error,*999)
4167  IF(ASSOCIATED(state_field)) CALL field_parameter_set_data_restore(state_field,field_u_variable_type, &
4168  & field_values_set_type,state_data,err,error,*999)
4169  IF(ASSOCIATED(parameters_field)) CALL field_parameter_set_data_restore(parameters_field, &
4170  & field_u_variable_type,field_values_set_type,parameters_data,err,error,*999)
4171  IF(ASSOCIATED(intermediate_field)) CALL field_parameter_set_data_restore(intermediate_field, &
4172  & field_u_variable_type,field_values_set_type,intermediate_data,err,error,*999)
4173 
4174  !Make sure fields have been updated to the current value of any mapped CellML fields
4175  CALL cellml_cellml_to_field_update(cellml_environment,err,error,*999)
4176 
4177  ELSE
4178  local_error="The CellML models field is not associated for CellML index "// &
4179  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
4180  CALL flagerror(local_error,err,error,*999)
4181  ENDIF
4182  ELSE
4183  local_error="The CellML models field is not associated for CellML index "// &
4184  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
4185  CALL flagerror(local_error,err,error,*999)
4186  ENDIF
4187  ELSE
4188  local_error="The CellML enviroment is not associated for for CellML index "// &
4189  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
4190  CALL flagerror(local_error,err,error,*999)
4191  ENDIF
4192  ENDDO !cellml_idx
4193  ELSE
4194  CALL flagerror("Solver solver equations is not associated.",err,error,*999)
4195  ENDIF
4196  ELSE
4197  CALL flagerror("Solver is not associated.",err,error,*999)
4198  ENDIF
4199  ELSE
4200  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*999)
4201  ENDIF
4202  ELSE
4203  CALL flagerror("External Euler differential-algebraic equation solver is not associated.",err,error,*999)
4204  ENDIF
4205 
4206  exits("SOLVER_DAE_EXTERNAL_SOLVE")
4207  RETURN
4208 999 errorsexits("SOLVER_DAE_EXTERNAL_SOLVE",err,error)
4209  RETURN 1
4210 
4211  END SUBROUTINE solver_dae_external_solve
4212 
4213  !
4214  !================================================================================================================================
4215  !
4216 
4218  SUBROUTINE solver_daecellmlrhsevaluate(model,time,stateStartIdx,stateDataOffset,stateData,parameterStartIdx,parameterDataOffset, &
4219  & parameterdata,intermediatestartidx,intermediatedataoffset,intermediatedata,ratestartidx,ratedataoffset,ratedata,err,error,*)
4221  !Argument variables
4222  TYPE(cellml_model_type), POINTER :: model
4223  REAL(DP), INTENT(IN) :: time
4224  INTEGER(INTG), INTENT(IN) :: stateStartIdx
4225  INTEGER(INTG), INTENT(IN) :: stateDataOffset
4226  REAL(DP), POINTER :: stateData(:)
4227  INTEGER(INTG), INTENT(IN) :: parameterStartIdx
4228  INTEGER(INTG), INTENT(IN) :: parameterDataOffset
4229  REAL(DP), POINTER :: parameterData(:)
4230  INTEGER(INTG), INTENT(IN) :: intermediateStartIdx
4231  INTEGER(INTG), INTENT(IN) :: intermediateDataOffset
4232  REAL(DP), POINTER :: intermediateData(:)
4233  INTEGER(INTG), INTENT(IN) :: rateStartIdx
4234  INTEGER(INTG), INTENT(IN) :: rateDataOffset
4235  REAL(DP), POINTER :: rateData(:)
4236  INTEGER(INTG), INTENT(OUT) :: err
4237  TYPE(varying_string), INTENT(OUT) :: error
4238  !Local Variables
4239  INTEGER(INTG) :: intermediateIdx,intermediateEndDOF,intermediateStartDOF,numberOfIntermediates,numberOfParameters, &
4240  & numberOfStates,parameterIdx,parameterEndDOF,parameterStartDOF,rateIdx,rateEndDOF,rateStartDOF,stateIdx,stateEndDOF, &
4241  & stateStartDOF
4242  REAL(DP) :: intermediates(max(1,intermediatedataoffset)),parameters(max(1,parameterdataoffset)),rates(max(1,ratedataoffset)), &
4243  & states(MAX(1,stateDataOffset))
4244 
4245  enters("Solver_DAECellMLRHSEvaluate",err,error,*999)
4246 
4247 #ifdef WITH_CELLML
4248 
4249  IF(ASSOCIATED(model)) THEN
4250  numberofstates=model%NUMBER_OF_STATE
4251  numberofintermediates=model%NUMBER_OF_INTERMEDIATE
4252  numberofparameters=model%NUMBER_OF_PARAMETERS
4253  IF(numberofstates>0) THEN
4254  IF(.NOT.ASSOCIATED(statedata)) CALL flagerror("State data is not associated.",err,error,*999)
4255  IF(.NOT.ASSOCIATED(ratedata)) CALL flagerror("Rate data is not associated.",err,error,*999)
4256  ENDIF
4257  IF(numberofparameters>0) THEN
4258  IF(.NOT.ASSOCIATED(parameterdata)) CALL flagerror("Parameter data is not associated.",err,error,*999)
4259  ENDIF
4260  IF(numberofintermediates>0) THEN
4261  IF(.NOT.ASSOCIATED(intermediatedata)) CALL flagerror("Intermediate data is not associated.",err,error,*999)
4262  ENDIF
4263  IF(statedataoffset>1.OR.numberofstates==0) THEN
4264  !State data is not contiguous or there are no states
4265 
4266  !Copy state data to temporary array
4267  DO stateidx=1,numberofstates
4268  states(stateidx)=statedata((statestartidx-1)*statedataoffset+stateidx)
4269  ENDDO !stateIdx
4270 
4271  IF(parameterdataoffset>1.OR.numberofparameters==0) THEN
4272  !Parameter data is not contiguous or there are no parameters
4273 
4274  !Copy parameter data to temporary array
4275  DO parameteridx=1,numberofparameters
4276  parameters(parameteridx)=parameterdata((parameterstartidx-1)*parameterdataoffset+parameteridx)
4277  ENDDO !parameterIdx
4278 
4279  IF(intermediatedataoffset>1.OR.numberofintermediates==0) THEN
4280  !Intermediate data is not contiguous or there are no intermediates
4281 
4282  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4283  !Rates data is not contiguous or there are no rates
4284 
4285  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,rates,intermediates,parameters)
4286 
4287  !Copy intermediate data from temporary array
4288  DO intermediateidx=1,numberofintermediates
4289  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4290  ENDDO !intermediateIdx
4291 
4292  !Copy rate data from temporary array
4293  DO rateidx=1,numberofstates
4294  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4295  ENDDO !rateIdx
4296 
4297  ELSE
4298  !Rates data is contiguous
4299 
4300  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4301  rateenddof=ratestartdof+numberofstates-1
4302 
4303  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,ratedata(ratestartdof:rateenddof), &
4304  & intermediates,parameters)
4305 
4306  !Copy intermediate data from temporary array
4307  DO intermediateidx=1,numberofintermediates
4308  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4309  ENDDO !intermediateIdx
4310 
4311  ENDIF
4312 
4313  ELSE
4314  !Intermediate data is contiguous
4315 
4316  intermediatestartdof=(intermediatestartidx-1)*intermediatedataoffset+1
4317  intermediateenddof=intermediatestartdof+numberofintermediates-1
4318 
4319  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4320  !Rates data is not contiguous or there are no rates
4321 
4322  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,rates, &
4323  & intermediatedata(intermediatestartdof:intermediateenddof),parameters)
4324 
4325  !Copy rate data from temporary array
4326  DO rateidx=1,numberofstates
4327  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4328  ENDDO !rateIdx
4329 
4330  ELSE
4331  !Rates data is contiguous
4332 
4333  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4334  rateenddof=ratestartdof+numberofstates-1
4335 
4336  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,ratedata(ratestartdof:rateenddof), &
4337  & intermediatedata(intermediatestartdof:intermediateenddof),parameters)
4338 
4339  ENDIF
4340  ENDIF
4341  ELSE
4342  !Parameters data is contiguous
4343 
4344  parameterstartdof=(parameterstartidx-1)*parameterdataoffset+1
4345  parameterenddof=parameterstartdof+numberofparameters-1
4346 
4347  IF(intermediatedataoffset>1.OR.numberofintermediates==0) THEN
4348  !Intermediate data is not contiguous or there are no intermediates
4349 
4350  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4351  !Rates data is not contiguous or there are no rates
4352 
4353  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,rates,intermediates, &
4354  & parameters(parameterstartdof:parameterenddof))
4355 
4356  !Copy intermediate data from temporary array
4357  DO intermediateidx=1,numberofintermediates
4358  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4359  ENDDO !intermediateIdx
4360 
4361  !Copy rate data from temporary array
4362  DO rateidx=1,numberofstates
4363  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4364  ENDDO !rateIdx
4365 
4366  ELSE
4367  !Rates data is contiguous
4368 
4369  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4370  rateenddof=ratestartdof+numberofstates-1
4371 
4372  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,ratedata(ratestartdof:rateenddof), &
4373  & intermediates,parameters(parameterstartdof:parameterenddof))
4374 
4375  !Copy intermediate data from temporary array
4376  DO intermediateidx=1,numberofintermediates
4377  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4378  ENDDO !intermediateIdx
4379 
4380  ENDIF
4381 
4382  ELSE
4383  !Intermediate data is contiguous
4384 
4385  intermediatestartdof=(intermediatestartidx-1)*intermediatedataoffset+1
4386  intermediateenddof=intermediatestartdof+numberofintermediates-1
4387 
4388  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4389  !Rates data is not contiguous or there are no rates
4390 
4391  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,rates, &
4392  & intermediatedata(intermediatestartdof:intermediateenddof), &
4393  & parameters(parameterstartdof:parameterenddof))
4394 
4395  !Copy rate data from temporary array
4396  DO rateidx=1,numberofstates
4397  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4398  ENDDO !rateIdx
4399 
4400  ELSE
4401  !Rates data is contiguous
4402 
4403  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4404  rateenddof=ratestartdof+numberofstates-1
4405 
4406  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states,ratedata(ratestartdof:rateenddof), &
4407  & intermediatedata(intermediatestartdof:intermediateenddof), &
4408  & parameters(parameterstartdof:parameterenddof))
4409 
4410  ENDIF
4411  ENDIF
4412  ENDIF
4413  ELSE
4414  !State data is contiguous
4415 
4416  statestartdof=(statestartidx-1)*statedataoffset+1
4417  stateenddof=statestartdof+numberofstates-1
4418 
4419  IF(parameterdataoffset>1.OR.numberofparameters==0) THEN
4420  !Parameter data is not contiguous or there are no parameters
4421 
4422  !Copy parameter data to temporary array
4423  DO parameteridx=1,numberofparameters
4424  parameters(parameteridx)=parameterdata((parameterstartidx-1)*parameterdataoffset+parameteridx)
4425  ENDDO !parameterIdx
4426 
4427  IF(intermediatedataoffset>1.OR.numberofintermediates==0) THEN
4428  !Intermediate data is not contiguous or there are no intermediates
4429 
4430  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4431  !Rates data is not contiguous or there are no rates
4432 
4433  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof), &
4434  & rates,intermediates,parameters)
4435 
4436  !Copy intermediate data from temporary array
4437  DO intermediateidx=1,numberofintermediates
4438  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4439  ENDDO !intermediateIdx
4440 
4441  !Copy rate data from temporary array
4442  DO rateidx=1,numberofstates
4443  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4444  ENDDO !rateIdx
4445 
4446  ELSE
4447  !Rates data is contiguous
4448 
4449  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4450  rateenddof=ratestartdof+numberofstates-1
4451 
4452  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof), &
4453  & ratedata(ratestartdof:rateenddof),intermediates,parameters)
4454 
4455  !Copy intermediate data from temporary array
4456  DO intermediateidx=1,numberofintermediates
4457  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4458  ENDDO !intermediateIdx
4459 
4460  ENDIF
4461 
4462  ELSE
4463  !Intermediate data is contiguous
4464 
4465  intermediatestartdof=(intermediatestartidx-1)*intermediatedataoffset+1
4466  intermediateenddof=intermediatestartdof+numberofintermediates-1
4467 
4468  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4469  !Rates data is not contiguous or there are no rates
4470 
4471  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof),rates, &
4472  & intermediatedata(intermediatestartdof:intermediateenddof),parameters)
4473 
4474  !Copy rate data from temporary array
4475  DO rateidx=1,numberofstates
4476  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4477  ENDDO !rateIdx
4478 
4479  ELSE
4480  !Rates data is contiguous
4481 
4482  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4483  rateenddof=ratestartdof+numberofstates-1
4484 
4485  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof), &
4486  & ratedata(ratestartdof:rateenddof),intermediatedata(intermediatestartdof:intermediateenddof), &
4487  & parameters)
4488 
4489  ENDIF
4490  ENDIF
4491  ELSE
4492  !Parameters data is contiguous
4493 
4494  parameterstartdof=(parameterstartidx-1)*parameterdataoffset+1
4495  parameterenddof=parameterstartdof+numberofparameters-1
4496 
4497  IF(intermediatedataoffset>1.OR.numberofintermediates==0) THEN
4498  !Intermediate data is not contiguous or there are no intermediates
4499 
4500  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4501  !Rates data is not contiguous or there are no rates
4502 
4503  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof), &
4504  & rates,intermediates,parameters(parameterstartdof:parameterenddof))
4505 
4506  !Copy intermediate data from temporary array
4507  DO intermediateidx=1,numberofintermediates
4508  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4509  ENDDO !intermediateIdx
4510 
4511  !Copy rate data from temporary array
4512  DO rateidx=1,numberofstates
4513  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4514  ENDDO !rateIdx
4515 
4516  ELSE
4517  !Rates data is contiguous
4518 
4519  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4520  rateenddof=ratestartdof+numberofstates-1
4521 
4522  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof), &
4523  & ratedata(ratestartdof:rateenddof),intermediates,parameters(parameterstartdof:parameterenddof))
4524 
4525  !Copy intermediate data from temporary array
4526  DO intermediateidx=1,numberofintermediates
4527  intermediatedata((intermediatestartidx-1)*intermediatedataoffset+intermediateidx)=intermediates(intermediateidx)
4528  ENDDO !intermediateIdx
4529 
4530  ENDIF
4531 
4532  ELSE
4533  !Intermediate data is contiguous
4534 
4535  intermediatestartdof=(intermediatestartidx-1)*intermediatedataoffset+1
4536  intermediateenddof=intermediatestartdof+numberofintermediates-1
4537 
4538  IF(ratedataoffset>1.OR.numberofstates==0) THEN
4539  !Rates data is not contiguous or there are no rates
4540 
4541  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof), &
4542  & rates,intermediatedata(intermediatestartdof:intermediateenddof), &
4543  & parameters(parameterstartdof:parameterenddof))
4544 
4545  !Copy rate data from temporary array
4546  DO rateidx=1,numberofstates
4547  ratedata((ratestartidx-1)*ratedataoffset+rateidx)=rates(rateidx)
4548  ENDDO !rateIdx
4549 
4550  ELSE
4551  !Rates data is contiguous
4552 
4553  ratestartdof=(ratestartidx-1)*ratedataoffset+1
4554  rateenddof=ratestartdof+numberofstates-1
4555 
4556  CALL cellml_model_definition_call_rhs_routine(model%ptr,time,states(statestartdof:stateenddof), &
4557  & ratedata(ratestartdof:rateenddof),intermediatedata(intermediatestartdof:intermediateenddof), &
4558  & parameters(parameterstartdof:parameterenddof))
4559 
4560  ENDIF
4561  ENDIF
4562  ENDIF
4563  ENDIF
4564  ELSE
4565  CALL flagerror("Model is not associated.",err,error,*999)
4566  ENDIF
4567 
4568 #else
4569  CALL flagerror("Must compile with WITH_CELLML ON to use CellML functionality.",err,error,*999)
4570 #endif
4571 
4572  exits("Solver_DAECellMLRHSEvaluate")
4573  RETURN
4574 999 errorsexits("Solver_DAECellMLRHSEvaluate",err,error)
4575  RETURN 1
4576 
4577  END SUBROUTINE solver_daecellmlrhsevaluate
4578 
4579  !
4580  !================================================================================================================================
4581  !
4582 
4584  SUBROUTINE solver_dae_runge_kutta_finalise(RUNGE_KUTTA_SOLVER,ERR,ERROR,*)
4586  !Argument variables
4587  TYPE(runge_kutta_dae_solver_type), POINTER :: RUNGE_KUTTA_SOLVER
4588  INTEGER(INTG), INTENT(OUT) :: ERR
4589  TYPE(varying_string), INTENT(OUT) :: ERROR
4590  !Local Variables
4591 
4592  enters("SOLVER_DAE_RUNGE_KUTTA_FINALISE",err,error,*999)
4593 
4594  IF(ASSOCIATED(runge_kutta_solver)) THEN
4595  DEALLOCATE(runge_kutta_solver)
4596  ENDIF
4597 
4598  exits("SOLVER_DAE_RUNGE_KUTTA_FINALISE")
4599  RETURN
4600 999 errorsexits("SOLVER_DAE_RUNGE_KUTTA_FINALISE",err,error)
4601  RETURN 1
4602 
4603  END SUBROUTINE solver_dae_runge_kutta_finalise
4604 
4605  !
4606  !================================================================================================================================
4607  !
4608 
4610  SUBROUTINE solver_dae_runge_kutta_initialise(DAE_SOLVER,ERR,ERROR,*)
4612  !Argument variables
4613  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
4614  INTEGER(INTG), INTENT(OUT) :: ERR
4615  TYPE(varying_string), INTENT(OUT) :: ERROR
4616  !Local Variables
4617  INTEGER(INTG) :: DUMMY_ERR
4618  TYPE(varying_string) :: DUMMY_ERROR
4619 
4620  enters("SOLVER_DAE_RUNGE_KUTTA_INITIALISE",err,error,*998)
4621 
4622  IF(ASSOCIATED(dae_solver)) THEN
4623  IF(ASSOCIATED(dae_solver%RUNGE_KUTTA_SOLVER)) THEN
4624  CALL flagerror("Runge-Kutta solver is already associated for this differential-algebraic equation solver.",err,error,*998)
4625  ELSE
4626  !Allocate the Runge-Kutta solver
4627  ALLOCATE(dae_solver%RUNGE_KUTTA_SOLVER,stat=err)
4628  IF(err/=0) CALL flagerror("Could not allocate Runge-Kutta solver.",err,error,*999)
4629  !Initialise
4630  dae_solver%RUNGE_KUTTA_SOLVER%DAE_SOLVER=>dae_solver
4631  dae_solver%RUNGE_KUTTA_SOLVER%SOLVER_LIBRARY=0
4632  !Defaults
4633  ENDIF
4634  ELSE
4635  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*998)
4636  ENDIF
4637 
4638  exits("SOLVER_DAE_RUNGE_KUTTA_INITIALISE")
4639  RETURN
4640 999 CALL solver_dae_runge_kutta_finalise(dae_solver%RUNGE_KUTTA_SOLVER,dummy_err,dummy_error,*998)
4641 998 errorsexits("SOLVER_DAE_RUNGE_KUTTA_INITIALISE",err,error)
4642  RETURN 1
4643 
4644  END SUBROUTINE solver_dae_runge_kutta_initialise
4645 
4646  !
4647  !================================================================================================================================
4648  !
4649 
4651  SUBROUTINE solver_dae_runge_kutta_solve(RUNGE_KUTTA_SOLVER,ERR,ERROR,*)
4653  !Argument variables
4654  TYPE(runge_kutta_dae_solver_type), POINTER :: RUNGE_KUTTA_SOLVER
4655  INTEGER(INTG), INTENT(OUT) :: ERR
4656  TYPE(varying_string), INTENT(OUT) :: ERROR
4657  !Local Variables
4658 
4659  enters("SOLVER_DAE_RUNGE_KUTTA_SOLVE",err,error,*999)
4660 
4661  IF(ASSOCIATED(runge_kutta_solver)) THEN
4662  CALL flagerror("Not implemented.",err,error,*999)
4663  ELSE
4664  CALL flagerror("Runge-Kutta differential-algebraic equation solver is not associated.",err,error,*999)
4665  ENDIF
4666 
4667  exits("SOLVER_DAE_RUNGE_KUTTA_SOLVE")
4668  RETURN
4669 999 errorsexits("SOLVER_DAE_RUNGE_KUTTA_SOLVE",err,error)
4670  RETURN 1
4671 
4672  END SUBROUTINE solver_dae_runge_kutta_solve
4673 
4674  !
4675  !================================================================================================================================
4676  !
4677 
4679  SUBROUTINE solver_dae_rush_larson_finalise(RUSH_LARSON_SOLVER,ERR,ERROR,*)
4681  !Argument variables
4682  TYPE(rush_larson_dae_solver_type), POINTER :: RUSH_LARSON_SOLVER
4683  INTEGER(INTG), INTENT(OUT) :: ERR
4684  TYPE(varying_string), INTENT(OUT) :: ERROR
4685  !Local Variables
4686 
4687  enters("SOLVER_DAE_RUSH_LARSON_FINALISE",err,error,*999)
4688 
4689  IF(ASSOCIATED(rush_larson_solver)) THEN
4690  DEALLOCATE(rush_larson_solver)
4691  ENDIF
4692 
4693  exits("SOLVER_DAE_RUSH_LARSON_FINALISE")
4694  RETURN
4695 999 errorsexits("SOLVER_DAE_RUSH_LARSON_FINALISE",err,error)
4696  RETURN 1
4697 
4698  END SUBROUTINE solver_dae_rush_larson_finalise
4699 
4700  !
4701  !================================================================================================================================
4702  !
4703 
4705  SUBROUTINE solver_dae_rush_larson_initialise(DAE_SOLVER,ERR,ERROR,*)
4707  !Argument variables
4708  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
4709  INTEGER(INTG), INTENT(OUT) :: ERR
4710  TYPE(varying_string), INTENT(OUT) :: ERROR
4711  !Local Variables
4712  INTEGER(INTG) :: DUMMY_ERR
4713  TYPE(varying_string) :: DUMMY_ERROR
4714 
4715  enters("SOLVER_DAE_RUSH_LARSON_INITIALISE",err,error,*998)
4716 
4717  IF(ASSOCIATED(dae_solver)) THEN
4718  IF(ASSOCIATED(dae_solver%RUSH_LARSON_SOLVER)) THEN
4719  CALL flagerror("Rush-Larson solver is already associated for this differential-algebraic equation solver.",err,error,*998)
4720  ELSE
4721  !Allocate the Rush-Larson solver
4722  ALLOCATE(dae_solver%RUSH_LARSON_SOLVER,stat=err)
4723  IF(err/=0) CALL flagerror("Could not allocate Rush-Larson solver.",err,error,*999)
4724  !Initialise
4725  dae_solver%RUSH_LARSON_SOLVER%DAE_SOLVER=>dae_solver
4726  dae_solver%RUSH_LARSON_SOLVER%SOLVER_LIBRARY=0
4727  !Defaults
4728  ENDIF
4729  ELSE
4730  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*998)
4731  ENDIF
4732 
4733  exits("SOLVER_DAE_RUSH_LARSON_INITIALISE")
4734  RETURN
4735 999 CALL solver_dae_rush_larson_finalise(dae_solver%RUSH_LARSON_SOLVER,dummy_err,dummy_error,*998)
4736 998 errorsexits("SOLVER_DAE_RUSH_LARSON_INITIALISE",err,error)
4737  RETURN 1
4738 
4739  END SUBROUTINE solver_dae_rush_larson_initialise
4740 
4741  !
4742  !================================================================================================================================
4743  !
4744 
4746  SUBROUTINE solver_dae_rush_larson_solve(RUSH_LARSON_SOLVER,ERR,ERROR,*)
4748  !Argument variables
4749  TYPE(rush_larson_dae_solver_type), POINTER :: RUSH_LARSON_SOLVER
4750  INTEGER(INTG), INTENT(OUT) :: ERR
4751  TYPE(varying_string), INTENT(OUT) :: ERROR
4752  !Local Variables
4753 
4754  enters("SOLVER_DAE_RUSH_LARSON_SOLVE",err,error,*999)
4755 
4756  IF(ASSOCIATED(rush_larson_solver)) THEN
4757  CALL flagerror("Not implemented.",err,error,*999)
4758  ELSE
4759  CALL flagerror("Rush-Larson differential-algebraic equation solver is not associated.",err,error,*999)
4760  ENDIF
4761 
4762  exits("SOLVER_DAE_RUSH_LARSON_SOLVE")
4763  RETURN
4764 999 errorsexits("SOLVER_DAE_RUSH_LARSON_SOLVE",err,error)
4765  RETURN 1
4766 
4767  END SUBROUTINE solver_dae_rush_larson_solve
4768 
4769  !
4770  !================================================================================================================================
4771  !
4772 
4774  SUBROUTINE solver_dae_solve(DAE_SOLVER,ERR,ERROR,*)
4776  !Argument variables
4777  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
4778  INTEGER(INTG), INTENT(OUT) :: ERR
4779  TYPE(varying_string), INTENT(OUT) :: ERROR
4780  !Local Variables
4781  INTEGER(INTG) :: cellml_idx
4782  TYPE(cellml_type), POINTER :: CELLML
4783  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
4784  TYPE(cellml_state_field_type), POINTER :: CELLML_STATE_FIELD
4785  TYPE(solver_type), POINTER :: SOLVER
4786  TYPE(varying_string) :: LOCAL_ERROR
4787 
4788  enters("SOLVER_DAE_SOLVE",err,error,*999)
4789 
4790  IF(ASSOCIATED(dae_solver)) THEN
4791  solver=>dae_solver%SOLVER
4792  IF(ASSOCIATED(solver)) THEN
4793  SELECT CASE(dae_solver%DAE_SOLVE_TYPE)
4794  CASE(solver_dae_euler)
4795  CALL solver_dae_euler_solve(dae_solver%EULER_SOLVER,err,error,*999)
4797  CALL solver_dae_crank_nicolson_solve(dae_solver%CRANK_NICOLSON_SOLVER,err,error,*999)
4799  CALL solver_dae_runge_kutta_solve(dae_solver%RUNGE_KUTTA_SOLVER,err,error,*999)
4801  CALL solver_dae_adams_moulton_solve(dae_solver%ADAMS_MOULTON_SOLVER,err,error,*999)
4802  CASE(solver_dae_bdf)
4803  CALL solver_dae_bdf_solve(dae_solver%BDF_SOLVER,err,error,*999)
4805  CALL solver_dae_rush_larson_solve(dae_solver%RUSH_LARSON_SOLVER,err,error,*999)
4806  CASE(solver_dae_external)
4807  CALL solver_dae_external_solve(dae_solver%EXTERNAL_SOLVER,err,error,*999)
4808  CASE DEFAULT
4809  local_error="The differential-algebraic equation solver solve type of "// &
4810  & trim(numbertovstring(dae_solver%DAE_SOLVE_TYPE,"*",err,error))//" is invalid."
4811  CALL flagerror(local_error,err,error,*999)
4812  END SELECT
4813  IF(solver%OUTPUT_TYPE>solver_solver_output) THEN
4814 #ifdef TAUPROF
4815  CALL tau_static_phase_start("Solution Output Phase")
4816 #endif
4817  cellml_equations=>solver%CELLML_EQUATIONS
4818  IF(ASSOCIATED(cellml_equations)) THEN
4819  CALL write_string(general_output_type,"",err,error,*999)
4820  CALL write_string(general_output_type,"Solver State vectors:",err,error,*999)
4821  CALL write_string_value(general_output_type,"Number of CellML environments = ",cellml_equations% &
4822  & number_of_cellml_environments,err,error,*999)
4823  DO cellml_idx=1,cellml_equations%NUMBER_OF_CELLML_ENVIRONMENTS
4824  cellml=>cellml_equations%CELLML_ENVIRONMENTS(cellml_idx)%PTR
4825  IF(ASSOCIATED(cellml)) THEN
4826  cellml_state_field=>cellml%STATE_FIELD
4827  IF(ASSOCIATED(cellml_state_field)) THEN
4828  CALL write_string_value(general_output_type,"CellML index : ",cellml_idx,err,error,*999)
4829  CALL field_parameter_set_output(general_output_type,cellml_state_field%STATE_FIELD,field_u_variable_type, &
4830  & field_values_set_type,err,error,*999)
4831  ELSE
4832  CALL flagerror("CellML environment state field is not associated.",err,error,*999)
4833  ENDIF
4834  ELSE
4835  local_error="CellML environment is not associated for CellML index "// &
4836  & trim(numbertovstring(cellml_idx,"*",err,error))//"."
4837  CALL flagerror(local_error,err,error,*999)
4838  ENDIF
4839  ENDDO !cellml_idx
4840 
4841  ELSE
4842  CALL flagerror("Solver CellML equations is not associated.",err,error,*999)
4843  ENDIF
4844 #ifdef TAUPROF
4845  CALL tau_static_phase_stop("Solution Output Phase")
4846 #endif
4847  ENDIF
4848  ELSE
4849  CALL flagerror("Differential-algebraic solver solver is not associated.",err,error,*999)
4850  ENDIF
4851  ELSE
4852  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*999)
4853  ENDIF
4854 
4855  exits("SOLVER_DAE_SOLVE")
4856  RETURN
4857 999 errorsexits("SOLVER_DAE_SOLVE",err,error)
4858  RETURN 1
4859 
4860  END SUBROUTINE solver_dae_solve
4861 
4862  !
4863  !================================================================================================================================
4864  !
4865 
4867  SUBROUTINE solver_dae_solver_type_get(SOLVER,DAE_SOLVE_TYPE,ERR,ERROR,*)
4869  !Argument variables
4870  TYPE(solver_type), POINTER :: SOLVER
4871  INTEGER(INTG), INTENT(OUT) :: DAE_SOLVE_TYPE
4872  INTEGER(INTG), INTENT(OUT) :: ERR
4873  TYPE(varying_string), INTENT(OUT) :: ERROR
4874  !Local Variables
4875  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
4876 
4877  enters("SOLVER_DAE_SOLVER_TYPE_GET",err,error,*999)
4878 
4879  IF(ASSOCIATED(solver)) THEN
4880  IF(solver%SOLVER_FINISHED) THEN
4881  IF(solver%SOLVE_TYPE==solver_dae_type) THEN
4882  dae_solver=>solver%DAE_SOLVER
4883  IF(ASSOCIATED(dae_solver)) THEN
4884  dae_solve_type=dae_solver%DAE_SOLVE_TYPE
4885  ELSE
4886  CALL flagerror("The solver differential-algebraic equation solver is not associated.",err,error,*999)
4887  ENDIF
4888  ELSE
4889  CALL flagerror("The solver is not a differential-algebraic equation solver.",err,error,*999)
4890  ENDIF
4891  ELSE
4892  CALL flagerror("Solver has not been finished.",err,error,*999)
4893  ENDIF
4894  ELSE
4895  CALL flagerror("Solver is not associated.",err,error,*999)
4896  ENDIF
4897 
4898  exits("SOLVER_DAE_SOLVER_TYPE_GET")
4899  RETURN
4900 999 errorsexits("SOLVER_DAE_SOLVER_TYPE_GET",err,error)
4901  RETURN 1
4902 
4903  END SUBROUTINE solver_dae_solver_type_get
4904 
4905  !
4906  !================================================================================================================================
4907  !
4908 
4910  SUBROUTINE solver_dae_solver_type_set(SOLVER,DAE_SOLVE_TYPE,ERR,ERROR,*)
4912  !Argument variables
4913  TYPE(solver_type), POINTER :: SOLVER
4914  INTEGER(INTG), INTENT(IN) :: DAE_SOLVE_TYPE
4915  INTEGER(INTG), INTENT(OUT) :: ERR
4916  TYPE(varying_string), INTENT(OUT) :: ERROR
4917  !Local Variables
4918  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
4919  TYPE(varying_string) :: LOCAL_ERROR
4920 
4921  enters("SOLVER_DAE_SOLVER_TYPE_SET",err,error,*999)
4922 
4923  IF(ASSOCIATED(solver)) THEN
4924  IF(solver%SOLVER_FINISHED) THEN
4925  CALL flagerror("Solver has already been finished.",err,error,*999)
4926  ELSE
4927  IF(solver%SOLVE_TYPE==solver_dae_type) THEN
4928  dae_solver=>solver%DAE_SOLVER
4929  IF(ASSOCIATED(dae_solver)) THEN
4930  IF(dae_solve_type/=dae_solver%DAE_SOLVE_TYPE) THEN
4931  !Intialise the new differential-algebraic equation solver type
4932  SELECT CASE(dae_solve_type)
4933  CASE(solver_dae_euler)
4934  CALL solver_dae_euler_initialise(dae_solver,err,error,*999)
4936  CALL solver_dae_crank_nicolson_initialise(dae_solver,err,error,*999)
4938  CALL solver_dae_runge_kutta_initialise(dae_solver,err,error,*999)
4940  CALL solver_dae_adams_moulton_initialise(dae_solver,err,error,*999)
4941  CASE(solver_dae_bdf)
4942  CALL solver_dae_bdf_initialise(dae_solver,err,error,*999)
4944  CALL solver_dae_rush_larson_initialise(dae_solver,err,error,*999)
4945  CASE(solver_dae_external)
4946  CALL solver_dae_external_initialise(dae_solver,err,error,*999)
4947  CASE DEFAULT
4948  local_error="The specified differential-algebraic equation solver type of "// &
4949  & trim(numbertovstring(dae_solve_type,"*",err,error))//" is invalid."
4950  CALL flagerror(local_error,err,error,*999)
4951  END SELECT
4952  !Finalise the old differential-algebraic equation solver type
4953  SELECT CASE(dae_solver%DAE_SOLVE_TYPE)
4954  CASE(solver_dae_euler)
4955  CALL solver_dae_euler_finalise(dae_solver%EULER_SOLVER,err,error,*999)
4957  CALL solver_dae_crank_nicolson_finalise(dae_solver%CRANK_NICOLSON_SOLVER,err,error,*999)
4959  CALL solver_dae_runge_kutta_finalise(dae_solver%RUNGE_KUTTA_SOLVER,err,error,*999)
4961  CALL solver_dae_adams_moulton_finalise(dae_solver%ADAMS_MOULTON_SOLVER,err,error,*999)
4962  CASE(solver_dae_bdf)
4963  CALL solver_dae_bdf_finalise(dae_solver%BDF_SOLVER,err,error,*999)
4965  CALL solver_dae_rush_larson_finalise(dae_solver%RUSH_LARSON_SOLVER,err,error,*999)
4966  CASE(solver_dae_external)
4967  CALL solver_dae_external_finalise(dae_solver%EXTERNAL_SOLVER,err,error,*999)
4968  CASE DEFAULT
4969  local_error="The differential-algebraic equation solve type of "// &
4970  & trim(numbertovstring(dae_solver%DAE_SOLVE_TYPE,"*",err,error))//" is invalid."
4971  CALL flagerror(local_error,err,error,*999)
4972  END SELECT
4973  dae_solver%DAE_SOLVE_TYPE=dae_solve_type
4974  ENDIF
4975  ELSE
4976  CALL flagerror("The solver differential-algebraic equation solver is not associated.",err,error,*999)
4977  ENDIF
4978  ELSE
4979  CALL flagerror("The solver is not a differential-algebraic equation solver.",err,error,*999)
4980  ENDIF
4981  ENDIF
4982  ELSE
4983  CALL flagerror("Solver is not associated.",err,error,*999)
4984  ENDIF
4985 
4986  exits("SOLVER_DAE_SOLVER_TYPE_SET")
4987  RETURN
4988 999 errorsexits("SOLVER_DAE_SOLVER_TYPE_SET",err,error)
4989  RETURN 1
4990 
4991  END SUBROUTINE solver_dae_solver_type_set
4992 
4993  !
4994  !================================================================================================================================
4995  !
4996 
4998  SUBROUTINE solver_dae_times_set(SOLVER,START_TIME,END_TIME,ERR,ERROR,*)
5000  !Argument variables
5001  TYPE(solver_type), POINTER :: SOLVER
5002  REAL(DP), INTENT(IN) :: START_TIME
5003  REAL(DP), INTENT(IN) :: END_TIME
5004  INTEGER(INTG), INTENT(OUT) :: ERR
5005  TYPE(varying_string), INTENT(OUT) :: ERROR
5006  !Local Variables
5007  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
5008  TYPE(varying_string) :: LOCAL_ERROR
5009 
5010  enters("SOLVER_DAE_TIMES_SET",err,error,*999)
5011 
5012  IF(ASSOCIATED(solver)) THEN
5013  IF(solver%SOLVE_TYPE==solver_dae_type) THEN
5014  dae_solver=>solver%DAE_SOLVER
5015  IF(ASSOCIATED(dae_solver)) THEN
5016  IF(end_time>start_time) THEN
5017  dae_solver%START_TIME=start_time
5018  dae_solver%END_TIME=end_time
5019  ELSE
5020  local_error="The specified end time of "//trim(numbertovstring(end_time,"*",err,error))// &
5021  & " is not > than the specified start time of "//trim(numbertovstring(start_time,"*",err,error))//"."
5022  CALL flagerror(local_error,err,error,*999)
5023  ENDIF
5024  ELSE
5025  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*999)
5026  ENDIF
5027  ELSE
5028  CALL flagerror("The solver is not a differential-algebraic equation solver.",err,error,*999)
5029  ENDIF
5030  ELSE
5031  CALL flagerror("Solver is not associated.",err,error,*999)
5032  ENDIF
5033 
5034  exits("SOLVER_DAE_TIMES_SET")
5035  RETURN
5036 999 errorsexits("SOLVER_DAE_TIMES_SET",err,error)
5037  RETURN 1
5038 
5039  END SUBROUTINE solver_dae_times_set
5040 
5041  !
5042  !================================================================================================================================
5043  !
5044 
5046  SUBROUTINE solver_dae_time_step_set(SOLVER,TIME_STEP,ERR,ERROR,*)
5048  !Argument variables
5049  TYPE(solver_type), POINTER :: SOLVER
5050  REAL(DP), INTENT(IN) :: TIME_STEP
5051  INTEGER(INTG), INTENT(OUT) :: ERR
5052  TYPE(varying_string), INTENT(OUT) :: ERROR
5053  !Local Variables
5054  TYPE(dae_solver_type), POINTER :: DAE_SOLVER
5055  TYPE(varying_string) :: LOCAL_ERROR
5056 
5057  enters("SOLVER_DAE_TIME_STEP_SET",err,error,*999)
5058 
5059  IF(ASSOCIATED(solver)) THEN
5060  IF(solver%SOLVE_TYPE==solver_dae_type) THEN
5061  dae_solver=>solver%DAE_SOLVER
5062  IF(ASSOCIATED(dae_solver)) THEN
5063  IF(abs(time_step)<=zero_tolerance) THEN
5064  local_error="The specified time step of "//trim(numbertovstring(time_step,"*",err,error))// &
5065  & " is invalid. The time step must not be zero."
5066  CALL flagerror(local_error,err,error,*999)
5067  ELSE
5068  dae_solver%INITIAL_STEP=time_step
5069  ENDIF
5070  ELSE
5071  CALL flagerror("Differential-algebraic equation solver is not associated.",err,error,*999)
5072  ENDIF
5073  ELSE
5074  CALL flagerror("The solver is not a differential-algebraic equation solver.",err,error,*999)
5075  ENDIF
5076  ELSE
5077  CALL flagerror("Solver is not associated.",err,error,*999)
5078  ENDIF
5079 
5080  exits("SOLVER_DAE_TIME_STEP_SET")
5081  RETURN
5082 999 errorsexits("SOLVER_DAE_TIME_STEP_SET",err,error)
5083  RETURN 1
5084 
5085  END SUBROUTINE solver_dae_time_step_set
5086 
5087  !
5088  !================================================================================================================================
5089  !
5090 
5092  SUBROUTINE solver_destroy(SOLVER,ERR,ERROR,*)
5094  !Argument variables
5095  TYPE(solver_type), POINTER :: SOLVER
5096  INTEGER(INTG), INTENT(OUT) :: ERR
5097  TYPE(varying_string), INTENT(OUT) :: ERROR
5098  !Local Variables
5099 
5100  enters("SOLVER_DESTROY",err,error,*999)
5101 
5102  IF(ASSOCIATED(solver)) THEN
5103  CALL flagerror("Not implemented.",err,error,*999)
5104  ELSE
5105  CALL flagerror("Solver is not associated.",err,error,*999)
5106  ENDIF
5107 
5108  exits("SOLVER_DESTROY")
5109  RETURN
5110 999 errorsexits("SOLVER_DESTROY",err,error)
5111  RETURN 1
5112 
5113  END SUBROUTINE solver_destroy
5114 
5115  !
5116  !================================================================================================================================
5117  !
5118 
5120  SUBROUTINE solver_dynamic_create_finish(DYNAMIC_SOLVER,ERR,ERROR,*)
5122  !Argument variables
5123  TYPE(dynamic_solver_type), POINTER :: DYNAMIC_SOLVER
5124  INTEGER(INTG), INTENT(OUT) :: ERR
5125  TYPE(varying_string), INTENT(OUT) :: ERROR
5126  !Local Variables
5127  INTEGER(INTG) :: DYNAMIC_VARIABLE_TYPE,equations_matrix_idx,equations_set_idx,LINEAR_LIBRARY_TYPE,NONLINEAR_LIBRARY_TYPE
5128  INTEGER(INTG) :: VariableType=0
5129  TYPE(equations_type), POINTER :: EQUATIONS
5130  TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
5131  TYPE(equations_mapping_dynamic_type), POINTER :: DYNAMIC_MAPPING
5132  TYPE(equations_mapping_nonlinear_type), POINTER :: NonlinearMapping
5133  TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
5134  TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
5135  TYPE(equations_matrices_dynamic_type), POINTER :: DYNAMIC_MATRICES
5136  TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
5137  TYPE(equations_matrix_type), POINTER :: DAMPING_MATRIX,EQUATIONS_MATRIX,MASS_MATRIX
5138  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
5139  TYPE(field_type), POINTER :: DEPENDENT_FIELD !, INDEPENDENT_FIELD
5140  TYPE(field_variable_type), POINTER :: DYNAMIC_VARIABLE,LINEAR_VARIABLE,ResidualVariable
5141  TYPE(solver_type), POINTER :: SOLVER,LINEAR_SOLVER,NONLINEAR_SOLVER
5142  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
5143  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
5144  TYPE(solver_matrices_type), POINTER :: SOLVER_MATRICES
5145  TYPE(varying_string) :: LOCAL_ERROR
5146 
5147  enters("SOLVER_DYNAMIC_CREATE_FINISH",err,error,*999)
5148 
5149  IF(ASSOCIATED(dynamic_solver)) THEN
5150  solver=>dynamic_solver%SOLVER
5151  IF(ASSOCIATED(solver)) THEN
5152  solver_equations=>solver%SOLVER_EQUATIONS
5153  IF(ASSOCIATED(solver_equations)) THEN
5154  SELECT CASE(dynamic_solver%SOLVER_LIBRARY)
5155  CASE(solver_cmiss_library)
5156  !Create the parameter sets required for the solver
5157  solver_equations=>solver%SOLVER_EQUATIONS
5158  IF(ASSOCIATED(solver_equations)) THEN
5159  solver_mapping=>solver_equations%SOLVER_MAPPING
5160  IF(ASSOCIATED(solver_mapping)) THEN
5161  !Initialise for explicit solve
5162  dynamic_solver%EXPLICIT=abs(dynamic_solver%THETA(dynamic_solver%DEGREE))<zero_tolerance
5163  !Loop over the equations set in the solver equations
5164  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
5165  equations=>solver_mapping%EQUATIONS_SET_TO_SOLVER_MAP(equations_set_idx)%EQUATIONS
5166  IF(ASSOCIATED(equations)) THEN
5167  equations_set=>equations%EQUATIONS_SET
5168  IF(ASSOCIATED(equations_set)) THEN
5169  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
5170  IF(ASSOCIATED(dependent_field)) THEN
5171  equations_mapping=>equations%EQUATIONS_MAPPING
5172  IF(ASSOCIATED(equations_mapping)) THEN
5173  dynamic_mapping=>equations_mapping%DYNAMIC_MAPPING
5174  IF(ASSOCIATED(dynamic_mapping)) THEN
5175  dynamic_variable=>dynamic_mapping%DYNAMIC_VARIABLE
5176  dynamic_variable_type=dynamic_mapping%DYNAMIC_VARIABLE_TYPE
5177  IF(ASSOCIATED(dynamic_variable)) THEN
5178  !Set up the parameter sets to hold the required solver parameters
5179  !1st degree or higher so set up displacement parameter sets
5180 
5181 
5182 
5183  IF(dynamic_solver%DEGREE>=solver_dynamic_second_degree) THEN
5184  !2nd degree or higher so set up velocity parameter sets
5185  CALL field_parametersetensurecreated(dependent_field,dynamic_variable_type, &
5186  & field_velocity_values_set_type,err,error,*999)
5187  CALL field_parametersetensurecreated(dependent_field,dynamic_variable_type, &
5188  & field_previous_velocity_set_type,err,error,*999)
5189  CALL field_parametersetensurecreated(dependent_field, &
5190  & dynamic_variable_type,field_mean_predicted_velocity_set_type,err,error,*999)
5191  IF(dynamic_solver%DEGREE>=solver_dynamic_third_degree) THEN
5192  !3rd degree or higher so set up acceleration parameter sets
5193  CALL field_parametersetensurecreated(dependent_field,dynamic_variable_type, &
5194  & field_acceleration_values_set_type,err,error,*999)
5195  CALL field_parametersetensurecreated(dependent_field,dynamic_variable_type, &
5196  & field_previous_acceleration_set_type,err,error,*999)
5197  CALL field_parametersetensurecreated( &
5198  & dependent_field,dynamic_variable_type,field_mean_predicted_acceleration_set_type, &
5199  & err,error,*999)
5200  ENDIF
5201  ENDIF
5202 
5203 
5204 
5205  !Create the dynamic matrices temporary vector for matrix-vector products
5206  equations_matrices=>equations%EQUATIONS_MATRICES
5207  IF(ASSOCIATED(equations_matrices)) THEN
5208  dynamic_matrices=>equations_matrices%DYNAMIC_MATRICES
5209  IF(ASSOCIATED(dynamic_matrices)) THEN
5210  IF(.NOT.ASSOCIATED(dynamic_matrices%TEMP_VECTOR)) THEN
5211  CALL distributed_vector_create_start(dynamic_variable%DOMAIN_MAPPING, &
5212  & dynamic_matrices%TEMP_VECTOR,err,error,*999)
5213  CALL distributed_vector_data_type_set(dynamic_matrices%TEMP_VECTOR, &
5214  & distributed_matrix_vector_dp_type,err,error,*999)
5215  CALL distributed_vector_create_finish(dynamic_matrices%TEMP_VECTOR,err,error,*999)
5216  ENDIF
5217  !Check to see if we have an explicit solve
5218  IF(abs(dynamic_solver%THETA(dynamic_solver%DEGREE))<zero_tolerance) THEN
5219  IF(dynamic_mapping%DAMPING_MATRIX_NUMBER/=0) THEN
5220  damping_matrix=>dynamic_matrices%MATRICES(dynamic_mapping%DAMPING_MATRIX_NUMBER)%PTR
5221  IF(ASSOCIATED(damping_matrix)) THEN
5222  dynamic_solver%EXPLICIT=dynamic_solver%EXPLICIT.AND.damping_matrix%LUMPED
5223  ELSE
5224  CALL flagerror("Damping matrix is not associated.",err,error,*999)
5225  ENDIF
5226  ENDIF
5227  IF(dynamic_mapping%MASS_MATRIX_NUMBER/=0) THEN
5228  mass_matrix=>dynamic_matrices%MATRICES(dynamic_mapping%MASS_MATRIX_NUMBER)%PTR
5229  IF(ASSOCIATED(mass_matrix)) THEN
5230  dynamic_solver%EXPLICIT=dynamic_solver%EXPLICIT.AND.mass_matrix%LUMPED
5231  ELSE
5232  CALL flagerror("Mass matrix is not associated.",err,error,*999)
5233  ENDIF
5234  ENDIF
5235  ENDIF
5236  ELSE
5237  CALL flagerror("Equations matrices dynamic matrices are not associated.",err,error,*999)
5238  ENDIF
5239  ELSE
5240  CALL flagerror("Equations equations matrices is not associated.",err,error,*999)
5241  ENDIF
5242  variabletype=dynamic_variable_type
5243  ELSE
5244  CALL flagerror("Dynamic mapping dynamic variable is not associated.",err,error,*999)
5245  ENDIF
5246  ENDIF
5247 
5248 
5249  IF(variabletype==0) THEN
5250  !We now allow for static equation sets for dynamic solvers to be able to couple static eqs - dynamic eqs
5251  nonlinearmapping=>equations_mapping%NONLINEAR_MAPPING
5252  IF(ASSOCIATED(nonlinearmapping)) THEN
5253  IF(dynamic_solver%LINEARITY==solver_dynamic_nonlinear) THEN
5254  !Default to first variable type for now
5255  residualvariable=>nonlinearmapping%RESIDUAL_VARIABLES(1)%PTR
5256  IF(ASSOCIATED(residualvariable)) THEN
5257  variabletype=residualvariable%VARIABLE_TYPE
5258  ELSE
5259  CALL flagerror("Residual variable is not associated.",err,error,*999)
5260  ENDIF
5261  ELSE
5262  local_error="The specified dynamic solver linearity type of "// &
5263  & trim(numbertovstring(dynamic_solver%LINEARITY,"*",err,error))// &
5264  & " is invalid for a nonlinear equations mapping."
5265  CALL flagerror(local_error,err,error,*999)
5266  ENDIF
5267  ENDIF
5268  ENDIF
5269  CALL field_parametersetensurecreated(dependent_field,variabletype, &
5270  & field_previous_values_set_type,err,error,*999)
5271  CALL field_parametersetensurecreated(dependent_field,variabletype, &
5272  & field_mean_predicted_displacement_set_type,err,error,*999)
5273 
5274  CALL field_parametersetensurecreated(dependent_field,variabletype, &
5275  & field_incremental_values_set_type,err,error,*999)
5276  CALL field_parametersetensurecreated(dependent_field,variabletype, &
5277  & field_predicted_displacement_set_type,err,error,*999)
5278  CALL field_parametersetensurecreated(dependent_field,variabletype, &
5279  & field_residual_set_type,err,error,*999)
5280  CALL field_parametersetensurecreated(dependent_field,variabletype, &
5281  & field_previous_residual_set_type,err,error,*999)
5282 
5283 
5284  !Check if there are any linear mappings
5285  linear_mapping=>equations_mapping%LINEAR_MAPPING
5286  IF(ASSOCIATED(linear_mapping)) THEN
5287  !If there are any linear matrices create temporary vector for matrix-vector products
5288  equations_matrices=>equations%EQUATIONS_MATRICES
5289  IF(ASSOCIATED(equations_matrices)) THEN
5290  linear_matrices=>equations_matrices%LINEAR_MATRICES
5291  IF(ASSOCIATED(linear_matrices)) THEN
5292  DO equations_matrix_idx=1,linear_matrices%NUMBER_OF_LINEAR_MATRICES
5293  equations_matrix=>linear_matrices%MATRICES(equations_matrix_idx)%PTR
5294  IF(ASSOCIATED(equations_matrix)) THEN
5295  IF(.NOT.ASSOCIATED(equations_matrix%TEMP_VECTOR)) THEN
5296  linear_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(equations_matrix_idx)%VARIABLE
5297  IF(ASSOCIATED(linear_variable)) THEN
5298  CALL distributed_vector_create_start(linear_variable%DOMAIN_MAPPING, &
5299  & equations_matrix%TEMP_VECTOR,err,error,*999)
5300  CALL distributed_vector_data_type_set(equations_matrix%TEMP_VECTOR, &
5301  & distributed_matrix_vector_dp_type,err,error,*999)
5302  CALL distributed_vector_create_finish(equations_matrix%TEMP_VECTOR,err,error,*999)
5303  ELSE
5304  CALL flagerror("Linear mapping linear variable is not associated.",err,error,*999)
5305  ENDIF
5306  ENDIF
5307  ELSE
5308  CALL flagerror("Equations matrix is not associated.",err,error,*999)
5309  ENDIF
5310  ENDDO !equations_matrix_idx
5311  ELSE
5312  CALL flagerror("Equations matrices linear matrices is not associated.",err,error,*999)
5313  ENDIF
5314  ELSE
5315  CALL flagerror("Equations equations matrices is not associated.",err,error,*999)
5316  ENDIF
5317  ENDIF
5318  ELSE
5319  CALL flagerror("Equations equations mapping is not associated.",err,error,*999)
5320  ENDIF
5321  ELSE
5322  local_error="Equations set dependent field is not associated for equations set index "// &
5323  & trim(numbertovstring(equations_set_idx,"*",err,error))//"."
5324  CALL flagerror(local_error,err,error,*999)
5325  ENDIF
5326  ELSE
5327  local_error="Equations equations set is not associated for equations set index "// &
5328  & trim(numbertovstring(equations_set_idx,"*",err,error))//"."
5329  CALL flagerror(local_error,err,error,*999)
5330  ENDIF
5331  ELSE
5332  local_error="Equations is not associated for equations set index "// &
5333  & trim(numbertovstring(equations_set_idx,"*",err,error))//"."
5334  CALL flagerror(local_error,err,error,*999)
5335  ENDIF
5336  ENDDO !equations_set_idx
5337  !Create the solver matrices and vectors
5338  IF(dynamic_solver%LINEARITY==solver_dynamic_linear) THEN
5339  linear_solver=>dynamic_solver%LINEAR_SOLVER
5340  IF(ASSOCIATED(linear_solver)) THEN
5341  NULLIFY(solver_matrices)
5342  CALL solver_matrices_create_start(solver_equations,solver_matrices,err,error,*999)
5343  CALL solver_matrices_library_type_get(linear_solver,linear_library_type,err,error,*999)
5344  CALL solver_matrices_library_type_set(solver_matrices,linear_library_type,err,error,*999)
5345  IF(dynamic_solver%EXPLICIT) THEN
5347  & err,error,*999)
5348  ELSE
5349  SELECT CASE(solver_equations%SPARSITY_TYPE)
5352  & err,error,*999)
5353  CASE(solver_full_matrices)
5355  & err,error,*999)
5356  CASE DEFAULT
5357  local_error="The specified solver equations sparsity type of "// &
5358  & trim(numbertovstring(solver_equations%SPARSITY_TYPE,"*",err,error))// &
5359  & " is invalid."
5360  CALL flagerror(local_error,err,error,*999)
5361  END SELECT
5362  ENDIF
5363  CALL solver_matrices_create_finish(solver_matrices,err,error,*999)
5364  !Link linear solver
5365  linear_solver%SOLVER_EQUATIONS=>solver%SOLVER_EQUATIONS
5366  !Finish the creation of the linear solver
5367  CALL solver_linear_create_finish(linear_solver%LINEAR_SOLVER,err,error,*999)
5368  ELSE
5369  CALL flagerror("Dynamic solver linear solver is not associated.",err,error,*999)
5370  ENDIF
5371  ELSE IF(dynamic_solver%LINEARITY==solver_dynamic_nonlinear) THEN
5372  nonlinear_solver=>dynamic_solver%NONLINEAR_SOLVER
5373  IF(ASSOCIATED(nonlinear_solver)) THEN
5374  NULLIFY(solver_matrices)
5375  CALL solver_matrices_create_start(solver_equations,solver_matrices,err,error,*999)
5376  CALL solver_matrices_library_type_get(nonlinear_solver,nonlinear_library_type,err,error,*999)
5377  CALL solver_matrices_library_type_set(solver_matrices,nonlinear_library_type,err,error,*999)
5378  IF(dynamic_solver%EXPLICIT) THEN
5380  & err,error,*999)
5381  ELSE
5382  SELECT CASE(solver_equations%SPARSITY_TYPE)
5385  & err,error,*999)
5386  CASE(solver_full_matrices)
5388  & err,error,*999)
5389  CASE DEFAULT
5390  local_error="The specified solver equations sparsity type of "// &
5391  & trim(numbertovstring(solver_equations%SPARSITY_TYPE,"*",err,error))// &
5392  & " is invalid."
5393  CALL flagerror(local_error,err,error,*999)
5394  END SELECT
5395  ENDIF
5396  CALL solver_matrices_create_finish(solver_matrices,err,error,*999)
5397  !Link nonlinear solver
5398  nonlinear_solver%SOLVER_EQUATIONS=>solver%SOLVER_EQUATIONS
5399  !Finish the creation of the nonlinear solver
5400  CALL solver_nonlinear_create_finish(nonlinear_solver%NONLINEAR_SOLVER,err,error,*999)
5401  ELSE
5402  CALL flagerror("Dynamic solver linear solver is not associated.",err,error,*999)
5403  ENDIF
5404  ENDIF
5405  ELSE
5406  CALL flagerror("Solver equations solver mapping is not associated.",err,error,*999)
5407  ENDIF
5408  ELSE
5409  CALL flagerror("Solver solver equations is not associated.",err,error,*999)
5410  ENDIF
5411  CASE(solver_petsc_library)
5412  CALL flagerror("Not implemented.",err,error,*999)
5413  CASE DEFAULT
5414  local_error="The solver library type of "// &
5415  & trim(numbertovstring(dynamic_solver%SOLVER_LIBRARY,"*",err,error))//" is invalid."
5416  CALL flagerror(local_error,err,error,*999)
5417  END SELECT
5418  ELSE
5419  CALL flagerror("Solver solver equations is not associated.",err,error,*999)
5420  ENDIF
5421  ELSE
5422  CALL flagerror("Dynamic solver solver is not associated.",err,error,*999)
5423  ENDIF
5424  ELSE
5425  CALL flagerror("Dynamic solver is not associated.",err,error,*999)
5426  ENDIF
5427 
5428  exits("SOLVER_DYNAMIC_CREATE_FINISH")
5429  RETURN
5430 999 errorsexits("SOLVER_DYNAMIC_CREATE_FINISH",err,error)
5431  RETURN 1
5432 
5433  END SUBROUTINE solver_dynamic_create_finish
5434 
5435  !
5436  !================================================================================================================================
5437  !
5438 
5440  SUBROUTINE solver_dynamic_degree_get(SOLVER,DEGREE,ERR,ERROR,*)
5442  !Argument variables
5443  TYPE(solver_type), POINTER :: SOLVER
5444  INTEGER(INTG), INTENT(OUT) :: DEGREE
5445  INTEGER(INTG), INTENT(OUT) :: ERR
5446  TYPE(varying_string), INTENT(OUT) :: ERROR
5447  !Local Variables
5448  TYPE(dynamic_solver_type), POINTER :: DYNAMIC_SOLVER
5449 
5450  enters("SOLVER_DYNAMIC_DEGREE_GET",err,error,*999)
5451 
5452  IF(ASSOCIATED(solver)) THEN
5453  IF(solver%SOLVER_FINISHED) THEN
5454  IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
5455  dynamic_solver=>solver%DYNAMIC_SOLVER
5456  IF(ASSOCIATED(dynamic_solver)) THEN
5457  degree=dynamic_solver%DEGREE
5458  ELSE
5459  CALL flagerror("Dynamic solver is not associated.",err,error,*999)
5460  ENDIF
5461  ELSE
5462  CALL flagerror("The specified solver is not a dynamic solver.",err,error,*999)
5463  ENDIF
5464  ELSE
5465  CALL flagerror("The solver has not been finished.",err,error,*999)
5466  ENDIF
5467  ELSE
5468  CALL flagerror("Solver is not associated.",err,error,*999)
5469  ENDIF
5470 
5471  exits("SOLVER_DYNAMIC_DEGREE_GET")
5472  RETURN
5473 999 errorsexits("SOLVER_DYNAMIC_DEGREE_GET",err,error)
5474  RETURN 1
5475 
5476  END SUBROUTINE solver_dynamic_degree_get
5477 
5478  !
5479  !================================================================================================================================
5480  !
5481 
5483  SUBROUTINE solver_dynamic_degree_set(SOLVER,DEGREE,ERR,ERROR,*)
5485  !Argument variables
5486  TYPE(solver_type), POINTER :: SOLVER
5487  INTEGER(INTG), INTENT(IN) :: DEGREE
5488  INTEGER(INTG), INTENT(OUT) :: ERR
5489  TYPE(varying_string), INTENT(OUT) :: ERROR
5490  !Local Variables
5491  INTEGER(INTG) :: degree_idx
5492  REAL(DP), ALLOCATABLE :: OLD_THETA(:)
5493  TYPE(dynamic_solver_type), POINTER :: DYNAMIC_SOLVER
5494  TYPE(varying_string) :: LOCAL_ERROR
5495 
5496  enters("SOLVER_DYNAMIC_DEGREE_SET",err,error,*999)
5497 
5498  IF(ASSOCIATED(solver)) THEN
5499  IF(solver%SOLVER_FINISHED) THEN
5500  CALL flagerror("The solver has already been finished.",err,error,*999)
5501  ELSE
5502  IF(solver%SOLVE_TYPE==solver_dynamic_type) THEN
5503  dynamic_solver=>solver%DYNAMIC_SOLVER
5504  IF(ASSOCIATED(dynamic_solver)) THEN
5505  IF(degree/=dynamic_solver%DEGREE) THEN
5506  IF(degree>=dynamic_solver%ORDER) THEN
5507  SELECT CASE(degree)
5509  ALLOCATE(old_theta(dynamic_solver%DEGREE),stat=err)
5510  IF(err/=0) CALL flagerror("Could not allocate old theta.",err,error,*999)
5511  old_theta(1:dynamic_solver%DEGREE)=dynamic_solver%THETA(1:dynamic_solver%DEGREE)
5512  IF(ALLOCATED(dynamic_solver%THETA)) DEALLOCATE(dynamic_solver%THETA)
5513  ALLOCATE(dynamic_solver%THETA(degree),stat=err)
5514  IF(err/=0) CALL flagerror("Could not allocate theta.",err,error,*999)
5515  IF(degree>dynamic_solver%DEGREE) THEN
5516  DO degree_idx=1,dynamic_solver%DEGREE
5517  dynamic_solver%THETA(degree_idx)=old_theta(degree_idx)
5518  ENDDO !degree_idx
5519  DO degree_idx=dynamic_solver%DEGREE+1,degree
5520  dynamic_solver%THETA(degree_idx)=1.0_dp
5521  ENDDO !degree_idx
5522  ELSE
5523  DO degree_idx=1,degree
5524  dynamic_solver%THETA(degree_idx)=old_theta(degree_idx)
5525  ENDDO !degree_idx
5526  ENDIF
5527  IF(ALLOCATED(old_theta)) DEALLOCATE(old_theta)
5528  dynamic_solver%DEGREE=degree
5529  CASE DEFAULT
5530  local_error="The specified degree of "//trim(numbertovstring(degree,"*",err,error))//" is invalid."
5531  CALL flagerror(local_error,err,error,*999)
5532  END SELECT
5533  ELSE
5534  local_error="Invalid dynamic solver setup. The specfied degree of "// &
5535  & trim(numbertovstring(degree,"*",err,error))//" must be >= the current dynamic order of "// &
5536  & trim(numbertovstring(dynamic_solver%ORDER,"*",err,error))//"."
5537  CALL flagerror(local_error,err,error,*999)
5538  ENDIF
5539  ENDIF
5540  ELSE
5541  CALL flagerror("Dynamic solver is not associated.",err,error,*999)
5542  ENDIF
5543  ELSE
5544  CALL flagerror("The specified solver is not a dynamic solver.",err,error,*999)
5545  ENDIF
5546  ENDIF
5547  ELSE
5548  CALL flagerror("Solver is not associated.",err,error,*999)
5549  ENDIF
5550 
5551  exits("SOLVER_DYNAMIC_DEGREE_SET")
5552  RETURN
5553 999 IF(ALLOCATED(old_theta)) DEALLOCATE(old_theta)
5554  errorsexits("SOLVER_DYNAMIC_DEGREE_SET",err,error)
5555  RETURN 1
5556  END SUBROUTINE solver_dynamic_degree_set
5557 
5558  !
5559  !================================================================================================================================
5560  !
5561 
5563  RECURSIVE SUBROUTINE solver_dynamic_finalise(DYNAMIC_SOLVER,ERR,ERROR,*)
5565  !Argument variables
5566  TYPE(dynamic_solver_type), POINTER :: DYNAMIC_SOLVER
5567  INTEGER(INTG), INTENT(OUT) :: ERR
5568  TYPE(varying_string), INTENT(OUT) :: ERROR
5569  !Local Variables
5570 
5571  enters("SOLVER_DYNAMIC_FINALISE",err,error,*999)
5572  IF(ASSOCIATED(dynamic_solver)) THEN
5573  IF(ALLOCATED(dynamic_solver%THETA)) THEN
5574 ! CALL WRITE_STRING_VALUE(GENERAL_OUTPUT_TYPE," Dynamic solver - theta = ",DYNAMIC_SOLVER%THETA(1), &
5575 ! & ERR,ERROR,*999)
5576  DEALLOCATE(dynamic_solver%THETA)
5577  ENDIF
5578  CALL solver_finalise(dynamic_solver%LINEAR_SOLVER,err,error,*999)
5579  CALL solver_finalise(dynamic_solver%NONLINEAR_SOLVER,err,error,*999)
5580  DEALLOCATE(dynamic_solver)
5581  ENDIF
5582 
5583  exits("SOLVER_DYNAMIC_FINALISE")
5584  RETURN
5585 999 errorsexits("SOLVER_DYNAMIC_FINALISE",err,error)
5586  RETURN 1
5587 
5588  END SUBROUTINE solver_dynamic_finalise
5589 
5590  !
5591  !================================================================================================================================
5592  !
5593 
5595  SUBROUTINE solver_dynamic_initialise(SOLVER,ERR,ERROR,*)
5597  !Argument variables
5598  TYPE(solver_type), POINTER :: SOLVER
5599  INTEGER(INTG), INTENT(OUT) :: ERR
5600  TYPE(varying_string), INTENT(OUT) :: ERROR
5601  !Local Variables
5602  TYPE(dynamic_solver_type), POINTER :: DYNAMIC_SOLVER
5603 
5604 
5605  enters(