diff --git a/ExampleMakefile b/ExampleMakefile index fe7eb595..83e28c71 100644 --- a/ExampleMakefile +++ b/ExampleMakefile @@ -54,7 +54,7 @@ ifndef GLOBAL_ROOT $(error GLOBAL_ROOT is not setup) else - GLOBAL_CM_ROOT := $(GLOBAL_ROOT) + GLOBAL_CM_ROOT := $(GLOBAL_ROOT) endif ifndef OPENCMISSEXAMPLES_ROOT @@ -142,11 +142,11 @@ $(EXE_DIR) : $(EXECUTABLE) : $(OBJECTS) $(OPENCMISS_LIBRARY) $(EXE_LINK) -o $@ $(OBJECTS) $(OPENCMISS_LIBRARY) $(ELFLAGS) $(EXTERNAL_LIBRARIES) +# $(EXE_LINK) -o $@ $(OBJECTS) $(OPENCMISS_LIBRARY) $(ELFLAGS) $(EXTERNAL_LIBRARIES) -L/people/vbud003/NVIDIA_GPU_Computing_SDK/C/lib -L/people/vbud003/NVIDIA_GPU_Computing_SDK/C/common/lib/linux -L/people/vbud003/NVIDIA_GPU_Computing_SDK/shared//lib -lcudart -L/people/vbud003/NVIDIA_GPU_Computing_SDK/C/lib -L/people/vbud003/NVIDIA_GPU_Computing_SDK/C/common/lib/linux -L/people/vbud003/NVIDIA_GPU_Computing_SDK/shared//lib -L/usr/local/cuda/lib64/ -lcudart -lcutil_x86_64 -lshrutil_x86_64 $(OPENCMISS_LIBRARY) : ( cd $(GLOBAL_CM_ROOT); $(MAKE) ) - # Place the list of dependencies for the objects here. # # ---------------------------------------------------------------------------- @@ -215,13 +215,14 @@ help: @echo "Options: (The former is the default unless specified.)" @echo @echo " (DEBUG=|OPT=)" - @echo " MPI=(mpich2|intel|openmpi|mvapich2|cray)" + @echo " MPI=(mpich2|intel|openmpi)" @echo " PROF=(true|)" @echo " MPIPROF=(true|)" @echo " ABI=(32|64)" @echo " COMPILER=(intel|gnu|ibm|cray)" @echo " USECELLML=(false|true)" - @echo " USEFIELDML=(false|true)" + @echo " USEFIELDML=(true|false)" + @echo " USECUDA=(false|true)" @echo @echo "Available targets: " @echo diff --git a/Makefile b/Makefile index f1db49f9..addbebaf 100644 --- a/Makefile +++ b/Makefile @@ -102,12 +102,17 @@ main: preliminaries \ PREPROCESSED_OBJECTS = +.SUFFIXES: .f90 .c .cu + $(OBJECT_DIR)/%.o : $(SOURCE_DIR)/%.f90 $(OBJECT_DIR)/.directory ( cd $(OBJECT_DIR) && $(FC) -o $@ $(FFLAGS) $(FPPFLAGS) -c $< ) $(OBJECT_DIR)/%.o : $(SOURCE_DIR)/%.c $(OBJECT_DIR)/.directory ( cd $(OBJECT_DIR) && $(CC) -o $@ $(CFLAGS) $(CPPFLAGS) -c $< ) +$(OBJECT_DIR)/%.o : $(SOURCE_DIR)/%.cu + ( cd $(OBJECT_DIR) ; $(NVCC) -o $@ $(NVCCFLAGS) $(NVCCPPFLAGS) -c $< ) + $(PREPROCESSED_OBJECTS) : $(OBJECT_DIR)/%.o : $(SOURCE_DIR)/%.f90 $(OBJECT_DIR)/.directory ( m4 --prefix-builtins $< > $(subst .o,-expanded.f90,$@) && cd $(OBJECT_DIR) && $(FC) -o $@ $(FFLAGS) $(FPPFLAGS) -c $(subst .o,-expanded.f90,$@) ) @@ -125,6 +130,12 @@ else FIELDML_OBJECT = # endif +ifeq ($(USECUDA),true) + CUDA_OBJECT = $(OBJECT_DIR)/external_dae_solver_routines.o +else + CUDA_OBJECT = $(OBJECT_DIR)/external_dae_solver_routines_dummy.o +endif + #ifeq ($(COMPILER),intel) # TODO: temporarily disable intel build for opencmiss.f90 and etc. # FIELDML_OBJECT = # # MOD_INCLUDE := # @@ -180,7 +191,6 @@ OBJECTS = $(OBJECT_DIR)/advection_diffusion_equation_routines.o \ $(OBJECT_DIR)/equations_matrices_routines.o \ $(OBJECT_DIR)/equations_set_constants.o \ $(OBJECT_DIR)/equations_set_routines.o \ - $(OBJECT_DIR)/external_dae_solver_routines.o \ $(OBJECT_DIR)/field_routines.o \ $(OBJECT_DIR)/field_IO_routines.o \ $(OBJECT_DIR)/finite_elasticity_routines.o \ @@ -233,6 +243,8 @@ OBJECTS = $(OBJECT_DIR)/advection_diffusion_equation_routines.o \ $(OBJECT_DIR)/trees.o \ $(OBJECT_DIR)/types.o \ $(OBJECT_DIR)/util_array.o \ + $(OBJECT_DIR)/vtk_import_routines.o \ + $(CUDA_OBJECT) \ $(FIELDML_OBJECT) \ $(PREPROCESSED_OBJECTS) @@ -444,6 +456,8 @@ $(OBJECT_DIR)/cmiss_cellml_dummy.o : $(SOURCE_DIR)/cmiss_cellml_dummy.f90 $(OBJECT_DIR)/cmiss_fortran_c.o : $(SOURCE_DIR)/cmiss_fortran_c.f90 +$(OBJECT_DIR)/cmiss_fortran_c.o : $(SOURCE_DIR)/cmiss_fortran_c.f90 + $(OBJECT_DIR)/cmiss_mpi.o : $(SOURCE_DIR)/cmiss_mpi.f90 \ $(OBJECT_DIR)/base_routines.o \ $(OBJECT_DIR)/constants.o \ @@ -874,8 +888,10 @@ $(OBJECT_DIR)/equations_set_routines.o : $(SOURCE_DIR)/equations_set_routines.f9 $(OBJECT_DIR)/timer_f.o \ $(OBJECT_DIR)/types.o -$(OBJECT_DIR)/external_dae_solver_routines.o : $(SOURCE_DIR)/external_dae_solver_routines.c \ - $(SOURCE_DIR)/external_dae_solver_routines.h +$(OBJECT_DIR)/external_dae_solver_routines.o : $(SOURCE_DIR)/external_dae_solver_routines.cu \ + $(SOURCE_DIR)/external_dae_solver_routines.h + +$(OBJECT_DIR)/external_dae_solver_routines_dummy.o : $(SOURCE_DIR)/external_dae_solver_routines_dummy.c $(OBJECT_DIR)/field_routines.o : $(SOURCE_DIR)/field_routines.f90 \ $(OBJECT_DIR)/base_routines.o \ @@ -912,7 +928,8 @@ $(OBJECT_DIR)/field_IO_routines.o : $(SOURCE_DIR)/field_IO_routines.f90 \ $(OBJECT_DIR)/mesh_routines.o \ $(OBJECT_DIR)/node_routines.o \ $(OBJECT_DIR)/strings.o \ - $(OBJECT_DIR)/types.o + $(OBJECT_DIR)/types.o \ + $(OBJECT_DIR)/vtk_import_routines.o $(OBJECT_DIR)/finite_elasticity_routines.o : $(SOURCE_DIR)/finite_elasticity_routines.f90 \ $(OBJECT_DIR)/base_routines.o \ @@ -1404,7 +1421,6 @@ $(OBJECT_DIR)/Poiseuille_equations_routines.o : $(SOURCE_DIR)/Poiseuille_equatio $(OBJECT_DIR)/timer_f.o \ $(OBJECT_DIR)/types.o - $(OBJECT_DIR)/Poisson_equations_routines.o : $(SOURCE_DIR)/Poisson_equations_routines.f90 \ $(OBJECT_DIR)/base_routines.o \ $(OBJECT_DIR)/basis_routines.o \ @@ -1511,7 +1527,7 @@ $(OBJECT_DIR)/solver_routines.o : $(SOURCE_DIR)/solver_routines.f90 \ $(OBJECT_DIR)/constants.o \ $(OBJECT_DIR)/distributed_matrix_vector.o \ $(OBJECT_DIR)/equations_set_constants.o \ - $(OBJECT_DIR)/external_dae_solver_routines.o \ + $(CUDA_OBJECT) \ $(OBJECT_DIR)/field_routines.o \ $(OBJECT_DIR)/kinds.o \ $(OBJECT_DIR)/input_output.o \ @@ -1622,6 +1638,9 @@ $(OBJECT_DIR)/util_array.o : $(SOURCE_DIR)/util_array.f90 \ $(OBJECT_DIR)/base_routines.o \ $(OBJECT_DIR)/types.o +$(OBJECT_DIR)/vtk_import_routines.o : $(SOURCE_DIR)/vtk_import_routines.c \ + $(SOURCE_DIR)/vtk_import_routines.h + # ---------------------------------------------------------------------------- # # SWIG bindings to other languages diff --git a/src/analytic_analysis_routines.f90 b/src/analytic_analysis_routines.f90 index d151c506..bcc887a6 100644 --- a/src/analytic_analysis_routines.f90 +++ b/src/analytic_analysis_routines.f90 @@ -161,6 +161,8 @@ SUBROUTINE ANALYTIC_ANALYSIS_OUTPUT(FIELD,FILENAME,ERR,ERROR,*) LOCAL_STRING="Field "//TRIM(NUMBER_TO_VSTRING(FIELD%USER_NUMBER,"*",ERR,ERROR))//" : "//FIELD%LABEL IF(ERR/=0) GOTO 999 CALL WRITE_STRING(OUTPUT_ID,LOCAL_STRING,ERR,ERROR,*999) + NULLIFY(NUMERICAL_VALUES) + NULLIFY(ANALYTIC_VALUES) !Loop over the variables DO var_idx=1,FIELD%NUMBER_OF_VARIABLES variable_type=FIELD%VARIABLES(var_idx)%VARIABLE_TYPE diff --git a/src/cmiss_cellml.f90 b/src/cmiss_cellml.f90 index 49dd5b35..51cf8542 100644 --- a/src/cmiss_cellml.f90 +++ b/src/cmiss_cellml.f90 @@ -438,7 +438,9 @@ SUBROUTINE CELLML_CREATE_FINISH(CELLML,ERR,ERROR,*) !Check that we have set up the models IF(CELLML%NUMBER_OF_MODELS>0) THEN DO model_idx=1,CELLML%NUMBER_OF_MODELS - write(*,*) 'model_idx = ',model_idx + IF(DIAGNOSTICS1) THEN + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'model_idx = ',model_idx,ERR,ERROR,*999) + ENDIF CELLML_MODEL => CELLML%MODELS(model_idx)%PTR !CALL CELLML_MODEL_DEFINITION_SET_SAVE_TEMP_FILES(CELLML_MODEL%PTR,1) ERROR_CODE = CELLML_MODEL_DEFINITION_INSTANTIATE(CELLML_MODEL%PTR) @@ -2734,8 +2736,11 @@ SUBROUTINE CELLML_STATE_FIELD_CREATE_FINISH(CELLML,ERR,ERROR,*) & TRIM(NUMBER_TO_VSTRING(state_component_idx,"*",ERR,ERROR))//"." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF - WRITE(*,*) '(single model) Initial value for state variable: ',state_component_idx,'; type: ',& - & CELLML_VARIABLE_TYPE,'; value = ',INITIAL_VALUE + IF(DIAGNOSTICS1) THEN + CALL WRITE_STRING_TWO_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'(single model) Initial value for state variable: ', & + & state_component_idx,'; type: ',CELLML_VARIABLE_TYPE,ERR,ERROR,*999) + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'value = ',INITIAL_VALUE,ERR,ERROR,*999) + ENDIF CALL FIELD_COMPONENT_VALUES_INITIALISE(CELLML%STATE_FIELD%STATE_FIELD,FIELD_U_VARIABLE_TYPE, & & FIELD_VALUES_SET_TYPE,state_component_idx,INITIAL_VALUE,ERR,ERROR,*999) ENDDO !state_component_idx @@ -2765,8 +2770,11 @@ SUBROUTINE CELLML_STATE_FIELD_CREATE_FINISH(CELLML,ERR,ERROR,*) & TRIM(NUMBER_TO_VSTRING(state_component_idx,"*",ERR,ERROR))//"." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF - WRITE(*,*) '(multiple models) Initial value for state variable: ',state_component_idx,'; type: ',& - & CELLML_VARIABLE_TYPE,'; value = ',INITIAL_VALUE + IF(DIAGNOSTICS1) THEN + CALL WRITE_STRING_TWO_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'(multiple models) Initial value for state variable: '& + &,state_component_idx,'; type: ', CELLML_VARIABLE_TYPE,ERR,ERROR,*999) + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'value = ',INITIAL_VALUE,ERR,ERROR,*999) + ENDIF !\todo make diagnostic output CALL CELLML_FIELD_VARIABLE_SOURCE_DOF_SET_CONSTANT(STATE_VARIABLE,FIELD_VALUES_SET_TYPE,source_dof_idx, & & state_component_idx,INITIAL_VALUE,ERR,ERROR,*999) ENDDO !state_component_idx @@ -3607,8 +3615,11 @@ SUBROUTINE CELLML_PARAMETERS_FIELD_CREATE_FINISH(CELLML,ERR,ERROR,*) & TRIM(NUMBER_TO_VSTRING(parameter_component_idx,"*",ERR,ERROR))//"." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF - WRITE(*,*) '(single model) Initial value for parameter variable: ',parameter_component_idx,'; type: ',& - & CELLML_VARIABLE_TYPE,'; value = ',INITIAL_VALUE + IF(DIAGNOSTICS1) THEN + CALL WRITE_STRING_TWO_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'(single model) Initial value for state variable: ', & + & parameter_component_idx,'; type: ', CELLML_VARIABLE_TYPE,ERR,ERROR,*999) + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'value = ',INITIAL_VALUE,ERR,ERROR,*999) + ENDIF CALL FIELD_COMPONENT_VALUES_INITIALISE(CELLML%PARAMETERS_FIELD%PARAMETERS_FIELD,FIELD_U_VARIABLE_TYPE, & & FIELD_VALUES_SET_TYPE,parameter_component_idx,INITIAL_VALUE,ERR,ERROR,*999) ENDDO !parameter_component_idx @@ -3638,8 +3649,11 @@ SUBROUTINE CELLML_PARAMETERS_FIELD_CREATE_FINISH(CELLML,ERR,ERROR,*) & TRIM(NUMBER_TO_VSTRING(parameter_component_idx,"*",ERR,ERROR))//"." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF - WRITE(*,*) '(multiple models) Initial value for parameter variable: ',parameter_component_idx,'; type: ',& - & CELLML_VARIABLE_TYPE,'; value = ',INITIAL_VALUE + IF(DIAGNOSTICS1) THEN + CALL WRITE_STRING_TWO_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'(multiple models) Initial value for state variable: '& + &, parameter_component_idx,'; type: ', CELLML_VARIABLE_TYPE,ERR,ERROR,*999) + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,'value = ',INITIAL_VALUE,ERR,ERROR,*999) + ENDIF !\todo make diagnostic output CALL CELLML_FIELD_VARIABLE_SOURCE_DOF_SET_CONSTANT(PARAMETERS_VARIABLE,FIELD_VALUES_SET_TYPE, & & source_dof_idx, parameter_component_idx,INITIAL_VALUE,ERR,ERROR,*999) ENDDO !parameter_component_idx diff --git a/src/cuda_solver_routines.cu b/src/cuda_solver_routines.cu new file mode 100644 index 00000000..5d174a90 --- /dev/null +++ b/src/cuda_solver_routines.cu @@ -0,0 +1,721 @@ +/* \file + * \author Jin Budlemann + * \brief This file provides the routines for solving differential-algebraic equations with an external solver. + *. + * \section LICENSE + * + * Version: MPL 1.1/GPL 2.0/LGPL 2.1 + * + * The contents of this file are subject to the Mozilla Public License + * Version 1.1 (the "License"); you may not use this file except in + * compliance with the License. You may obtain a copy of the License at + * http://www.mozilla.org/MPL/ + * + * Software distributed under the License is distributed on an "AS IS" + * basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the + * License for the specific language governing rights and limitations + * under the License. + * + * The Original Code is OpenCMISS + * + * The Initial Developer of the Original Code is University of Auckland, + * Auckland, New Zealand, the University of Oxford, Oxford, United + * Kingdom and King's College, London, United Kingdom. Portions created + * by the University of Auckland, the University of Oxford and King's + * College, London are Copyright (C) 2007-2010 by the University of + * Auckland, the University of Oxford and King's College, London. + * All Rights Reserved. + * + * Contributor(s): Chris Bradley + * + * Alternatively, the contents of this file may be used under the terms of + * either the GNU General Public License Version 2 or later (the "GPL"), or + * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), + * in which case the provisions of the GPL or the LGPL are applicable instead + * of those above. If you wish to allow use of your version of this file only + * under the terms of either the GPL or the LGPL, and not to allow others to + * use your version of this file under the terms of the MPL, indicate your + * decision by deleting the provisions above and replace them with the notice + * and other provisions required by the GPL or the LGPL. If you do not delete + * the provisions above, a recipient may use your version of this file under + * the terms of any one of the MPL, the GPL or the LGPL. + * + */ + +/* +File: cuda_solver_routines.cu +============================= + +This file provides provides the routines for cuda acceleration with an external solver. + +Functions included: + + +*/ + +#include +#include + +#include + +//extern __shared__ double shared_array[]; +//const int algebraicCount = 25; +//const int rateStateCount = 8; +//const int constantCount = 0; +//const double FLOPSPerFunction = 193.0f; +//const int DEFAULT_TESTING_THREADS = 2000000; +//const int sharedMemoryCellModel = 0; +//const char* cellModelName = "LR R3"; +// +////////////////////////////////////////////////////////////////////////////////// +//// Cell Model Device Functions +////////////////////////////////////////////////////////////////////////////////// +//__device__ void computeRates(float VOI, double* DUMMY, double* STATES, double* ALGEBRAIC, double* RATES) +//{ +// ALGEBRAIC[0] = -25.5; // Add stimulus in proper +// ALGEBRAIC[1] = (0.32f*STATES[0]+15.0816f)/(1.0f - (expf(- 0.1f*STATES[0]-4.713f))); // 7 ops +// +// if (STATES[0] < -40.0f) { +// ALGEBRAIC[2] = 0.135f*(expf(((80.0f+STATES[0])/- 6.8f))); // 4 ops +// ALGEBRAIC[3] = (( - 127140*(expf(0.2444f*STATES[0])) - 3.47400e-05*(expf(-0.04391f*STATES[0])))*(STATES[0]+37.78f))/(1.0f+(expf(0.311f*STATES[0]+24.64053))); // 14 ops +// ALGEBRAIC[9] = 3.56f*(expf(0.079f*STATES[0]))+ 310000*(expf(0.35f*STATES[0])); // 7 ops +// ALGEBRAIC[10] = 0.1212f*(expf(-0.01052f*STATES[0]))/(1.0f+(expf(-0.1378f*STATES[0]-5.531292f))); // 8 ops +// } else { +// ALGEBRAIC[2] = 0.00000; +// ALGEBRAIC[3] = 0.00000; +// ALGEBRAIC[9] = 1.0f/( 0.13f*(1.0f+(expf(((STATES[0]+10.66f)/- 11.1f))))); +// ALGEBRAIC[10] = ( 0.3f*(expf(-2.53500e-07*STATES[0])))/(1.0f+(expf(-0.1f*STATES[0]-3.2f))); +// } +// if (STATES[0] < -100.0f) { +// ALGEBRAIC[16] = 2.837f*(expf(0.04f*STATES[0]+3.08f) - 1.0f)/((STATES[0]+77.0f)*(expf(0.04f*STATES[0]+1.4f))); // 11 ops +// } else { +// ALGEBRAIC[16] = 1.0f; +// } +// +// ALGEBRAIC[4] = (0.095f*(expf(-0.01f*STATES[0] + 0.5f)))/(1.0f+(expf(-0.072*STATES[0] + 0.36f))); // 9 ops +// ALGEBRAIC[5] = (0.012f*(expf(-0.008f*STATES[0]-0.224f)))/(1.0f+(expf(0.15f*STATES[0]+4.2f))); // 9 ops +// ALGEBRAIC[6] = (0.0005f*(expf(0.083f*STATES[0]+4.15f)))/(1.0f+(expf(0.057f*STATES[0]+2.85f))); // 9 ops +// ALGEBRAIC[7] = 23*(powf(STATES[1], 3.0f))*STATES[2]*STATES[3]*(STATES[0] - 54.794463f); // 6 ops +// ALGEBRAIC[8] = 0.08f*(expf(-STATES[0]/11.0000)); // 3 ops +// ALGEBRAIC[11] = (0.07f*(expf(-0.017f*STATES[0]-0.748f)))/(1.0f+(expf(0.05f*STATES[0]+2.2f))); // 9 ops +// ALGEBRAIC[12] = (0.0065f*(expf(-0.02f*STATES[0]-0.6f)))/(1.0f+(expf(-0.2f*STATES[0]-6.0f))); // 9 ops +// ALGEBRAIC[13] = (0.0013f*(expf(-0.06f*STATES[0]-1.2f)))/(1.0f+(expf(-0.04f*STATES[0]-0.8f))); // 9 ops +// ALGEBRAIC[14] = 7.7f - 13.0287f*logf(STATES[4]); // 3 ops +// ALGEBRAIC[15] = 0.09f*STATES[5]*STATES[6]*(STATES[0] - ALGEBRAIC[14]); // 4 ops +// ALGEBRAIC[17] = 0.282f*STATES[7]*ALGEBRAIC[16]*(STATES[0] + 77.56758f); // 4 ops +// ALGEBRAIC[18] = 1.02f/(1.0f+(expf(0.2385f*STATES[0] + 6.83967915f))); // 4 ops +// ALGEBRAIC[19] = (0.49124f*(expf( 0.08032f *STATES[0] + 7.49939f) + expf(0.06175f*STATES[0] - 31.271255925f)))/(1.00000+expf(-0.514300*STATES[0] - 214.85137268791f)); // 13 ops +// ALGEBRAIC[20] = ALGEBRAIC[18]/(ALGEBRAIC[18]+ALGEBRAIC[19]); // 2 ops +// ALGEBRAIC[21] = 0.6047f*ALGEBRAIC[20]*(STATES[0] + 87.89290f); // 3 ops +// ALGEBRAIC[22] = 1.0f/(1.0f+(exp(((7.488f - STATES[0])/5.98f)))); // 5 ops +// ALGEBRAIC[23] = 0.0183f*ALGEBRAIC[22]*(STATES[0] + 87.89290f); // 3 ops +// ALGEBRAIC[24] = 0.03921f*STATES[0] +2.3475027f; // 3 ops +// +// RATES[0] = -(ALGEBRAIC[0]+ALGEBRAIC[7]+ALGEBRAIC[15]+ALGEBRAIC[17]+ALGEBRAIC[21]+ALGEBRAIC[23]+ALGEBRAIC[24]); // 7 ops +// RATES[1] = ALGEBRAIC[1]*(1.00000 - STATES[1]) - ALGEBRAIC[8]*STATES[1]; // 4 ops +// RATES[2] = ALGEBRAIC[2]*(1.00000 - STATES[2]) - ALGEBRAIC[9]*STATES[2]; // 4 ops +// RATES[3] = ALGEBRAIC[3]*(1.00000 - STATES[3]) - ALGEBRAIC[10]*STATES[3]; // 4 ops +// RATES[4] = - 0.0001f*ALGEBRAIC[15]+ 0.000007f - 0.07f*STATES[4]; // 4 ops +// RATES[5] = ALGEBRAIC[4]*(1.00000 - STATES[5]) - ALGEBRAIC[11]*STATES[5]; // 4 ops +// RATES[6] = ALGEBRAIC[5]*(1.00000 - STATES[6]) - ALGEBRAIC[12]*STATES[6]; // 4 ops +// RATES[7] = ALGEBRAIC[6]*(1.00000 - STATES[7]) - ALGEBRAIC[13]*STATES[7]; // 4 ops +//} + +extern __shared__ double shared_array[]; +const int rateStateCount = 2; +const int constantCount = 0; +const int algebraicCount = 1; +const double FLOPSPerFunction = 9.0f; +const int DEFAULT_TESTING_THREADS = 20000000; +const int sharedMemoryCellModel = 0; +const char* cellModelName = "FN R1"; + +//////////////////////////////////////////////////////////////////////////////// +// Cell Model Device Functions +//////////////////////////////////////////////////////////////////////////////// +__device__ void computeRates(double time, double* constants, double* states, double* algebraic, double* rates) +{ + rates[1] = 0.005f*(states[0] - 3.0f*states[1]); + rates[0] = ((states[0]*(states[0] - -0.08f)*(1.0f - states[0]) - states[1])+ algebraic[0]); +} + +void initProblem(int num_threads, double* STATES) { + printf("test\n\n\n\n\n\n"); +} + +////////////////////////////////////////////////////////////////////////////////// +//// Cell Model Host Functions ////////// Should Not Be Needed Later ///////////// +////////////////////////////////////////////////////////////////////////////////// +//void initProblem(int num_threads, double* STATES) +//{ +// int i; +// +// STATES[0] = -84.3801107371; +// STATES[1] = 0.00171338077730188; +// STATES[2] = 0.982660523699656; +// STATES[3] = 0.989108212766685; +// STATES[4] = 0.00017948816388306; +// STATES[5] = 0.00302126301779861; +// STATES[6] = 0.999967936476325; +// STATES[7] = 0.0417603108167287; +// +// for (i=1; i 0 + double threadConstants[constantCount]; + + intitialiseConstants(threadConstants); +#endif + + for (i=0; i 0 + integrator(timeSteps, stepSize, threadConstants, threadStates + threadIdx.x*rateStateCount, threadAlgebraic); +#else + integrator(timeSteps, stepSize, NULL, threadStates + threadIdx.x*rateStateCount, threadAlgebraic); +#endif + + for (i=0; i deviceProp.maxThreadsPerBlock ? deviceProp.maxThreadsPerBlock : (*threads_per_block); + + // Adjust threads per blocks so that states variables fit in shared memory + count=0; + (*sharedMem) = (size_t) (sharedMemoryIntegrator + sharedMemoryDevice + sharedMemoryCellModel) * (*threads_per_block) * sizeof(double); + while ( (*sharedMem) > deviceProp.sharedMemPerBlock*0.5 && count < 200 ) { + if ((*threads_per_block) > 32 ) { + (*threads_per_block)-=32; + } else { + fprintf(stderr, "Cannot fit variables in shared memory"); + exit(EXIT_FAILURE); + } + count++; + (*sharedMem) = (size_t) (sharedMemoryIntegrator + sharedMemoryDevice + sharedMemoryCellModel) * (*threads_per_block) * sizeof(double); + } + (*num_blocks)=num_threads/(*threads_per_block); + spill[2]=num_threads % (*threads_per_block); + //if (spill[2] !=0) (*num_blocks)++; // Round up calculation + + (*blocksPerDomain) = (*num_blocks)/(num_streams*(*num_partitions)); + if ((*blocksPerDomain) == 0) { + fprintf(stderr, "Too many streams and partitions for your problem size\nReduce these and try again"); + exit(EXIT_FAILURE); + } else if ((*blocksPerDomain) > deviceProp.maxGridSize[0]) { + fprintf(stderr, "Too many blocks allocated to each domain for your problem size\nIncrease number of partitions and try again"); + exit(EXIT_FAILURE); + } + remainingBlocks = (*num_blocks)%(num_streams*(*num_partitions)); + if (remainingBlocks > 0 && remainingBlocks <= num_streams*(*blocksPerDomain)) { + (*num_partitions)++; + } else if (remainingBlocks > num_streams*(*blocksPerDomain)) { + numberRemainingDomains = ceil((float)remainingBlocks/(*blocksPerDomain)); + (*num_partitions) += ceil((float)numberRemainingDomains/num_streams); + } + + if ((*num_blocks)/(num_streams*(*num_partitions))==0) { + (*num_partitions)--; + } + + remainingBlocks = (*num_blocks)%(*num_partitions); + + // Adjust number of partitions so that data fits in GPU memory + count = 0; + cutilSafeCall( cudaMemGetInfo(&freeGPUMemory, &test) ); + *pagedMemorySize = sizeof(double)*rateStateCount*(*blocksPerDomain)*num_streams*(*threads_per_block); + while ( *pagedMemorySize > freeGPUMemory*0.9 && count < 200 ) { + if ((*blocksPerDomain) > 1 ) { + (*num_partitions) = (*num_blocks)/num_streams/--(*blocksPerDomain); + remainingBlocks = (*num_blocks)%(num_streams*(*num_partitions)); + if (remainingBlocks != 0) { + (*num_partitions)++; + } else if (remainingBlocks > num_streams*(*blocksPerDomain)) { + numberRemainingDomains = ceil((float)remainingBlocks/(*blocksPerDomain)); + (*num_partitions) += ceil((float)numberRemainingDomains/num_streams); + } + remainingBlocks = (*num_blocks)%(*num_partitions); + } else { + fprintf(stderr, "Cannot fit variables in GPU device memory\nReduce threads per block and try again"); + exit(EXIT_FAILURE); + } + count++; + + *pagedMemorySize = sizeof(double)*rateStateCount*(*blocksPerDomain)*num_streams*(*threads_per_block); + } + + if (*pagedMemorySize > freeGPUMemory*0.9) { + fprintf(stderr, "Cannot fit variables in GPU device memory\nReduce threads per block and try again"); + exit(EXIT_FAILURE); + } + + (*threads_per_domain) = (*threads_per_block)*(*blocksPerDomain); + spill[0] = remainingBlocks/(*blocksPerDomain); + spill[1] = remainingBlocks%(*blocksPerDomain); + + // if (renainingBlocks != 0) { + // blocksLastStream = + // } + // remainingBlocks = num_threads%(num_blocks*(*num_partitions)) + // for (i=0; i<(*num_partitions); i++) { + // if (i == (*num_partitions) && remainingBlocks != 0) { + // + // for (j=0; j Device name : %s\n", deviceProp.name ); + // printf("> CUDA Capable SM %d.%d hardware with %d multi-processors\n", + // deviceProp.major, deviceProp.minor, deviceProp.multiProcessorCount); + // printf("> Cell Model = %s, Integrator = %s\n", cellModelName, integratorName); + // printf("> num_threads = %d, num_blocks = %d, threads_per_block = %d, num_partitions = %d, timeSteps = %d, num_streams = %d\n", + // num_threads, num_blocks, threads_per_block, num_partitions, timeSteps, num_streams); + // //#endif + // printf("grid.x %d threads.x %d sharedMem %d\n", grid.x, threads.x, sharedMem); + // printf("Spills %d %d %d\n", spill[0], spill[1], spill[2]); + + // Setup and start global timer + timer = 0; + cutCreateTimer(&timer); + cutilCheckError(cutStartTimer(timer)); + + // // Test kernel speed in default stream (timing is more accurate in default stream) + // memcpy(h_paged_states, h_states, pagedMemorySize/num_streams); + // cutilSafeCall( cudaMemcpy(d_states, h_paged_states, pagedMemorySize/num_streams, + // cudaMemcpyHostToDevice) ); + // // Start kernel timer + // kernel_timer = 0; + // cutCreateTimer(&kernel_timer); + // cutilCheckError(cutStartTimer(kernel_timer)); + // // Start kernel + // solveSystem<<>>(timeSteps, stepSize, d_states); + // checkCUDAError("Single Kernel Execution"); + // cutilSafeCall( cudaThreadSynchronize() ); + // // Stop kernel Timer + // cutilCheckError(cutStopTimer(kernel_timer)); + // cutilSafeCall( cudaMemcpy(h_paged_states, d_states, pagedMemorySize/num_streams, + // cudaMemcpyDeviceToHost) ); + // memcpy(h_states, h_paged_states, pagedMemorySize/num_streams); + + // Prefetch data for next partition in first stream + // if (num_partitions>1) { + // global_offset = rateStateCount * num_streams * grid.x * threads.x; + // memcpy(h_paged_states, h_states + global_offset, pagedMemorySize/num_streams); + // } + // } else { + memcpy(h_paged_states, h_states, pagedMemorySize); + } + + // Queue kernel calls into streams to hide memory transfers (num_partitions sets of kernel calls in each stream) + for(i = 0; i < num_partitions+1; i++) { + // Asynchronously launch num_streams memcopies + for(j = 0; j < num_streams; j++) { //((timing_file!=NULL && i==0) ? 1 : 0); j < num_streams; j++){ + domainIndex = j + i*num_streams; + if (domainIndex <= lastFullDomain) { + local_offset = j * rateStateCount * grid.x * threads.x ; + if (domainIndex == lastFullDomain && (spill[1]!=0 || spill[2]!=0)) { + //printf("last async in %d, size %d\n", domainIndex, lastStreamMemorySize); + cutilSafeCall( cudaMemcpyAsync(d_states + local_offset, h_paged_states + local_offset, + lastStreamMemorySize, cudaMemcpyHostToDevice, streams[j]) ); + } else { + //printf("normal async in %d, size %d\n", domainIndex, pagedMemorySize/num_streams); + cutilSafeCall( cudaMemcpyAsync(d_states + local_offset, h_paged_states + local_offset, + pagedMemorySize/num_streams, cudaMemcpyHostToDevice, streams[j]) ); + } + } + } + // Execute the kernel + // Asynchronously launch num_streams kernels, each operating on its own portion of data + for(j = 0; j < num_streams; j++) { //((timing_file!=NULL && i==0) ? 1 : 0); j < num_streams; j++){ + domainIndex = j + i*num_streams; + if (domainIndex <= lastFullDomain) { + local_offset = j * grid.x * threads.x ; + if (domainIndex == lastFullDomain && (spill[1]!=0 || spill[2]!=0)) { + grid.x = spill[1]+1; + solveSystem<<>>(timeSteps, stepSize, + d_states + rateStateCount*local_offset); + } else { + solveSystem<<>>(timeSteps, stepSize, + d_states + rateStateCount*local_offset); + } + } + } + + // Asynchronoously launch num_streams memcopies + for(j = 0; j < num_streams; j++) { //((timing_file!=NULL && i==0) ? 1 : 0); j < num_streams; j++){ + domainIndex = j + i*num_streams; + if (domainIndex <= lastFullDomain) { + local_offset = j * rateStateCount * grid.x * threads.x ; + if (domainIndex == lastFullDomain && (spill[1]!=0 || spill[2]!=0)) { + //printf("last async out %d, size %d\n", domainIndex, lastStreamMemorySize); + cutilSafeCall( cudaMemcpyAsync(h_paged_states + local_offset, d_states + local_offset, + lastStreamMemorySize, cudaMemcpyDeviceToHost, streams[j]) ); + } else { + //printf("normal async out %d, size %d\n", domainIndex, pagedMemorySize/num_streams); + cutilSafeCall( cudaMemcpyAsync(h_paged_states + local_offset, d_states + local_offset, + pagedMemorySize/num_streams, cudaMemcpyDeviceToHost, streams[j]) ); + } + } + } + + // Execute memcpys in and out of paged memory when CUDA calls in the streams have finished + for(j = 0; j < num_streams; j++) { //((timing_file!=NULL && i==0) ? 1 : 0); j < num_streams; j++){ + domainIndex = j + i*num_streams; + if (domainIndex <= lastFullDomain) { + cudaStreamSynchronize(streams[j]); + + local_offset = j * rateStateCount * grid.x * threads.x ; + global_offset = i * num_streams * grid.x * threads.x; + + if (domainIndex == lastFullDomain && (spill[1]!=0 || spill[2]!=0)) { + //printf("last memcpy out %d\n", domainIndex); + memcpy(h_states + rateStateCount * global_offset + local_offset, h_paged_states + local_offset, + lastStreamMemorySize); + } else { + //printf("normal memcpy out %d\n", domainIndex); + memcpy(h_states + rateStateCount * global_offset + local_offset, h_paged_states + local_offset, + pagedMemorySize/num_streams); + } + + global_offset = (i + 1) * num_streams * grid.x * threads.x; + if (domainIndex == lastFullDomain - num_streams && (spill[1]!=0 || spill[2]!=0)) { + //printf("last memcpy in %d\n", domainIndex); + lastStreamMemorySize = sizeof(double)*rateStateCount*(spill[1]*threads_per_block+spill[2]); + memcpy(h_paged_states + local_offset, h_states + rateStateCount * global_offset + local_offset, + lastStreamMemorySize); + } else if (domainIndex < lastFullDomain - num_streams) { + //printf("normal memcpy in %d\n", domainIndex); + memcpy(h_paged_states + local_offset, h_states + rateStateCount * global_offset + local_offset, + pagedMemorySize/num_streams); + } + } + } + } + + if (timing_file) { + // Stop global timer + cutilCheckError(cutStopTimer(timer)); + + // Calculate timing statistics + dSeconds = cutGetTimerValue(timer)/1000.0; + nFLOPS = (FLOPSPerTimeStep*rateStateCount + FLOPSPerFunction*FunctionEvals)*timeSteps*num_threads; + gflops = 1.0e-9 * nFLOPS/dSeconds; + + // kernel_dSeconds = cutGetTimerValue(kernel_timer)/1000.0; + // kernel_nFLOPS = (FLOPSPerTimeStep*rateStateCount + FLOPSPerFunction*FunctionEvals)*timeSteps*num_threads/num_streams/num_partitions; + // kernel_gflops = 1.0e-9 * kernel_nFLOPS/kernel_dSeconds; + + // Store Stats + fprintf(timing_file,"%s\t%s\t%d\t%d\t%d\t%d\t%d\t%f\t%f\n", cellModelName, integratorName, num_threads, num_blocks, + threads_per_block, num_partitions, num_streams, dSeconds, gflops); + // fprintf(timing_file,"%s\t%s\t%d\t%d\t%d\t%d\t%d\t%f\t%f\t%f\t%f\t%f\n", cellModelName, integratorName, num_threads, num_blocks, + // threads_per_block, num_partitions, num_streams, dSeconds, gflops, kernel_dSeconds, + // kernel_gflops, gflops/kernel_gflops*100); + } + + // Deallocate Memory and Release Threads + for(i = 0; i < num_streams; i++) + cutilSafeCall( cudaStreamDestroy(streams[i]) ); + cutilSafeCall( cudaFree(d_states) ); + cutilSafeCall( cudaFreeHost(h_paged_states) ); + cudaThreadExit(); +} + + +//////////////////////////////////////////////////////////////////////////////// +// Auxilliary Testing Functions +//////////////////////////////////////////////////////////////////////////////// + +void solveProblem(unsigned int timeSteps, unsigned int num_threads, unsigned int threads_per_block, + unsigned int num_partitions, unsigned int num_streams, FILE *timing_file, FILE *results_file) +{ + unsigned int i, j; + float startTime = 0.0f; + float endTime = 0.2f; + + double* h_states = NULL; + + h_states = (double *) safeMalloc(sizeof(double)*rateStateCount*num_threads); + + initProblem(num_threads, h_states); + + solve(h_states, startTime, endTime, timeSteps, num_threads, threads_per_block, num_partitions, num_streams, timing_file); + + if (results_file) { + fprintf(results_file,"\n\n"); + + for (i=0; i +#include +#include + + +#include "external_dae_solver_routines.h" +#include "cuda_solver_routines.cu" + +/* Type definitions */ + +/* Function Definitions */ +extern "C" +{ +void SolverDAEExternalIntegrate(const int NumberOfDofs, + const double StartTime, + const double EndTime, + double *InitialStep, + const int ThreadsPerBlock, + const int NumberOfPartitions, + const int NumberOfStreams, + const int OnlyOneModelIndex, + int *ModelsData, + int NumberOfState, + double *StateData, + int NumberOfParameters, + double *ParametersData, + int NumberOfIntermediate, + double *IntermediateData, + int *err) +{ + FILE* timing_file = NULL; + char *filename = NULL; + + asprintf(&filename,"Results/MonodomainExample-CUDAON-%d-%d-%d-%d.txt",(NumberOfDofs/101) - 1,ThreadsPerBlock,NumberOfPartitions,NumberOfStreams); + timing_file = fopen(filename, "rt"); + + if (!timing_file) { + timing_file = fopen(filename, "wt"); + if (!timing_file) { + fprintf(stderr,"File error: %s\n",strerror(errno)); + fprintf(stderr, "Timing file could not be opened or created."); + exit(EXIT_FAILURE); + } + fprintf(timing_file,"Cell Model\tIntegrator\tNumber of Threads\tNumber 0f Blocks\tThreads Per Block\tNumber of Partitions\tNumber of Streams\tTotal Computational Time(s)\tTotal GFLOPS\n"); +// fprintf(timing_file,"Cell Model\tIntegrator\tNumber of Threads\tNumber 0f Blocks\tThreads Per Block\tNumber of Partitions\tNumber of Streams\tTotal Computational Time(s)\tTotal GFLOPS\tSingle Kernel Computaional Time(s)\tKernel GFLOPS\tDevice Utilisation\n"); + } else { + fclose(timing_file); + timing_file = fopen(filename, "at"); + if (!timing_file) { + fprintf(stderr,"File error: %s\n",strerror(errno)); + fprintf(stderr, "Timing file could not be opened or created."); + exit(EXIT_FAILURE); + } + } + + //printf("start %f end %f steps %d\n", StartTime, EndTime, (int)((EndTime-StartTime)/InitialStep[0])); + // timeSteps = (int)ceil(((EndTime-StartTime)/InitialStep[0])); + + solve(StateData, StartTime, EndTime, InitialStep[0], NumberOfDofs, ThreadsPerBlock, NumberOfPartitions, NumberOfStreams, timing_file); + + if (timing_file != NULL ) fclose(timing_file); + free(filename); +} +} + diff --git a/src/external_dae_solver_routines.h b/src/external_dae_solver_routines.h index 3137484d..574b5aaf 100644 --- a/src/external_dae_solver_routines.h +++ b/src/external_dae_solver_routines.h @@ -42,3 +42,22 @@ * */ +extern "C" +{ + void SolverDAEExternalIntegrate(const int NumberOfDofs, + const double StartTime, + const double EndTime, + double *InitialStep, + const int ThreadsPerBlock, + const int NumberOfPartitions, + const int NumberOfStreams, + const int OnlyOneModelIndex, + int *ModelsData, + int NumberOfState, + double *StateData, + int NumberOfParameters, + double *ParametersData, + int NumberOfIntermediate, + double *IntermediateData, + int *err); +} diff --git a/src/external_dae_solver_routines.c b/src/external_dae_solver_routines_dummy.c similarity index 85% rename from src/external_dae_solver_routines.c rename to src/external_dae_solver_routines_dummy.c index 199fc82b..28238937 100644 --- a/src/external_dae_solver_routines.c +++ b/src/external_dae_solver_routines_dummy.c @@ -1,6 +1,6 @@ /* \file * \author Chris Bradley - * \brief This file provides the routines for solving differential-algebraic equations with an external solver. + * \brief This file provides the dummy routines for solving differential-algebraic equations with an external solver. *. * \section LICENSE * @@ -43,10 +43,10 @@ */ /* -File: external_dae_solver_routines.c -=================== +File: external_dae_solver_routines_dummy.c +========================================== -This file provides provides the routines for solving differential-algebraic equations with an external solver. +This file provides provides the dummy routines for solving differential-algebraic equations with an external solver. Functions included: @@ -56,11 +56,6 @@ SolverDAEExternalIntegrate Solves the differential-algebraic equation. /* Included files */ -#include -#include - -#include "external_dae_solver_routines.h" - /* Type definitions */ /* Function Definitions */ @@ -69,6 +64,9 @@ void SolverDAEExternalIntegrate(const int NumberOfDofs, const double StartTime, const double EndTime, double *InitialStep, + const int ThreadsPerBlock, + const int NumberOfPartitions, + const int NumberOfStreams, const int OnlyOneModelIndex, int *ModelsData, int NumberOfState, @@ -80,8 +78,7 @@ void SolverDAEExternalIntegrate(const int NumberOfDofs, int *err) { - printf("Hello World!\n"); - -} + // Do Nothing +} diff --git a/src/field_IO_routines.f90 b/src/field_IO_routines.f90 index ab977e7c..13657a37 100644 --- a/src/field_IO_routines.f90 +++ b/src/field_IO_routines.f90 @@ -352,6 +352,19 @@ FUNCTION FieldExport_DerivativeIndices( handle, componentNumber, fieldType, vari INTEGER(C_INT), VALUE :: valueIndex INTEGER(C_INT) :: FieldExport_DerivativeIndices END FUNCTION FieldExport_DerivativeIndices + + SUBROUTINE READ_VTK(filePath, numberOfMeshDimensions, nodesPerElement, points, numPoints, cells, numCells) & + & BIND(C, NAME="readVTK") + USE TYPES + USE ISO_C_BINDING + CHARACTER(LEN=1,KIND=C_CHAR), INTENT(IN) :: filePath(*) + INTEGER(C_INT), VALUE, INTENT(IN) :: numberOfMeshDimensions + INTEGER(C_INT), VALUE, INTENT(IN) :: nodesPerElement + TYPE(C_PTR), INTENT(OUT):: points + INTEGER(C_INT), INTENT(OUT) :: numPoints + TYPE(C_PTR), INTENT(OUT) :: cells + INTEGER(C_INT), INTENT(OUT) :: numCells + END SUBROUTINE READ_VTK END INTERFACE @@ -382,7 +395,7 @@ END FUNCTION FieldExport_DerivativeIndices MODULE PROCEDURE CHECKED_DEALLOCATE_FIELD END INTERFACE !CHECKED_DEALLOCATE - PUBLIC :: FIELD_IO_FIELDS_IMPORT, FIELD_IO_NODES_EXPORT, FIELD_IO_ELEMENTS_EXPORT + PUBLIC :: FIELD_IO_FIELDS_IMPORT, FIELD_IO_NODES_EXPORT, FIELD_IO_ELEMENTS_EXPORT, READ_VTK CONTAINS @@ -529,7 +542,6 @@ END SUBROUTINE REALLOCATE_BASIS ! !================================================================================================================================ ! - SUBROUTINE REALLOCATE_FIELD( array, newSize, errorMessage, ERR, ERROR, * ) TYPE(FIELD_PTR_TYPE), ALLOCATABLE, INTENT(INOUT) :: array(:) INTEGER(INTG), INTENT(IN) :: newSize @@ -821,6 +833,7 @@ SUBROUTINE CHECKED_DEALLOCATE_BASIS( array ) END SUBROUTINE CHECKED_DEALLOCATE_BASIS + ! !================================================================================================================================ ! diff --git a/src/fieldml_input_routines.f90 b/src/fieldml_input_routines.f90 index 674588ec..548cd0a7 100644 --- a/src/fieldml_input_routines.f90 +++ b/src/fieldml_input_routines.f90 @@ -73,7 +73,7 @@ MODULE FIELDML_INPUT_ROUTINES PUBLIC :: FIELDML_INPUT_INITIALISE_FROM_FILE, FIELDML_INPUT_MESH_CREATE_START, & & FIELDML_INPUT_COORDINATE_SYSTEM_CREATE_START, FIELDML_INPUT_BASIS_CREATE_START, FIELDML_INPUT_CREATE_MESH_COMPONENT, & - & FIELDML_INPUT_FIELD_CREATE_START, FIELDML_INPUT_FIELD_NODAL_PARAMETERS_UPDATE, FIELDML_INPUT_NODES_CREATE_START + & FIELDML_INPUT_FIELD_CREATE_START, FIELDML_INPUT_FIELD_PARAMETERS_UPDATE, FIELDML_INPUT_NODES_CREATE_START CONTAINS @@ -1139,6 +1139,68 @@ END SUBROUTINE FIELDML_INPUT_FIELD_CREATE_START !================================================================================================================================ ! + !>Inputs from a FieldML file the parameters for a field variable parameter set. + SUBROUTINE FIELDML_INPUT_FIELD_PARAMETERS_UPDATE(FIELDML_INFO, EVALUATOR_NAME, FIELD, VARIABLE_TYPE, SET_TYPE, & + & ERR, ERROR, * ) + !Argument variables + TYPE(FIELDML_IO_TYPE), INTENT(INOUT) :: FIELDML_INFO !1) THEN + CALL FIELD_COMPONENT_INTERPOLATION_GET(FIELD,VARIABLE_TYPE,1,INTERPOLATION_TYPE,ERR,ERROR,*999) + CALL FIELD_COMPONENT_MESH_COMPONENT_GET(FIELD,VARIABLE_TYPE,1,MESH_COMPONENT1,ERR,ERROR,*999) + IS_ALL_NODAL_INTERPOLATION=INTERPOLATION_TYPE==FIELD_NODE_BASED_INTERPOLATION + IS_SAME_MESH_COMPONENTS=.TRUE. + DO component_idx=2,NUMBER_OF_COMPONENTS + CALL FIELD_COMPONENT_INTERPOLATION_GET(FIELD,VARIABLE_TYPE,component_idx,INTERPOLATION_TYPE,ERR,ERROR,*999) + CALL FIELD_COMPONENT_MESH_COMPONENT_GET(FIELD,VARIABLE_TYPE,component_idx,MESH_COMPONENT2,ERR,ERROR,*999) + IS_ALL_NODAL_INTERPOLATION=IS_ALL_NODAL_INTERPOLATION.AND.INTERPOLATION_TYPE==FIELD_NODE_BASED_INTERPOLATION + IS_SAME_MESH_COMPONENTS=IS_SAME_MESH_COMPONENTS.AND.MESH_COMPONENT2==MESH_COMPONENT1 + ENDDO !component_idx + IF(IS_ALL_NODAL_INTERPOLATION) THEN + IF(IS_SAME_MESH_COMPONENTS) THEN + CALL FIELDML_INPUT_FIELD_NODAL_PARAMETERS_UPDATE(FIELDML_INFO,EVALUATOR_NAME,FIELD,VARIABLE_TYPE,SET_TYPE, & + & ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR( & + & "FieldML input parameters only implemented for fields where all components have the same mesh component.", & + & ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("FieldML input parameters only implemented for fields where all components are nodally interpolated.", & + & ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("Field does not have any components.",ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("Field is not associated.",ERR,ERROR,*999) + ENDIF + + CALL EXITS("FIELDML_INPUT_FIELD_PARAMETERS_UPDATE") + RETURN +999 CALL ERRORS("FIELDML_INPUT_FIELD_PARAMETERS_UPDATE",ERR,ERROR) + CALL EXITS("FIELDML_INPUT_FIELD_PARAMETERS_UPDATE") + RETURN 1 + END SUBROUTINE FIELDML_INPUT_FIELD_PARAMETERS_UPDATE + + ! + !================================================================================================================================ + ! + !>Update the given field's nodal parameters using the given parameter evaluator. SUBROUTINE FIELDML_INPUT_FIELD_NODAL_PARAMETERS_UPDATE( FIELDML_INFO, EVALUATOR_NAME, FIELD, VARIABLE_TYPE, SET_TYPE, & & ERR, ERROR, * ) diff --git a/src/mesh_routines.f90 b/src/mesh_routines.f90 index 25dac156..73a6b0dd 100644 --- a/src/mesh_routines.f90 +++ b/src/mesh_routines.f90 @@ -3235,16 +3235,28 @@ SUBROUTINE DECOMPOSITION_USER_NUMBER_FIND(USER_NUMBER,MESH,DECOMPOSITION,ERR,ERR TYPE(VARYING_STRING) :: LOCAL_ERROR CALL ENTERS("DECOMPOSITION_USER_NUMBER_FIND",ERR,ERROR,*999) + + PRINT *,MESH%DECOMPOSITIONS%NUMBER_OF_DECOMPOSITIONS + NULLIFY(DECOMPOSITION) IF(ASSOCIATED(MESH)) THEN IF(ASSOCIATED(MESH%DECOMPOSITIONS)) THEN decomposition_idx=1 DO WHILE(decomposition_idx<=MESH%DECOMPOSITIONS%NUMBER_OF_DECOMPOSITIONS.AND..NOT.ASSOCIATED(DECOMPOSITION)) - IF(MESH%DECOMPOSITIONS%DECOMPOSITIONS(decomposition_idx)%PTR%USER_NUMBER==USER_NUMBER) THEN - DECOMPOSITION=>MESH%DECOMPOSITIONS%DECOMPOSITIONS(decomposition_idx)%PTR + PRINT *,decomposition_idx + PRINT *,LBOUND(MESH%DECOMPOSITIONS%DECOMPOSITIONS) + PRINT *,SIZE(MESH%DECOMPOSITIONS%DECOMPOSITIONS) + PRINT *,UBOUND(MESH%DECOMPOSITIONS%DECOMPOSITIONS) + IF(ASSOCIATED(MESH%DECOMPOSITIONS%DECOMPOSITIONS)) THEN + IF(MESH%DECOMPOSITIONS%DECOMPOSITIONS(decomposition_idx)%PTR%USER_NUMBER==USER_NUMBER) THEN + PRINT *,"WTF" + DECOMPOSITION=>MESH%DECOMPOSITIONS%DECOMPOSITIONS(decomposition_idx)%PTR + ELSE + decomposition_idx=decomposition_idx+1 + ENDIF ELSE - decomposition_idx=decomposition_idx+1 + PRINT *,"WTF" ENDIF ENDDO ELSE diff --git a/src/opencmiss.f90 b/src/opencmiss.f90 index 639bea7e..50b9261f 100644 --- a/src/opencmiss.f90 +++ b/src/opencmiss.f90 @@ -2318,7 +2318,6 @@ MODULE OPENCMISS & EQUATIONS_SET_QUADRATIC_SOURCE_DIFFUSION_EQUATION_ONE_DIM_1 INTEGER(INTG), PARAMETER :: CMISSEquationsSetExponentialSourceDiffusionOneDim1 = & & EQUATIONS_SET_EXPONENTIAL_SOURCE_DIFFUSION_EQUATION_ONE_DIM_1 - INTEGER(INTG), PARAMETER :: CMISSEquationsSetMultiCompDiffusionTwoCompTwoDim = & & EQUATIONS_SET_MULTI_COMP_DIFFUSION_TWO_COMP_TWO_DIM !> + END INTERFACE !CMISSInterfaceMeshConnectivityElementXiSet !>Sets the number of elements coupled through a given interface mesh element INTERFACE CMISSInterfaceMeshConnectivityElementNumberSet @@ -5242,6 +5240,13 @@ MODULE OPENCMISS MODULE PROCEDURE CMISSSolverDAETimeStepSetObj END INTERFACE !CMISSSolverDAETimeStepSet + !>Sets/changes the parameters relating to an external differential-algebraic equation solver. + INTERFACE CMISSSolverExternalDAESolverParametersSet + MODULE PROCEDURE CMISSSolverExternalDAESolverParametersSetNumber0 + MODULE PROCEDURE CMISSSolverExternalDAESolverParametersSetNumber1 + MODULE PROCEDURE CMISSSolverExternalDAESolverParametersSetObj + END INTERFACE !CMISSSolverExternalDAESolverParametersSet + !>Returns the degree of the polynomial used to interpolate time for a dynamic solver. INTERFACE CMISSSolverDynamicDegreeGet MODULE PROCEDURE CMISSSolverDynamicDegreeGetNumber0 @@ -5629,7 +5634,7 @@ MODULE OPENCMISS PUBLIC CMISSSolverDAESolverTypeGet,CMISSSolverDAESolverTypeSet - PUBLIC CMISSSolverDAETimesSet,CMISSSolverDAETimeStepSet + PUBLIC CMISSSolverDAETimesSet,CMISSSolverDAETimeStepSet,CMISSSolverExternalDAESolverParametersSet PUBLIC CMISSSolverDynamicDegreeGet,CMISSSolverDynamicDegreeSet @@ -5739,13 +5744,13 @@ MODULE OPENCMISS MODULE PROCEDURE CMISSFieldMLInputFieldCreateStartNumberC END INTERFACE CMISSFieldMLInputFieldCreateStart - !> Updates the given field's nodal dofs using the given parameter evaluator. - INTERFACE CMISSFieldMLInputFieldNodalParametersUpdate - MODULE PROCEDURE CMISSFieldMLInputFieldNodalParametersUpdateObjVS - MODULE PROCEDURE CMISSFieldMLInputFieldNodalParametersUpdateNumberVS - MODULE PROCEDURE CMISSFieldMLInputFieldNodalParametersUpdateObjC - MODULE PROCEDURE CMISSFieldMLInputFieldNodalParametersUpdateNumberC - END INTERFACE CMISSFieldMLInputFieldNodalParametersUpdate + !> Updates the given field's dofs using the given parameter evaluator. + INTERFACE CMISSFieldMLInputFieldParametersUpdate + MODULE PROCEDURE CMISSFieldMLInputFieldParametersUpdateObjVS + MODULE PROCEDURE CMISSFieldMLInputFieldParametersUpdateNumberVS + MODULE PROCEDURE CMISSFieldMLInputFieldParametersUpdateObjC + MODULE PROCEDURE CMISSFieldMLInputFieldParametersUpdateNumberC + END INTERFACE CMISSFieldMLInputFieldParametersUpdate !> Creates a basis using the given FieldML evaluator. INTERFACE CMISSFieldMLInputBasisCreateStart @@ -5788,7 +5793,7 @@ MODULE OPENCMISS PUBLIC :: CMISSFieldMLInputCreateFromFile, CMISSFieldMLInputMeshCreateStart, & & CMISSFieldMLInputCoordinateSystemCreateStart, CMISSFieldMLInputCreateMeshComponent, & & CMISSFieldMLInputFieldCreateStart, CMISSFieldMLInputBasisCreateStart, CMISSFieldMLInputNodesCreateStart, & - & CMISSFieldMLInputFieldNodalParametersUpdate + & CMISSFieldMLInputFieldParametersUpdate PUBLIC :: CMISSFieldMLIOTypeFinalise, CMISSFieldMLIOTypeInitialise, CMISSFieldMLIOGetSession @@ -11364,7 +11369,7 @@ SUBROUTINE CMISSExtractErrorMessageC(ErrorMessage,Err) CALL EXTRACT_ERROR_MESSAGE(ErrorMessage,Err,ERROR,*999) RETURN -999 RETURN 1 +999 RETURN END SUBROUTINE CMISSExtractErrorMessageC @@ -11383,7 +11388,7 @@ SUBROUTINE CMISSExtractErrorMessageVS(ErrorMessage,Err) CALL EXTRACT_ERROR_MESSAGE(ErrorMessage,Err,ERROR,*999) RETURN -999 RETURN 1 +999 RETURN END SUBROUTINE CMISSExtractErrorMessageVS @@ -20328,7 +20333,6 @@ SUBROUTINE CMISSEquationsSetAnalyticDestroyObj(EquationsSet,Err) END SUBROUTINE CMISSEquationsSetAnalyticDestroyObj - ! !================================================================================================================================ ! @@ -20540,7 +20544,6 @@ SUBROUTINE CMISSEquationsSetAnalyticTimeSetObj(EquationsSet,Time,Err) END SUBROUTINE CMISSEquationsSetAnalyticTimeSetObj - ! !================================================================================================================================ ! @@ -30029,6 +30032,85 @@ SUBROUTINE CMISSFieldIONodesExportVSVSObj(Fields,FileName,Method,Err) END SUBROUTINE CMISSFieldIONodesExportVSVSObj + ! + !================================================================================================================================ + ! + + !> Import mesh from file \todo number method + SUBROUTINE CMISSFieldIOFieldsImport(FileName, Method, Region, Mesh, MeshUserNumber, Decomposition, & + & DecompositionUserNumber, DecompositionMethod, FieldValuesSetType, FieldScalingType, Err) + !Argument variables + CHARACTER(LEN=*), INTENT(IN):: FileName ! Import mesh from vtk file + SUBROUTINE CMISSReadVTK(filePath, numberOfMeshDimensions, nodesPerElement, points, numPoints, cells, numCells, Err) + !Argument variables + USE CMISS_FORTRAN_C + USE ISO_C_BINDING + CHARACTER(LEN=*), INTENT(IN) :: filePath + INTEGER(INTG), INTENT(IN) :: numberOfMeshDimensions + INTEGER(INTG), INTENT(IN) :: nodesPerElement + REAL(DP), POINTER, INTENT(OUT) :: points(:) + INTEGER(INTG), INTENT(OUT) :: numPoints + INTEGER(INTG), POINTER, INTENT(OUT) :: cells(:) + INTEGER(INTG), INTENT(OUT) :: numCells + INTEGER(INTG), INTENT(OUT) :: Err !Sets/changes the parameters realting to an external differential-algebraic equation solver identified by an user number. + SUBROUTINE CMISSSolverExternalDAESolverParametersSetNumber0(ProblemUserNumber,ControlLoopIdentifier,SolverIndex, & + & ThreadsPerBlock,NumberOfPartitons,NumberOfStreams,Err) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: ProblemUserNumber !Sets/changes the parameters realting to an external differential-algebraic equation solver identified by an user number. + SUBROUTINE CMISSSolverExternalDAESolverParametersSetNumber1(ProblemUserNumber,ControlLoopIdentifiers,SolverIndex, & + & ThreadsPerBlock,NumberOfPartitons,NumberOfStreams,Err) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: ProblemUserNumber !Sets/changes the parameters realting to an external differential-algebraic equation solver identified by an object. + SUBROUTINE CMISSSolverExternalDAESolverParametersSetObj(Solver,ThreadsPerBlock,NumberOfPartitons,NumberOfStreams,Err) + + !Argument variables + TYPE(CMISSSolverType), INTENT(IN) :: Solver !Returns the degree of the polynomial used to interpolate time for a dynamic solver identified by an user number. SUBROUTINE CMISSSolverDynamicDegreeGetNumber0(ProblemUserNumber,ControlLoopIdentifier,SolverIndex,Degree,Err) @@ -47414,8 +47612,8 @@ END SUBROUTINE CMISSFieldMLInputFieldCreateStartNumberC !================================================================================================================================ ! - !> Update the nodal parameters of the given field, using the given FieldML evaluator. - SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateObjVS( fieldml, field, evaluatorName, variableType, & + !> Update the DOF parameters of the given field, using the given FieldML evaluator. + SUBROUTINE CMISSFieldMLInputFieldParametersUpdateObjVS( fieldml, field, evaluatorName, variableType, & & setType, err ) !Arguments TYPE(CMISSFieldMLIOType), INTENT(INOUT) :: FIELDML !< The FieldML context containing the evaluator to use. @@ -47425,26 +47623,26 @@ SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateObjVS( fieldml, field, eva INTEGER(INTG), INTENT(IN) :: setType ! Update the nodal parameters of field with the given user number, using the given FieldML evaluator. - SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateNumberVS( fieldml, regionNumber, fieldNumber, & + !> Update the DOF parameters of field with the given user number, using the given FieldML evaluator. + SUBROUTINE CMISSFieldMLInputFieldParametersUpdateNumberVS( fieldml, regionNumber, fieldNumber, & & evaluatorName, variableType, setType, err ) !Arguments TYPE(CMISSFieldMLIOType), INTENT(INOUT) :: FIELDML !< The FieldML context containing the evaluator to use. @@ -47459,29 +47657,29 @@ SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateNumberVS( fieldml, regionN TYPE(REGION_TYPE), POINTER :: region TYPE(FIELD_TYPE), POINTER :: field - CALL ENTERS("CMISSFieldMLInputFieldNodalParametersUpdateNumberVS",Err,ERROR,*999) + CALL ENTERS("CMISSFieldMLInputFieldParametersUpdateNumberVS",Err,ERROR,*999) CALL REGION_USER_NUMBER_TO_REGION( regionNumber, region, err, error, *999 ) CALL FIELD_USER_NUMBER_TO_FIELD( fieldNumber, region, field, err, error, *999 ) - CALL FIELDML_INPUT_FIELD_NODAL_PARAMETERS_UPDATE( fieldml%fieldmlInfo, evaluatorName, field, variableType, setType, & + CALL FIELDML_INPUT_FIELD_PARAMETERS_UPDATE( fieldml%fieldmlInfo, evaluatorName, field, variableType, setType, & & err, error, *999 ) - CALL EXITS("CMISSFieldMLInputFieldNodalParametersUpdateNumberVS") + CALL EXITS("CMISSFieldMLInputFieldParametersUpdateNumberVS") RETURN -999 CALL ERRORS("CMISSFieldMLInputFieldNodalParametersUpdateNumberVS",Err,ERROR) - CALL EXITS("CMISSFieldMLInputFieldNodalParametersUpdateNumberVS") +999 CALL ERRORS("CMISSFieldMLInputFieldParametersUpdateNumberVS",Err,ERROR) + CALL EXITS("CMISSFieldMLInputFieldParametersUpdateNumberVS") CALL CMISS_HANDLE_ERROR(Err,ERROR) RETURN - END SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateNumberVS + END SUBROUTINE CMISSFieldMLInputFieldParametersUpdateNumberVS ! !================================================================================================================================ ! - !> Update the nodal parameters of the given field, using the given FieldML evaluator. - SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateObjC( fieldml, field, evaluatorName, & + !> Update the DOF parameters of the given field, using the given FieldML evaluator. + SUBROUTINE CMISSFieldMLInputFieldParametersUpdateObjC( fieldml, field, evaluatorName, & & variableType, setType, err ) !Arguments TYPE(CMISSFieldMLIOType), INTENT(INOUT) :: FIELDML !< The FieldML context containing the evaluator to use. @@ -47491,26 +47689,26 @@ SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateObjC( fieldml, field, eval INTEGER(INTG), INTENT(IN) :: setType ! Update the nodal parameters of field with the given user number, using the given FieldML evaluator. - SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateNumberC( fieldml, regionNumber, fieldNumber, & + !> Update the DOF parameters of field with the given user number, using the given FieldML evaluator. + SUBROUTINE CMISSFieldMLInputFieldParametersUpdateNumberC( fieldml, regionNumber, fieldNumber, & & evaluatorName, variableType, setType, err ) !Arguments TYPE(CMISSFieldMLIOType), INTENT(INOUT) :: FIELDML !< The FieldML context containing the evaluator to use. @@ -47525,22 +47723,22 @@ SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateNumberC( fieldml, regionNu TYPE(REGION_TYPE), POINTER :: region TYPE(FIELD_TYPE), POINTER :: field - CALL ENTERS("CMISSFieldMLInputFieldNodalParametersUpdateNumberC",Err,ERROR,*999) + CALL ENTERS("CMISSFieldMLInputFieldParametersUpdateNumberC",Err,ERROR,*999) CALL REGION_USER_NUMBER_TO_REGION( regionNumber, region, err, error, *999 ) CALL FIELD_USER_NUMBER_TO_FIELD( fieldNumber, region, field, err, error, *999 ) - CALL FIELDML_INPUT_FIELD_NODAL_PARAMETERS_UPDATE( fieldml%fieldmlInfo, var_str(evaluatorName), field, variableType, & + CALL FIELDML_INPUT_FIELD_PARAMETERS_UPDATE( fieldml%fieldmlInfo, var_str(evaluatorName), field, variableType, & & setType, err, error, *999 ) - CALL EXITS("CMISSFieldMLInputFieldNodalParametersUpdateNumberC") + CALL EXITS("CMISSFieldMLInputFieldParametersUpdateNumberC") RETURN -999 CALL ERRORS("CMISSFieldMLInputFieldNodalParametersUpdateNumberC",Err,ERROR) - CALL EXITS("CMISSFieldMLInputFieldNodalParametersUpdateNumberC") +999 CALL ERRORS("CMISSFieldMLInputFieldParametersUpdateNumberC",Err,ERROR) + CALL EXITS("CMISSFieldMLInputFieldParametersUpdateNumberC") CALL CMISS_HANDLE_ERROR(Err,ERROR) RETURN - END SUBROUTINE CMISSFieldMLInputFieldNodalParametersUpdateNumberC + END SUBROUTINE CMISSFieldMLInputFieldParametersUpdateNumberC ! !================================================================================================================================ @@ -47553,8 +47751,6 @@ SUBROUTINE CMISSFieldMLOutputWriteVS( fieldml, filename, err ) TYPE(VARYING_STRING), INTENT(IN) :: filename !< The name of the file to write the FieldML document to. INTEGER(INTG), INTENT(OUT) :: err !< The error code. - TYPE(VARYING_STRING) :: LOCAL_ERROR - CALL ENTERS("CMISSFieldMLOutputWriteVS",Err,ERROR,*999) IF( .NOT. fieldml%fieldmlInfo%IS_OUT ) THEN @@ -48245,7 +48441,7 @@ SUBROUTINE CMISSFieldMLIOGetSession( fieldml, sessionHandle, Err ) END SUBROUTINE CMISSFieldMLIOGetSession -#endif !USEFIELDML +#endif ! !================================================================================================================================ diff --git a/src/problem_routines.f90 b/src/problem_routines.f90 index e1ea439f..75442c2c 100644 --- a/src/problem_routines.f90 +++ b/src/problem_routines.f90 @@ -319,7 +319,6 @@ SUBROUTINE PROBLEM_CELLML_EQUATIONS_SOLVE(CELLML_EQUATIONS,ERR,ERROR,*) TYPE(VARYING_STRING), INTENT(OUT) :: ERROR !DAE_SOLVER !Defaults + DAE_SOLVER%EXTERNAL_SOLVER%THREADS_PER_BLOCK=512 + DAE_SOLVER%EXTERNAL_SOLVER%NUMBER_OF_PARTITIONS=1 + DAE_SOLVER%EXTERNAL_SOLVER%NUMBER_OF_STREAMS=1 ENDIF ELSE CALL FLAG_ERROR("Differential-algebraic equation solver is not associated.",ERR,ERROR,*998) @@ -3275,6 +3286,9 @@ SUBROUTINE SOLVER_DAE_EXTERNAL_SOLVE(EXTERNAL_SOLVER,ERR,ERROR,*) TYPE(SOLVER_TYPE), POINTER :: SOLVER TYPE(VARYING_STRING) :: LOCAL_ERROR + INTEGER(INTG) :: N + + CALL ENTERS("SOLVER_DAE_EXTERNAL_SOLVE",ERR,ERROR,*999) IF(ASSOCIATED(EXTERNAL_SOLVER)) THEN @@ -3336,11 +3350,20 @@ SUBROUTINE SOLVER_DAE_EXTERNAL_SOLVE(EXTERNAL_SOLVER,ERR,ERROR,*) ENDIF ELSE NULLIFY(INTERMEDIATE_DATA) + NULLIFY(INTERMEDIATE_FIELD) ENDIF + !PRINT *,STATE_DATA(4470) + !DO N=1,MODELS_VARIABLE%TOTAL_NUMBER_OF_DOFS + ! IF (IsNaN(STATE_DATA(N))) THEN + ! PRINT *,N,STATE_DATA(N) + ! ENDIF + !ENDDO + !Call the external solver to integrate these CellML equations CALL SOLVER_DAE_EXTERNAL_INTEGRATE(MODELS_VARIABLE%TOTAL_NUMBER_OF_DOFS,DAE_SOLVER%START_TIME, & - & DAE_SOLVER%END_TIME,DAE_SOLVER%INITIAL_STEP,CELLML_ENVIRONMENT%MODELS_FIELD% & + & DAE_SOLVER%END_TIME,DAE_SOLVER%INITIAL_STEP,EXTERNAL_SOLVER%THREADS_PER_BLOCK,EXTERNAL_SOLVER% & + & NUMBER_OF_PARTITIONS,EXTERNAL_SOLVER%NUMBER_OF_STREAMS,CELLML_ENVIRONMENT%MODELS_FIELD% & & ONLY_ONE_MODEL_INDEX,MODELS_DATA,CELLML_ENVIRONMENT%MAXIMUM_NUMBER_OF_STATE,STATE_DATA, & & CELLML_ENVIRONMENT%MAXIMUM_NUMBER_OF_PARAMETERS,PARAMETERS_DATA,CELLML_ENVIRONMENT% & & MAXIMUM_NUMBER_OF_INTERMEDIATE,INTERMEDIATE_DATA,ERR) @@ -3614,6 +3637,13 @@ SUBROUTINE SOLVER_DAE_SOLVE(DAE_SOLVER,ERR,ERROR,*) TYPE(SOLVER_TYPE), POINTER :: SOLVER TYPE(VARYING_STRING) :: LOCAL_ERROR + + ! TEMP ADDED BY JIN + REAL TEMP, ODE_TIME + REAL TIMEARRAY(2) + COMMON ODE_TIME + TEMP = DTIME(TIMEARRAY) + CALL ENTERS("SOLVER_DAE_SOLVE",ERR,ERROR,*999) IF(ASSOCIATED(DAE_SOLVER)) THEN @@ -3680,6 +3710,10 @@ SUBROUTINE SOLVER_DAE_SOLVE(DAE_SOLVER,ERR,ERROR,*) ELSE CALL FLAG_ERROR("Differential-algebraic equation solver is not associated.",ERR,ERROR,*999) ENDIF + + + ! TEMP CODE ADDED BY JIN + ODE_TIME = ODE_TIME + DTIME(TIMEARRAY) CALL EXITS("SOLVER_DAE_SOLVE") RETURN @@ -3918,6 +3952,56 @@ SUBROUTINE SOLVER_DAE_TIME_STEP_SET(SOLVER,TIME_STEP,ERR,ERROR,*) END SUBROUTINE SOLVER_DAE_TIME_STEP_SET + ! + !================================================================================================================================ + ! + + !>Set/change the parameters realting to the external differential-algebraic equation solver + SUBROUTINE SOLVER_EXTERNAL_DAE_PARAMETERS_SET(SOLVER,THREADS_PER_BLOCK,NUMBER_OF_PARTITIONS,NUMBER_OF_STREAMS,ERR,ERROR,*) + + !Argument variables + TYPE(SOLVER_TYPE), POINTER :: SOLVER !SOLVER%DAE_SOLVER + IF(ASSOCIATED(DAE_SOLVER)) THEN + EXTERNAL_DAE_SOLVER=>DAE_SOLVER%EXTERNAL_SOLVER + IF(ASSOCIATED(EXTERNAL_DAE_SOLVER)) THEN + EXTERNAL_DAE_SOLVER%THREADS_PER_BLOCK=THREADS_PER_BLOCK !\TODO validate data entry + EXTERNAL_DAE_SOLVER%NUMBER_OF_PARTITIONS=NUMBER_OF_PARTITIONS + EXTERNAL_DAE_SOLVER%NUMBER_OF_STREAMS=NUMBER_OF_STREAMS + ENDIF + ELSE + CALL FLAG_ERROR("Differential-algebraic equation solver is not associated.",ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("The solver is not a differential-algebraic equation solver.",ERR,ERROR,*999) + ENDIF + ELSE + CALL FLAG_ERROR("Solver is not associated.",ERR,ERROR,*999) + ENDIF + + CALL EXITS("SOLVER_EXTERNAL_DAE_PARAMETERS_SET") + RETURN +999 CALL ERRORS("SOLVER_EXTERNAL_DAE_PARAMETERS_SET",ERR,ERROR) + CALL EXITS("SOLVER_EXTERNAL_DAE_PARAMETERS_SET") + RETURN 1 + + END SUBROUTINE SOLVER_EXTERNAL_DAE_PARAMETERS_SET + + ! !================================================================================================================================ ! diff --git a/src/types.f90 b/src/types.f90 index d706775b..f531b100 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -2367,6 +2367,9 @@ MODULE TYPES !>Contains information for an external differential-algebraic equation solver TYPE EXTERNAL_DAE_SOLVER_TYPE TYPE(DAE_SOLVER_TYPE), POINTER :: DAE_SOLVER !Contains information for an differential-algebraic equation solver diff --git a/src/vtk_import_routines.c b/src/vtk_import_routines.c new file mode 100644 index 00000000..29e1ad10 --- /dev/null +++ b/src/vtk_import_routines.c @@ -0,0 +1,113 @@ +#include +#include +#include + +#ifdef _MSC_VER + #define FSCANF fscanf_s +#else + #define FSCANF fscanf +#endif + +extern int whiteSpaceCounter(FILE*); + +const char* cellString = "CELLS"; + +void readVTK(const char* filePath, int numberOfMeshDimensions, int nodesPerElement, double** points, int* numPoints, int** cells, int* numCells){ + int numLine = 0; + int index = 0; + int i; + int count = 0; + + char buffer[100]; + + FILE* file = NULL; + +#ifdef _MSC_VER + fopen_s(&file,filePath, "r"); +#else + file = fopen(filePath, "r"); +#endif + + // Trash first 4 lines on file + fgets(buffer, 100, file); + fgets(buffer, 100, file); + fgets(buffer, 100, file); + fgets(buffer, 100, file); + + FSCANF(file, "%*s %d %*s\n", numPoints); + + printf("Number of Points = %d\n", (*numPoints)); + + (*points) = (double *) malloc(sizeof(double)*(*numPoints)*numberOfMeshDimensions); + + while (count<(*numPoints)*numberOfMeshDimensions) { + numLine = whiteSpaceCounter(file); + for (i=0; i