OSCR

Diversity in emergent cell locomotion from the coupling cytosolic and cortical Marangoni flows with reaction-diffusion dynamics.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

C · 1,327 lines · 56 KB · MIT

  1. /* SPDX-License-Identifier: MIT */
  2. /*
  3. * Project: CellRDA-LS (Mechanochemical cell motility simulations)
  4. * File: main.c
  5. * Brief: Program entry point; parses JSON parameters and runs batch simulations.
  6. *
  7. * Author: Blaž Ivšić
  8. * Contact: [email hidden]
  9. * Copyright: (c) 2023–2025 Blaž Ivšić
  10. *
  11. * License: MIT — see LICENSE in the repository root.
  12. *
  13. * How to cite:
  14. * Article: Diversity in emergent cell locomotion from the coupling cytosolic and cortical Marangoni flows with reaction–diffusion dynamics, <Journal or preprint>, <Year>. DOI: <paper DOI>
  15. * Software: Cite this release via CITATION.cff (and Zenodo DOI if minted).
  16. *
  17. * Reproducibility:
  18. * - Language/Std: C11
  19. * - Compiler: gcc >= 12 and/or clang >= 16 (tested)
  20. * - Dependencies: NetCDF >= 4.9 (with HDF5), cJSON (vendored), -lm
  21. * - Optional: OpenMP, AVX2 (guard with feature flags)
  22. * - Inputs: SimulationParameters_X.json, InitialConditions_X.nc
  23. * - Outputs: <files>.nc/.bin/.txt (see README)
  24. * - Determinism: set Seed=<int> to reproduce stochastic runs
  25. *
  26. * Version: v1.0.0 | Commit: <git short SHA> | Date: <YYYY-MM-DD>
  27. * Build (ex): gcc -O3 -std=c11 -fopenmp -Ithird_party/cJSON \
  28. * src/*.c third_party/cJSON/cJSON.c \
  29. * -lnetcdf -lhdf5 -lm -o program
  30. */
  33. #include <stdio.h>
  34. #include <stdlib.h>
  35. #include <netcdf.h>
  36. #include <string.h>
  37. #include <math.h>
  38. #include <time.h>
  40. #include "cJSON.h"
  41. #include "simulation.h"
  43. #define chcJSONitem(item, name) \
  44. if ((item) == NULL) { \
  45. printf("Error: Missing or invalid value for key: %s\n", name); \
  46. return 11; \
  47. }
  48. #define ALLOC_ARRAY(name, type, rows, cols) \
  49. type *name = (type *)calloc((rows) * (cols), sizeof(type)); \
  50. if (name == NULL) \
  51. { \
  52. printf("Memory allocation failed for " #name "\n"); \
  53. \
  54. exit(EXIT_FAILURE); \
  55. }
  57. char *readJSON(const char *filename)
  58. {
  59. //Function for loading JSON file
  60. FILE *file = fopen(filename, "r");
  61. if(!file)
  62. {
  63. printf("Error: Failed to open the file\n");
  64. return NULL;
  65. }
  66. // Determining length of the file
  67. fseek(file, 0, SEEK_END);
  68. long length = ftell(file);
  69. fseek(file, 0, SEEK_SET);
  70. // Memory allocation
  71. char *content = (char *)malloc(length + 1); // While this makes the cast explicit, it's not required in C, and omitting the cast is typically the preferred approach in modern C programming.
  72. // ChatGPT: This was common in older versions of C or when using C++ compilers, which are stricter about type conversions. However, in modern C compilers, the cast is generally not necessary because malloc returns a void *, which can be assigned to any pointer type without a cast.
  73. fread(content, 1, length, file);
  74. fclose(file);
  75. content[length] = '\0';
  76. return content;
  77. }
  78. cJSON* loadJSON(const char* filename)
  79. {
  80. // Function for loading the JSON data
  81. char *JSONcontent = readJSON(filename);
  82. if (!JSONcontent)
  83. {
  84. printf("Error: Failed to read the file\n");
  85. return NULL;
  86. }
  87. // Parsing file content
  88. cJSON *json = cJSON_Parse(JSONcontent);
  89. free(JSONcontent);
  90. if (json == NULL) {
  91. printf("Error: Failed to parse JSON\n");
  92. return NULL;
  93. }
  94. return json;
  95. }
  96. int validateArray(cJSON *array, int expectedSize, const char *keyName) {
  97. if (array == NULL || array->type != cJSON_Array || cJSON_GetArraySize(array) != expectedSize) {
  98. printf("Error: Missing or invalid array for key: %s (expected size: %d)\n", keyName, expectedSize);
  99. return 11; // Error code for missing or invalid array
  100. }
  101. for (int i = 0; i < expectedSize; i++) {
  102. cJSON *item = cJSON_GetArrayItem(array, i);
  103. if (item == NULL || item->type != cJSON_Number) {
  104. printf("Error: Invalid value at index %d in array for key: %s\n", i, keyName);
  105. return 12; // Error code for invalid array value
  106. }
  107. }
  108. return 0; // Validation successful
  109. }
  110. int validateParameters(cJSON *json) {
  111. // Function for checking JSON keys
  112. // File parameters
  113. cJSON *FileParameters = cJSON_GetObjectItem(json, "FileParameters");
  114. chcJSONitem(FileParameters, "FileParameters");
  115. chcJSONitem(cJSON_GetObjectItem(FileParameters, "InitialFile"), "InitialFile");
  116. chcJSONitem(cJSON_GetObjectItem(FileParameters, "SaveFile"), "SaveFile");
  117. chcJSONitem(cJSON_GetObjectItem(FileParameters, "TracerFolder"), "TracerFolder");
  118. chcJSONitem(cJSON_GetObjectItem(FileParameters, "TracerPrefix"), "TracerPrefix");
  119. // Simulation parameters
  120. cJSON *SimulationParameters = cJSON_GetObjectItem(json, "SimulationParameters");
  121. chcJSONitem(SimulationParameters, "SimulationParameters");
  122. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "Cluster"), "Cluster");
  123. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "Seed"), "Seed");
  124. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "NumberOfTracersBulk"), "NumberOfTracersBulk");
  125. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "NumberOfTracersOut"), "NumberOfTracersOut");
  126. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "NumberOfSaveCycles"), "NumberOfSaveCycles");
  127. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "NumberOfDoubleStep"), "NumberOfDoubleStep");
  128. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "NumberOfPointsX"), "NumberOfPointsX");
  129. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "NumberOfPointsY"), "NumberOfPointsY");
  130. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "GridSpacingX"), "GridSpacingX");
  131. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "GridSpacingY"), "GridSpacingY");
  132. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "TimeStepSpacing"), "TimeStepSpacing");
  133. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "CutOff"), "CutOff");
  134. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "PotentialCutOffBulk"), "PotentialCutOffBulk");
  135. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "PotentialCutOffSurface"), "PotentialCutOffSurface");
  136. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "MaxIterations"), "MaxIterations");
  137. chcJSONitem(cJSON_GetObjectItem(SimulationParameters, "PressTol"), "PressTol");
  138. // Model Parameters
  139. cJSON *ModelParameters = cJSON_GetObjectItem(json, "ModelParameters");
  140. chcJSONitem(ModelParameters, "ModelParameters");
  141. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "Epsilon"), "Epsilon");
  142. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "Viscosity"), "Viscosity");
  143. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "Gamma"), "Gamma");
  144. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "SurfaceTension"), "SurfaceTension");
  145. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "SurfaceTensionCoefficient"), "SurfaceTensionCoefficient");
  146. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "SubstrateForceCoefficient"), "SubstrateForceCoefficient");
  147. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "Dissipation"), "Dissipation");
  148. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "AreaConstraint"), "AreaConstraint");
  149. chcJSONitem(cJSON_GetObjectItem(ModelParameters, "Rac1Max"), "Rac1Max");
  151. cJSON *DiffusionArray = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Diffusion"), "value");
  152. cJSON *ReactionsArray = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Reactions"), "value");
  153. int chcDiff = validateArray(DiffusionArray, 4, "Diffusion");
  154. if (chcDiff != 0) {return chcDiff;}
  155. int chcReac = validateArray(ReactionsArray, 5, "Reactions");
  156. if (chcReac != 0) {return chcReac;}
  158. return 0; // Validation successful
  159. }
  160. int loadParameters(cJSON *json,
  161. char *InitialConditionsFile, char *SaveFileName, char *TracerFolder, char * TracerPrefix,
  162. int *Seed, int *Cluster, int *nTracersB, int *nTracersO, int *NumberOfSaveCycles, int *NumberOfDoubleStep,
  163. size_t *I, size_t *J,
  164. double *dx, double *dy, double *dt, double *CO, double *vBB, double *vBS, int *MaxIter, double *PTol,
  165. double *eps, double *nu, double *cG, double *Si0, double *kSi, double *kC, double *Beta, double *Alpha,
  166. double *DRD, double *DRT, double *DG, double *DGc,
  167. double *k1, double *k2, double *k3, double *k11, double *k12, double *R_Max)
  168. {
  169. // Function for loading parameters from .json content
  170. // The Parameters are split into groups depending on the type of data
  171. // If new parameters are added the type of parameter to be loaded needs to be specified accordingly even dough the type is specified in .json file.
  172. // Checking the compatibility of JOSN file
  173. int chcVAL = validateParameters(json);
  174. if (chcVAL != 0) {return chcVAL;}
  175. // Loading JSON fields
  176. cJSON *FileParameters = cJSON_GetObjectItem(json, "FileParameters");
  177. cJSON *SimulationParameters = cJSON_GetObjectItem(json, "SimulationParameters");
  178. cJSON *ModelParameters = cJSON_GetObjectItem(json, "ModelParameters");
  180. // File parameters
  181. strcpy(InitialConditionsFile, cJSON_GetObjectItem(cJSON_GetObjectItem(FileParameters, "InitialFile"), "value")->valuestring);
  182. strcpy(SaveFileName, cJSON_GetObjectItem(cJSON_GetObjectItem(FileParameters, "SaveFile"), "value")->valuestring);
  183. strcpy(TracerFolder, cJSON_GetObjectItem(cJSON_GetObjectItem(FileParameters, "TracerFolder"), "value")->valuestring);
  184. strcpy(TracerPrefix, cJSON_GetObjectItem(cJSON_GetObjectItem(FileParameters, "TracerPrefix"), "value")->valuestring);
  185. // Simulation parameters
  186. *Cluster = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "Cluster"), "value")->valueint;
  187. *Seed = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "Seed"), "value")->valueint;
  188. *nTracersB = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "NumberOfTracersBulk"), "value")->valueint;
  189. *nTracersO = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "NumberOfTracersOut"), "value")->valueint;
  190. *NumberOfSaveCycles = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "NumberOfSaveCycles"), "value")->valueint;
  191. *NumberOfDoubleStep = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "NumberOfDoubleStep"), "value")->valueint;
  192. *I = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "NumberOfPointsX"), "value")->valueint;
  193. *J = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "NumberOfPointsY"), "value")->valueint;
  194. *dx = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "GridSpacingX"), "value")->valuedouble;
  195. *dy = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "GridSpacingY"), "value")->valuedouble;
  196. *dt = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "TimeStepSpacing"), "value")->valuedouble;
  197. *CO = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "CutOff"), "value")->valuedouble;
  198. *vBB = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "PotentialCutOffBulk"), "value")->valuedouble;
  199. *vBS = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "PotentialCutOffSurface"), "value")->valuedouble;
  200. *MaxIter = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "MaxIterations"), "value")->valueint;
  201. *PTol = cJSON_GetObjectItem(cJSON_GetObjectItem(SimulationParameters, "PressTol"), "value")->valuedouble;
  202. // Model Parameters
  203. *eps = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Epsilon"), "value")->valuedouble;
  204. *nu = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Viscosity"), "value")->valuedouble;
  205. *cG = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Gamma"), "value")->valuedouble;
  206. *Si0 = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "SurfaceTension"), "value")->valuedouble;
  207. *kSi = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "SurfaceTensionCoefficient"), "value")->valuedouble;
  208. *kC = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "SubstrateForceCoefficient"), "value")->valuedouble;
  209. *Beta = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Dissipation"), "value")->valuedouble;
  210. *Alpha = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "AreaConstraint"), "value")->valuedouble;
  211. cJSON *diffusion = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Diffusion"), "value");
  212. double Diffusion[4];
  213. for (int i = 0; i < 4; i++) {
  214. Diffusion[i] = cJSON_GetArrayItem(diffusion, i)->valuedouble;
  215. }
  216. *DRD = Diffusion[0];
  217. *DRT = Diffusion[1];
  218. *DG = Diffusion[2];
  219. *DGc = Diffusion[3];
  220. cJSON *reactions = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Reactions"), "value");
  221. double Reactions[5];
  222. for (int i = 0; i < 5; i++) {
  223. Reactions[i] = cJSON_GetArrayItem(reactions, i)->valuedouble;
  224. }
  225. *k1 = Reactions[0];
  226. *k2 = Reactions[1];
  227. *k3 = Reactions[2];
  228. *k11 = Reactions[3];
  229. *k12 = Reactions[4];
  230. *R_Max = cJSON_GetObjectItem(cJSON_GetObjectItem(ModelParameters, "Rac1Max"), "value")->valuedouble;
  231. return 0;
  232. }
  233. int checkSolutionFile(const char* SaveFileName)
  234. {
  235. // Function for checking if the solution exists and where has simulation terminated.
  236. int ncID, timeVarID;
  237. int RetVAL;
  238. size_t Len;
  240. // Try to open the solutions file
  241. RetVAL = nc_open(SaveFileName, NC_NOWRITE, &ncID);
  242. if (RetVAL == 2)
  243. {
  244. return -1; //The code that signals there is no Solution file
  245. }
  246. else if (RetVAL != NC_NOERR)
  247. {
  248. printf("Error opening file: %s\n", nc_strerror(RetVAL));
  249. return -30;
  250. }
  251. // Get the ID of the "Time" variable
  252. if (nc_inq_varid(ncID, "Time", &timeVarID) != NC_NOERR)
  253. {
  254. nc_close(ncID);
  255. return -31; // Time variable not found
  256. }
  257. // Get the length of the "Time" variable
  258. RetVAL = nc_inq_dimlen(ncID, timeVarID, &Len);
  259. if (RetVAL != NC_NOERR)
  260. {
  261. printf("Error: Time dimension not found in the file. %s\n", nc_strerror(RetVAL));
  262. nc_close(ncID);
  263. return -32; // Can't fetch the length of "Time" variable
  264. }
  265. // Allocate memory to hold the "Time" values
  266. double *TimeValues = (double *)malloc(Len * sizeof(double));
  267. if (TimeValues == NULL)
  268. {
  269. nc_close(ncID);
  270. return -98;
  271. }
  272. // Read the "Time" values
  273. if (nc_get_var_double(ncID, timeVarID, TimeValues) != NC_NOERR)
  274. {
  275. printf("Error reading Time vector: %s\n", nc_strerror(RetVAL));
  276. free(TimeValues);
  277. nc_close(ncID);
  278. return -33;
  279. }
  280. // Find the last index where the value is below the threshold
  281. int LastIndex = -1;
  282. for (size_t i = 0; i < Len; i++)
  283. {
  284. if (TimeValues[i] < 1e10)
  285. {
  286. LastIndex = i;
  287. }
  288. else
  289. {
  290. break;
  291. }
  292. }
  293. // Clean up
  294. free(TimeValues);
  295. nc_close(ncID);
  296. return LastIndex;
  297. }
  298. int readICNetCDF(const char *InitialConditionsFile)
  299. {
  300. // Function for opening NetCDF File
  301. int ncID, RetVAL;
  302. // Open the NetCDF file and assign the result to RetVAL
  303. RetVAL = nc_open(InitialConditionsFile, NC_NOWRITE, &ncID);
  304. if (RetVAL != NC_NOERR) // Check if an error occurred
  305. {
  306. printf("Error opening file: %s\n", nc_strerror(RetVAL));
  307. return -30;
  308. }
  309. // Return the NetCDF ID if successful
  310. return ncID;
  311. }
  312. int loadICNetCDF(int ncID, int *varID_Vx, int *varID_Vy, int *varID_Phi, int *varID_C_RD, int *varID_C_RT, int *varID_C_G, int *varID_C_Gc, int *varID_X, int *varID_Y)
  313. {
  314. // Function for loading NetCDF files for initial condition.
  315. // This function will return variable IDs and check if the needed initial condition matrices are included in the file.
  316. int RetVAL;
  317. RetVAL = nc_inq_varid(ncID, "Vx", varID_Vx);
  318. if (RetVAL != NC_NOERR)
  319. {
  320. printf("Error: Vx matrix not found in the file. %s\n", nc_strerror(RetVAL));
  321. return -31;
  322. }
  323. RetVAL = nc_inq_varid(ncID, "Vy", varID_Vy);
  324. if (RetVAL != NC_NOERR)
  325. {
  326. printf("Error: Vy matrix not found in the file. %s\n", nc_strerror(RetVAL));
  327. return -31;
  328. }
  329. RetVAL = nc_inq_varid(ncID, "Phi", varID_Phi);
  330. if (RetVAL != NC_NOERR)
  331. {
  332. printf("Error: Phi matrix not found in the file. %s\n", nc_strerror(RetVAL));
  333. return -31;
  334. }
  335. RetVAL = nc_inq_varid(ncID, "RD", varID_C_RD);
  336. if (RetVAL != NC_NOERR)
  337. {
  338. printf("Error: C_RD matrix not found in the file. %s\n", nc_strerror(RetVAL));
  339. return -31;
  340. }
  341. RetVAL = nc_inq_varid(ncID, "RT", varID_C_RT);
  342. if (RetVAL != NC_NOERR)
  343. {
  344. printf("Error: C_RT matrix not found in the file. %s\n", nc_strerror(RetVAL));
  345. return -31;
  346. }
  347. RetVAL = nc_inq_varid(ncID, "G", varID_C_G);
  348. if (RetVAL != NC_NOERR)
  349. {
  350. printf("Error: C_G matrix not found in the file. %s\n", nc_strerror(RetVAL));
  351. return -31;
  352. }
  353. RetVAL = nc_inq_varid(ncID, "Gc", varID_C_Gc);
  354. if (RetVAL != NC_NOERR)
  355. {
  356. printf("Error: C_Gc matrix not found in the file. %s\n", nc_strerror(RetVAL));
  357. return -31;
  358. }
  359. RetVAL = nc_inq_varid(ncID, "X", varID_X);
  360. if (RetVAL != NC_NOERR)
  361. {
  362. printf("Error: X vector not found in the file. %s\n", nc_strerror(RetVAL));
  363. return -31;
  364. }
  365. RetVAL = nc_inq_varid(ncID, "Y", varID_Y);
  366. if (RetVAL != NC_NOERR)
  367. {
  368. printf("Error: Y vector not found in the file. %s\n", nc_strerror(RetVAL));
  369. return -31;
  370. }
  371. return 0;
  372. }
  373. int checkICDimensions(int ncID, size_t I, size_t J)
  374. {
  375. // This function checks if the dimensions of initial conditions and ones specified in the parameters .json file are compatible
  376. int RetVAL;
  377. size_t FileI, FileJ;
  378. int xDim, yDim;
  380. RetVAL = nc_inq_dimid(ncID, "xDim", &xDim);
  381. if (RetVAL != NC_NOERR)
  382. {
  383. printf("Error: xDim dimenison not found in the file. %s\n", nc_strerror(RetVAL));
  384. return -32;
  385. }
  386. RetVAL = nc_inq_dimid(ncID, "yDim", &yDim);
  387. if (RetVAL != NC_NOERR)
  388. {
  389. printf("Error: yDim dimension not found in the file. %s\n", nc_strerror(RetVAL));
  390. return -32;
  391. }
  392. RetVAL = nc_inq_dimlen(ncID, xDim, &FileI);
  393. if (RetVAL != NC_NOERR)
  394. {
  395. printf("Error: Failed to obtain xDim %s\n", nc_strerror(RetVAL));
  396. return -32;
  397. }
  398. RetVAL = nc_inq_dimlen(ncID, yDim, &FileJ);
  399. if (RetVAL != NC_NOERR)
  400. {
  401. printf("Error: Failed to obtain yDim %s\n", nc_strerror(RetVAL));
  402. return -32;
  403. }
  404. if(I != FileI || J != FileJ)
  405. {
  406. printf("Error: Initial Conditions File dimensions not compatible with Parameters File.\n");
  407. return -97;
  408. }
  409. return 0;
  410. }
  411. int loadInitialConditions(const char *InitialConditionsFile, size_t I, size_t J,
  412. double *Vx, double *Vy, double *Phi, double *C_RD, double *C_RT, double *C_G, double *C_Gc,
  413. double *X, double *Y)
  414. {
  415. // Function for loading initial conditions.
  416. int varID_Vx, varID_Vy, varID_Phi, varID_C_RD, varID_C_RT, varID_C_G, varID_C_Gc, varID_X, varID_Y;
  417. int ncID, RetVAL;
  418. // Retrieving initial conditions file ID
  419. ncID = readICNetCDF(InitialConditionsFile);
  420. if (ncID == -30) { return -30; }
  421. // Retrieving variable IDs from initial conditions file
  422. RetVAL = loadICNetCDF(ncID, &varID_Vx, &varID_Vy, &varID_Phi, &varID_C_RD, &varID_C_RT, &varID_C_G, &varID_C_Gc, &varID_X, &varID_Y);
  423. if (RetVAL != 0) { nc_close(ncID); return RetVAL;}
  424. // Checking if dimensions in initial conditions file match those in parameters file
  425. RetVAL = checkICDimensions(ncID,I,J);
  426. if (RetVAL != 0) { nc_close(ncID); return RetVAL; }
  427. // Loading initial conditions matrices
  428. RetVAL = nc_get_var_double(ncID, varID_Vx, Vx);
  429. if (RetVAL != NC_NOERR)
  430. {
  431. printf("Error reading Vx matrix: %s\n", nc_strerror(RetVAL));
  432. nc_close(ncID);
  433. return -33;
  434. }
  435. RetVAL = nc_get_var_double(ncID, varID_Vy, Vy);
  436. if (RetVAL != NC_NOERR)
  437. {
  438. printf("Error reading Vy matrix: %s\n", nc_strerror(RetVAL));
  439. nc_close(ncID);
  440. return -33;
  441. }
  443. RetVAL = nc_get_var_double(ncID, varID_Phi, Phi);
  444. if (RetVAL != NC_NOERR)
  445. {
  446. printf("Error reading Phi matrix: %s\n", nc_strerror(RetVAL));
  447. nc_close(ncID);
  448. return -33;
  449. }
  451. RetVAL = nc_get_var_double(ncID, varID_C_RD, C_RD);
  452. if (RetVAL != NC_NOERR)
  453. {
  454. printf("Error reading C_RD matrix: %s\n", nc_strerror(RetVAL));
  455. nc_close(ncID);
  456. return -33;
  457. }
  459. RetVAL = nc_get_var_double(ncID, varID_C_RT, C_RT);
  460. if (RetVAL != NC_NOERR)
  461. {
  462. printf("Error reading C_RT matrix: %s\n", nc_strerror(RetVAL));
  463. nc_close(ncID);
  464. return -33;
  465. }
  467. RetVAL = nc_get_var_double(ncID, varID_C_G, C_G);
  468. if (RetVAL != NC_NOERR)
  469. {
  470. printf("Error reading C_G matrix: %s\n", nc_strerror(RetVAL));
  471. nc_close(ncID);
  472. return -33;
  473. }
  475. RetVAL = nc_get_var_double(ncID, varID_C_Gc, C_Gc);
  476. if (RetVAL != NC_NOERR)
  477. {
  478. printf("Error reading C_Gc matrix: %s\n", nc_strerror(RetVAL));
  479. nc_close(ncID);
  480. return -33;
  481. }
  483. RetVAL = nc_get_var_double(ncID, varID_X, X);
  484. if (RetVAL != NC_NOERR)
  485. {
  486. printf("Error reading X vector: %s\n", nc_strerror(RetVAL));
  487. nc_close(ncID);
  488. return -33;
  489. }
  491. RetVAL = nc_get_var_double(ncID, varID_Y, Y);
  492. if (RetVAL != NC_NOERR)
  493. {
  494. printf("Error reading Y vector: %s\n", nc_strerror(RetVAL));
  495. nc_close(ncID);
  496. return -33;
  497. }
  499. // Closing the file
  500. RetVAL = nc_close(ncID);
  501. if (RetVAL != NC_NOERR)
  502. {
  503. printf("Error closing file: %s\n", nc_strerror(RetVAL));
  504. return -38;
  505. }
  506. return 0;
  507. }
  508. int loadSolutionNetCDF(int ncID, int *varID_Vx, int *varID_Vy, int *varID_Phi, int *varID_C_RD, int *varID_C_RT, int *varID_C_G, int *varID_C_Gc)
  509. {
  510. // Function for loading NetCDF files for initial condition.
  511. // This function will return variable IDs and check if the needed initial condition matrices are included in the file.
  512. int RetVAL;
  513. RetVAL = nc_inq_varid(ncID, "VsX", varID_Vx);
  514. if (RetVAL != NC_NOERR)
  515. {
  516. printf("Error: Vx matrix not found in the file. %s\n", nc_strerror(RetVAL));
  517. return -31;
  518. }
  519. RetVAL = nc_inq_varid(ncID, "VsY", varID_Vy);
  520. if (RetVAL != NC_NOERR)
  521. {
  522. printf("Error: Vy matrix not found in the file. %s\n", nc_strerror(RetVAL));
  523. return -31;
  524. }
  525. RetVAL = nc_inq_varid(ncID, "Phi", varID_Phi);
  526. if (RetVAL != NC_NOERR)
  527. {
  528. printf("Error: Phi matrix not found in the file. %s\n", nc_strerror(RetVAL));
  529. return -31;
  530. }
  531. RetVAL = nc_inq_varid(ncID, "RD", varID_C_RD);
  532. if (RetVAL != NC_NOERR)
  533. {
  534. printf("Error: C_RD matrix not found in the file. %s\n", nc_strerror(RetVAL));
  535. return -31;
  536. }
  537. RetVAL = nc_inq_varid(ncID, "RT", varID_C_RT);
  538. if (RetVAL != NC_NOERR)
  539. {
  540. printf("Error: C_RT matrix not found in the file. %s\n", nc_strerror(RetVAL));
  541. return -31;
  542. }
  543. RetVAL = nc_inq_varid(ncID, "G", varID_C_G);
  544. if (RetVAL != NC_NOERR)
  545. {
  546. printf("Error: C_G matrix not found in the file. %s\n", nc_strerror(RetVAL));
  547. return -31;
  548. }
  549. RetVAL = nc_inq_varid(ncID, "Gc", varID_C_Gc);
  550. if (RetVAL != NC_NOERR)
  551. {
  552. printf("Error: C_Gc matrix not found in the file. %s\n", nc_strerror(RetVAL));
  553. return -31;
  554. }
  555. return 0;
  556. }
  557. int loadInitialConditionsFromSolution(const char *InitialConditionsFile, const char *SaveFileName, int SolutionIdx,
  558. size_t I, size_t J, double *X, double *Y,
  559. double *Vx, double *Vy, double *Phi, double *C_RD, double *C_RT, double *C_G, double *C_Gc, double *P,
  560. double *Time)
  561. {
  562. // Function for loading initial conditions for solution file.
  563. // In case the simulation was already run and terminated before finish.
  564. // Loading X and Y arrays
  565. int varID_Vx, varID_Vy, varID_Phi, varID_C_RD, varID_C_RT, varID_C_G, varID_C_Gc, varID_X, varID_Y;
  566. int ncID, RetVAL;
  567. // Retrieving initial conditions file ID
  568. ncID = readICNetCDF(InitialConditionsFile);
  569. if (ncID == -30) { return -30; }
  570. // Retrieving variable IDs from initial conditions file
  571. RetVAL = loadICNetCDF(ncID, &varID_Vx, &varID_Vy, &varID_Phi, &varID_C_RD, &varID_C_RT, &varID_C_G, &varID_C_Gc, &varID_X, &varID_Y);
  572. if (RetVAL != 0) { nc_close(ncID); return RetVAL;}
  573. RetVAL = nc_get_var_double(ncID, varID_X, X);
  574. if (RetVAL != NC_NOERR)
  575. {
  576. printf("Error reading X vector: %s\n", nc_strerror(RetVAL));
  577. nc_close(ncID);
  578. return -33;
  579. }
  581. RetVAL = nc_get_var_double(ncID, varID_Y, Y);
  582. if (RetVAL != NC_NOERR)
  583. {
  584. printf("Error reading Y vector: %s\n", nc_strerror(RetVAL));
  585. nc_close(ncID);
  586. return -33;
  587. }
  588. // Closing the file
  589. RetVAL = nc_close(ncID);
  590. if (RetVAL != NC_NOERR)
  591. {
  592. printf("Error closing Initial Conditions file: %s\n", nc_strerror(RetVAL));
  593. return -38;
  594. }
  595. // Clearing variables
  596. ncID, RetVAL = 0;
  597. varID_Vx, varID_Vy, varID_Phi, varID_C_RD, varID_C_RT, varID_C_G, varID_C_Gc = 0;
  598. int varID_Time, varID_P;
  599. // Loading initial conditions from Solution file.
  600. ncID = readICNetCDF(SaveFileName);
  601. // Retrieving variable IDs from Solution file.
  602. RetVAL = loadSolutionNetCDF(ncID, &varID_Vx, &varID_Vy, &varID_Phi, &varID_C_RD, &varID_C_RT, &varID_C_G, &varID_C_Gc);
  603. if (RetVAL != 0) { nc_close(ncID); return RetVAL;}
  604. // Open the solution file
  605. RetVAL = nc_open(SaveFileName, NC_NOWRITE, &ncID);
  606. if (RetVAL != NC_NOERR) {
  607. printf("Error opening solution file: %s\n", nc_strerror(RetVAL));
  608. return -30;
  609. }
  610. // Get variable IDs for required variables
  611. if ((RetVAL = nc_inq_varid(ncID, "VsX", &varID_Vx)) != NC_NOERR ||
  612. (RetVAL = nc_inq_varid(ncID, "VsY", &varID_Vy)) != NC_NOERR ||
  613. (RetVAL = nc_inq_varid(ncID, "Time", &varID_Time)) != NC_NOERR ||
  614. (RetVAL = nc_inq_varid(ncID, "Phi", &varID_Phi)) != NC_NOERR ||
  615. (RetVAL = nc_inq_varid(ncID, "RD", &varID_C_RD)) != NC_NOERR ||
  616. (RetVAL = nc_inq_varid(ncID, "RT", &varID_C_RT)) != NC_NOERR ||
  617. (RetVAL = nc_inq_varid(ncID, "G", &varID_C_G)) != NC_NOERR ||
  618. (RetVAL = nc_inq_varid(ncID, "Gc", &varID_C_Gc)) != NC_NOERR ||
  619. (RetVAL = nc_inq_varid(ncID, "P", &varID_P)) != NC_NOERR) {
  620. printf("Error: Required variable not found in file: %s\n", nc_strerror(RetVAL));
  621. nc_close(ncID);
  622. return -31;
  623. }
  624. // Define start and count for slicing data
  625. size_t start3D[3] = {SolutionIdx, 0, 0};
  626. size_t count3D_VsX[3] = {1, I - 1, J}; // For VsX
  627. size_t count3D_VsY[3] = {1, I, J - 1}; // For VsY
  628. size_t count3D[3] = {1, I, J}; // For 3D variables (Phi, RD, RT, G, Gc, P)
  629. size_t start1D[1] = {SolutionIdx};
  630. size_t count1D[1] = {1};
  631. // Read the variables
  632. if ((RetVAL = nc_get_vara_double(ncID, varID_Vx, start3D, count3D_VsX, Vx)) != NC_NOERR ||
  633. (RetVAL = nc_get_vara_double(ncID, varID_Vy, start3D, count3D_VsY, Vy)) != NC_NOERR ||
  634. (RetVAL = nc_get_vara_double(ncID, varID_Phi, start3D, count3D, Phi)) != NC_NOERR ||
  635. (RetVAL = nc_get_vara_double(ncID, varID_C_RD, start3D, count3D, C_RD)) != NC_NOERR ||
  636. (RetVAL = nc_get_vara_double(ncID, varID_C_RT, start3D, count3D, C_RT)) != NC_NOERR ||
  637. (RetVAL = nc_get_vara_double(ncID, varID_C_G, start3D, count3D, C_G)) != NC_NOERR ||
  638. (RetVAL = nc_get_vara_double(ncID, varID_C_Gc, start3D, count3D, C_Gc)) != NC_NOERR ||
  639. (RetVAL = nc_get_vara_double(ncID, varID_P, start3D, count3D, P)) != NC_NOERR) {
  640. printf("Error reading variable data: %s\n", nc_strerror(RetVAL));
  641. nc_close(ncID);
  642. return -33;
  643. }
  644. // Read the Time variable
  645. if ((RetVAL = nc_get_vara_double(ncID, varID_Time, start1D, count1D, Time)) != NC_NOERR) {
  646. printf("Error reading Time data: %s\n", nc_strerror(RetVAL));
  647. nc_close(ncID);
  648. return -33;
  649. }
  650. // Close the file
  651. if ((RetVAL = nc_close(ncID)) != NC_NOERR) {
  652. printf("Error closing file: %s\n", nc_strerror(RetVAL));
  653. return -38;
  654. }
  655. return 0;
  656. }
  657. int CreateSolutionsFile(const char *SaveFileName, size_t I, size_t J, int MaxTimeSteps)
  658. {
  659. // Function for creating solutions file.
  660. // Space needed is allocated before the start of simulation.
  661. // The empty matrix values are stored as biggest double value and not zero.
  662. int ncID, RetVAL;
  663. // Creating Solutions file
  664. int timeDimID, xDimID, yDimID, xDimVsXID, yDimVsYID, varID_Time, varID_VnsX, varID_VnsY, varID_VsX, varID_VsY, varID_VsX_1, varID_VsY_1, varID_Phi, varID_Phi_1, varID_RD, varID_RT, varID_G, varID_Gc, varID_P, varID_SimTime;
  665. RetVAL = nc_create(SaveFileName, NC_NETCDF4 | NC_CLOBBER, &ncID);
  666. if (RetVAL != NC_NOERR)
  667. {
  668. printf("Error creating file: %s\n", nc_strerror(RetVAL));
  669. return 34;
  670. }
  671. // Define dimensions
  672. RetVAL = nc_def_dim(ncID, "time", MaxTimeSteps, &timeDimID);
  673. if (RetVAL != NC_NOERR)
  674. {
  675. printf("Error defining time dimension: %s\n", nc_strerror(RetVAL));
  676. return 35;
  677. }
  678. RetVAL = nc_def_dim(ncID, "xDim", I, &xDimID);
  679. if (RetVAL != NC_NOERR)
  680. {
  681. printf("Error defining x dimension: %s\n", nc_strerror(RetVAL));
  682. return 35;
  683. }
  684. RetVAL = nc_def_dim(ncID, "yDim", J, &yDimID);
  685. if (RetVAL != NC_NOERR)
  686. {
  687. printf("Error defining y dimension: %s\n", nc_strerror(RetVAL));
  688. return 35;
  689. }
  690. // Define additional dimensions for VsX and VsY
  691. RetVAL = nc_def_dim(ncID, "xDimS", I - 1, &xDimVsXID);
  692. if (RetVAL != NC_NOERR)
  693. {
  694. printf("Error defining xDimS dimension: %s\n", nc_strerror(RetVAL));
  695. return 35;
  696. }
  698. RetVAL = nc_def_dim(ncID, "yDimS", J - 1, &yDimVsYID);
  699. if (RetVAL != NC_NOERR)
  700. {
  701. printf("Error defining yDimS dimension: %s\n", nc_strerror(RetVAL));
  702. return 35;
  703. }
  704. // Define variables
  705. int dims3D[3] = {timeDimID, xDimID, yDimID};
  706. int dimsTime[1] = {timeDimID};
  707. int dims2D[2] = {timeDimID, 1};
  708. // Time variable
  709. RetVAL = nc_def_var(ncID, "Time", NC_DOUBLE, 1, dimsTime, &varID_Time);
  710. if (RetVAL != NC_NOERR)
  711. {
  712. printf("Error defining Time variable: %s\n", nc_strerror(RetVAL));
  713. return 36;
  714. }
  715. // Simulation Time variable
  716. RetVAL = nc_def_var(ncID, "SimulationTime", NC_DOUBLE, 2, dims2D, &varID_SimTime);
  717. if (RetVAL != NC_NOERR)
  718. {
  719. printf("Error defining SimulationTime variable: %s\n", nc_strerror(RetVAL));
  720. return 36;
  721. }
  722. // Solution matrices
  723. RetVAL = nc_def_var(ncID, "VnsX", NC_DOUBLE, 3, dims3D, &varID_VnsX);
  724. if (RetVAL != NC_NOERR)
  725. {
  726. printf("Error defining VnsX variable: %s\n", nc_strerror(RetVAL));
  727. return 36;
  728. }
  729. RetVAL = nc_def_var(ncID, "VnsY", NC_DOUBLE, 3, dims3D, &varID_VnsY);
  730. if (RetVAL != NC_NOERR)
  731. {
  732. printf("Error defining VnsY variable: %s\n", nc_strerror(RetVAL));
  733. return 36;
  734. }
  735. // Redefine the dimensions for VsX (3D array with I-1 rows)
  736. dims3D[1] = xDimVsXID;
  737. RetVAL = nc_def_var(ncID, "VsX", NC_DOUBLE, 3, dims3D, &varID_VsX);
  738. if (RetVAL != NC_NOERR)
  739. {
  740. printf("Error defining VsX variable: %s\n", nc_strerror(RetVAL));
  741. return 36;
  742. }
  743. RetVAL = nc_def_var(ncID, "VsX_1", NC_DOUBLE, 3, dims3D, &varID_VsX_1);
  744. if (RetVAL != NC_NOERR)
  745. {
  746. printf("Error defining VsX_1 variable: %s\n", nc_strerror(RetVAL));
  747. return 36;
  748. }
  749. // Redefine the dimensions for VsY (3D array with J-1 columns)
  750. dims3D[1] = xDimID; // Reset to xDim for VsY
  751. dims3D[2] = yDimVsYID;
  752. RetVAL = nc_def_var(ncID, "VsY", NC_DOUBLE, 3, dims3D, &varID_VsY);
  753. if (RetVAL != NC_NOERR)
  754. {
  755. printf("Error defining VsY variable: %s\n", nc_strerror(RetVAL));
  756. return 36;
  757. }
  758. RetVAL = nc_def_var(ncID, "VsY_1", NC_DOUBLE, 3, dims3D, &varID_VsY_1);
  759. if (RetVAL != NC_NOERR)
  760. {
  761. printf("Error defining VsY_1 variable: %s\n", nc_strerror(RetVAL));
  762. return 36;
  763. }
  764. // Reset to original xDim and yDim for remaining matrices
  765. dims3D[1] = xDimID;
  766. dims3D[2] = yDimID;
  767. RetVAL = nc_def_var(ncID, "Phi", NC_DOUBLE, 3, dims3D, &varID_Phi);
  768. if (RetVAL != NC_NOERR)
  769. {
  770. printf("Error defining Phi variable: %s\n", nc_strerror(RetVAL));
  771. return 36;
  772. }
  773. RetVAL = nc_def_var(ncID, "Phi_1", NC_DOUBLE, 3, dims3D, &varID_Phi_1);
  774. if (RetVAL != NC_NOERR)
  775. {
  776. printf("Error defining Phi_1 variable: %s\n", nc_strerror(RetVAL));
  777. return 36;
  778. }
  779. RetVAL = nc_def_var(ncID, "RD", NC_DOUBLE, 3, dims3D, &varID_RD);
  780. if (RetVAL != NC_NOERR)
  781. {
  782. printf("Error defining RD variable: %s\n", nc_strerror(RetVAL));
  783. return 36;
  784. }
  785. RetVAL = nc_def_var(ncID, "RT", NC_DOUBLE, 3, dims3D, &varID_RT);
  786. if (RetVAL != NC_NOERR)
  787. {
  788. printf("Error defining RT variable: %s\n", nc_strerror(RetVAL));
  789. return 36;
  790. }
  791. RetVAL = nc_def_var(ncID, "G", NC_DOUBLE, 3, dims3D, &varID_G);
  792. if (RetVAL != NC_NOERR)
  793. {
  794. printf("Error defining G variable: %s\n", nc_strerror(RetVAL));
  795. return 36;
  796. }
  797. RetVAL = nc_def_var(ncID, "Gc", NC_DOUBLE, 3, dims3D, &varID_Gc);
  798. if (RetVAL != NC_NOERR)
  799. {
  800. printf("Error defining Gc variable: %s\n", nc_strerror(RetVAL));
  801. return 36;
  802. }
  803. RetVAL = nc_def_var(ncID, "P", NC_DOUBLE, 3, dims3D, &varID_P);
  804. if (RetVAL != NC_NOERR)
  805. {
  806. printf("Error defining P variable: %s\n", nc_strerror(RetVAL));
  807. return 36;
  808. }
  809. // End definitions and leave define mode
  810. RetVAL = nc_enddef(ncID);
  811. if (RetVAL != NC_NOERR)
  812. {
  813. printf("Error leaving define mode: %s\n", nc_strerror(RetVAL));
  814. return 37;
  815. }
  816. // Close the file (it will be reopened for writing during the simulation)
  817. RetVAL = nc_close(ncID);
  818. if (RetVAL != NC_NOERR)
  819. {
  820. printf("Error closing file: %s\n", nc_strerror(RetVAL));
  821. return 38;
  822. }
  823. return 0;
  824. }
  825. void GenerateTracerPosition(int Cluster, int Seed, int NumPoints, double Rmin, double Rmax, double DomainLength, double *VecX, double *VecY)
  826. {
  827. // Function for generating initial positions of Tracers.
  828. // Rmin and Rmax parameters define the minimum and maximum radius the point can have with respect to the domain center.
  829. double x,y,r,phi;
  830. // The seed is defined in Parameters file.
  831. srand(Seed);
  832. for (int i = 0; i < NumPoints; i++)
  833. {
  834. double u = (double)rand() / RAND_MAX;
  835. // the radius is proportional to sqrt(u) so that we get homogeneous distribution in x and y dimensions.
  836. r = (Rmax - Rmin) * sqrt(u) + Rmin;
  837. phi = 2 * M_PI * ((double)rand() / RAND_MAX);
  838. x = DomainLength/2 + r * cos(phi);
  839. y = DomainLength/2 + r * sin(phi);
  840. if(Cluster==0 || NumPoints != 0)
  841. {
  842. printf("x:%lf\ty:%lf\n",x,y);
  843. }
  844. VecX[i]=x;
  845. VecY[i]=y;
  846. }
  847. }
  848. int CrateTracerFiles(char *TracerFolder, int TracerNBulk, int TracerNOut,
  849. double *TracerPositionBx, double *TracerPositionBy, double *TracerPositionOx, double *TracerPositionOy,
  850. char TracersBulkFileNames[][100], char TracersOutFileNames[][100])
  851. {
  852. // Function for creating binary files tracer positions will be stored in.
  853. double time_stamp = 0.0;
  854. // Bulk Tracer Files
  855. for (int i = 0; i < TracerNBulk; i++)
  856. {
  857. snprintf(TracersBulkFileNames[i], 100, "%sBulkTracer_%d.bin", TracerFolder, i+1);
  858. FILE *file = fopen(TracersBulkFileNames[i], "wb");
  859. if (!file)
  860. {
  861. printf("Error opening file");
  862. return 40;
  863. }
  864. fwrite(&time_stamp, sizeof(double), 1, file); // This is time stamp equal to 0 for initial position.
  865. fwrite(&TracerPositionBx[i], sizeof(double), 1, file);
  866. fwrite(&TracerPositionBy[i], sizeof(double), 1, file);
  867. fclose(file);
  868. }
  869. // Tracers Out Files
  870. for (int i = 0; i < TracerNOut; i++)
  871. {
  872. snprintf(TracersOutFileNames[i], 100, "%sOutTracer_%d.bin", TracerFolder, i+1);
  873. FILE *file = fopen(TracersOutFileNames[i], "wb");
  874. if (!file)
  875. {
  876. printf("Error opening file");
  877. return 40;
  878. }
  879. fwrite(&time_stamp, sizeof(double), 1, file);
  880. fwrite(&TracerPositionOx[i], sizeof(double), 1, file);
  881. fwrite(&TracerPositionOy[i], sizeof(double), 1, file);
  882. fclose(file);
  883. }
  884. return 0;
  885. }
  886. void printParameters(const char *InitialConditionsFile, const char *SaveFileName,
  887. int nTracersB, int nTracersO,
  888. int NumberOfSaveCycles, int NumberOfDoubleStep, int MaxIter,
  889. size_t I, size_t J,
  890. double dx, double dy, double dt, double CO, double vBB, double vBS,
  891. double PTol, double eps, double nu, double cG, double Si0, double kSi, double kC, double Beta, double Alpha,
  892. double DRD, double DRT, double DG, double DGc, double k1, double k2, double k3, double k11, double k12, double R_Max)
  893. {
  894. // Function for displaying loaded Parameters
  895. printf("Original Initial conditions File is: '%s'.\n", InitialConditionsFile);
  896. printf("File for saving data with the name: '%s', will be created.\n\n", SaveFileName);
  898. printf("Solution Matrices will be saved every: %d steps, which is equivalent to saving every %2.3lf s.\n\n", 2*NumberOfDoubleStep, 2*NumberOfDoubleStep*dt);
  900. printf("Domain is split onto: %zu Points in x and %zu Points in y direction,\n", I, J);
  901. printf("With Grid spacing of: %1.3lf in x and %1.3lf in y direction,\n", dx, dy);
  902. printf("Which gives: %2.3lf by %2.3lf Domain.\n\n", dx*(I-1), dy*(J-1));
  904. printf("Time Step length is: %lf.\n\n", dt);
  906. printf("Gradient CutOff is: %1.1e, while Potential Saturation values are: %2.1lf for Bulk and %2.1lf for Surface.\n\n", CO, vBB, vBS);
  908. printf("Pressure Tolerance at which iterative process is stopped: %1.1e.\n\n", PTol);
  910. printf("Viscosity of the fluid is: %2.3lf, while dissipation coefficient is %2.3lf.\n\n", nu, Beta);
  912. printf("Level-Set Parameters are:\n \t Epsilon: %2.3lf,\n \t Gamma: %2.3lf. \n \t Alpha: %2.3lf.\n\n", eps, cG, Alpha);
  914. printf("Parameters for Force on the membrane:\n \t Surface tension Sigma: %2.3lf,\n \t Surface tension change rate: %2.3lf,\n \t Force from interaction with the substrate coefficient: %2.3lf.\n\n", Si0, kSi, kC);
  916. printf("Diffusion constants are:\n \t D_RD: %2.3lf,\n \t D_RT: %2.3lf,\n \t D_G: %2.3lf,\n \t D_Gc: %2.3lf.\n\n", DRD, DRT, DG, DGc);
  918. printf("Reaction constants are:\n \t k1: %2.3lf,\n \t k2: %2.3lf,\n \t k3: %2.3lf,\n \t k11: %2.3lf,\n \t k12: %2.3lf.\n", k1, k2, k3, k11, k12);
  919. printf("Concentration of binding spaces for Rac1T is: %3.1lf\n\n", R_Max);
  921. printf("Number of tracers inside cell: %d\n", nTracersB);
  922. printf("Number of tracers outside cell: %d\n\n", nTracersO);
  923. }
  925. // Error codes
  926. // 10 - Failed to load or parse the JSON file. (Check input file names and LoadJSON function)
  927. // 11 - Failed to find key in .json file
  928. // 12 - Array size in .json file is missing values or has extra values.
  929. // 20 - Cluster variable not set properly. Can't determine if the simulation is to be run on cluster without inputs or on personal computer with input and print options.
  930. // 30 - Failed to open .nc file
  931. // 31 - Variable not found in .nc file
  932. // 32 - Failed to get the dimension of variable
  933. // 33 - Failed to read the variable
  934. // 34 - Failed to create .nc file
  935. // 35 - Failed to define dimensions in .nc file
  936. // 36 - Failed to define variable in .nc file
  937. // 37 - Failed to exit variable definition mode
  938. // 38 - Failed to close .nc file
  939. // 40 - Failed to open Tracer files
  940. // 95 - Simulation aborted by the user
  941. // 96 - Error occurred during simulation.
  942. // 97 - Initial Conditions File and Parameters File dimensions are not compatible
  943. // 98 - Memory allocation failed.
  944. // 99 - The y/n question input is neither "y" nor "n"
  945. int main(int argc, char *argv[])
  946. {
  947. // Main function for simulating cell locomotion
  948. // Error code
  949. int ERROR;
  950. // y/n input
  951. char InputYN[3];
  952. // Simulation Parameters allocation
  953. char ParametersFile[100];
  954. // File
  955. char InitialConditionsFile[100], SaveFileName[100], TracerFolder[100], TracerPrefix[100];
  956. // Simulation
  957. int Seed, Cluster, nTracersB, nTracersO, NumberOfSaveCycles, NumberOfDoubleStep, MaxIter;
  958. size_t I, J;
  959. double dx, dy, dt, CO, vBB, vBS, PTol;
  960. // Model
  961. double eps, nu, cG, Si0, kSi, kC, Beta, Alpha;
  962. double DRD, DRT, DG, DGc;
  963. double k1, k2, k3, k11, k12;
  964. double R_Max;
  965. // Loading Parameters from .json file
  966. if (argc>=3) // If Parameters file and initial conditions file were provided
  967. {
  968. // Copying provided file names
  969. strcpy(ParametersFile, argv[1]);
  970. strcpy(InitialConditionsFile, argv[2]);
  971. // Loading parameters file
  972. cJSON *JSONFile = loadJSON(ParametersFile);
  973. if (!JSONFile)
  974. {
  975. printf("Error: Failed to load or parse the JSON file\n");
  976. cJSON_Delete(JSONFile); // Clean up cJSON object
  977. return 10;
  978. }
  979. // Loading parameters form .json
  980. ERROR = loadParameters(JSONFile,
  981. InitialConditionsFile, SaveFileName, TracerFolder, TracerPrefix,
  982. &Seed, &Cluster, &nTracersB, &nTracersO, &NumberOfSaveCycles, &NumberOfDoubleStep,
  983. &I, &J,
  984. &dx, &dy, &dt, &CO, &vBB, &vBS, &MaxIter, &PTol,
  985. &eps, &nu, &cG, &Si0, &kSi, &kC, &Beta, &Alpha,
  986. &DRD, &DRT, &DG, &DGc,
  987. &k1, &k2, &k3, &k11, &k12, &R_Max);
  988. if (ERROR != 0)
  989. {
  990. printf("Error: Failed to save parameter values from the JSON file\n");
  991. cJSON_Delete(JSONFile); // Clean up cJSON object
  992. return ERROR;
  993. }
  994. cJSON_Delete(JSONFile); // Clean up cJSON object
  995. }
  996. else // If Parameters file and initial conditions file were both NOT provided
  997. {
  998. // Declaration of variables
  999. char InitialConditionsFileInput[100];
  1000. // Starting simulation preparation
  1001. printf("Cell Locomotion Simulation\n\n");
  1002. // Asking to provide the parameters file name
  1003. printf("Please enter the Parameters File Name: ");
  1004. fgets(ParametersFile, sizeof(ParametersFile), stdin);
  1005. ParametersFile[strcspn(ParametersFile, "\n")] = 0;
  1006. // Loading parameters file
  1007. cJSON *JSONFile = loadJSON(ParametersFile);
  1008. if (!JSONFile)
  1009. {
  1010. printf("Error: Failed to load or parse the JSON file\n");
  1011. cJSON_Delete(JSONFile); // Clean up cJSON object
  1012. return 10;
  1013. }
  1014. // Loading parameters form .json
  1015. ERROR = loadParameters(JSONFile,
  1016. InitialConditionsFile, SaveFileName, TracerFolder, TracerPrefix,
  1017. &Seed, &Cluster, &nTracersB, &nTracersO, &NumberOfSaveCycles, &NumberOfDoubleStep,
  1018. &I, &J,
  1019. &dx, &dy, &dt, &CO, &vBB, &vBS, &MaxIter, &PTol,
  1020. &eps, &nu, &cG, &Si0, &kSi, &kC, &Beta, &Alpha,
  1021. &DRD, &DRT, &DG, &DGc,
  1022. &k1, &k2, &k3, &k11, &k12, &R_Max);
  1023. if (ERROR != 0)
  1024. {
  1025. printf("Error: Failed to load parameter values from the JSON file\n");
  1026. cJSON_Delete(JSONFile); // Clean up cJSON object
  1027. return ERROR;
  1028. }
  1029. cJSON_Delete(JSONFile); // Clean up cJSON object
  1030. // Asking to provide the initial conditions file name
  1031. printf("Do you want to use original Initial Conditions from Parameters File?: (y/n) \n");
  1032. fgets(InputYN, sizeof(InputYN), stdin);
  1033. InputYN[strcspn(InputYN, "\n")] = 0;
  1034. if(strcmp(InputYN,"n")==0)
  1035. {
  1036. // Asking to provide the initial conditions file name
  1037. printf("\nPlease provide the new Initial Conditions file in the same folder as the program before continuing!");
  1038. printf("\nPlease enter the Initial Conditions File Name: \n");
  1039. fgets(InitialConditionsFileInput, sizeof(InitialConditionsFileInput), stdin);
  1040. InitialConditionsFileInput[strcspn(InitialConditionsFileInput, "\n")] = 0;
  1041. // Copying initial conditions file name
  1042. for(size_t i=0; i<100; i++)
  1043. {
  1044. InitialConditionsFile[i]=InitialConditionsFileInput[i];
  1045. }
  1046. }
  1047. else if(strcmp(InputYN,"y")==0)
  1048. {
  1049. // Initial conditions file name was loaded from parameters file and not overwritten
  1050. printf("\nOriginal Initial Conditions will be loaded\n");
  1051. }
  1052. else
  1053. {
  1054. printf("\nError: Please enter y for YES and n for NO.\n");
  1055. return 99;
  1056. }
  1057. }
  1058. // Memory allocation for initial conditions arrays
  1059. ALLOC_ARRAY(Vx,double,I-1,J);
  1060. ALLOC_ARRAY(Vy,double,I,J-1);
  1061. ALLOC_ARRAY(Phi,double,I,J);
  1062. ALLOC_ARRAY(C_RD,double,I,J);
  1063. ALLOC_ARRAY(C_RT,double,I,J);
  1064. ALLOC_ARRAY(C_G,double,I,J);
  1065. ALLOC_ARRAY(C_Gc,double,I,J);
  1066. ALLOC_ARRAY(X,double,I,1);
  1067. ALLOC_ARRAY(Y,double,1,J);
  1068. ALLOC_ARRAY(P,double,I,J);
  1069. if (Vx == NULL) {printf("Error: Memory allocation failed for Vx\n"); return 98;}
  1070. if (Vy == NULL) {printf("Error: Memory allocation failed for Vy\n"); return 98;}
  1071. if (Phi == NULL) {printf("Error: Memory allocation failed for Phi\n"); return 98;}
  1072. if (C_RD == NULL) {printf("Error: Memory allocation failed for C_RD\n"); return 98;}
  1073. if (C_RT == NULL) {printf("Error: Memory allocation failed for C_RT\n"); return 98;}
  1074. if (C_G == NULL) {printf("Error: Memory allocation failed for C_G\n"); return 98;}
  1075. if (C_Gc == NULL) {printf("Error: Memory allocation failed for C_Gc\n"); return 98;}
  1076. if (X == NULL) {printf("Error: Memory allocation failed for X\n"); return 98;}
  1077. if (Y == NULL) {printf("Error: Memory allocation failed for Y\n"); return 98;}
  1078. if (P == NULL) {printf("Error: Memory allocation failed for P\n"); return 98;}
  1079. double Time = 0;
  1080. // Returning error value for loading initial conditions
  1081. int RetVAL;
  1082. // Memory allocation for tracers
  1083. double TracerPositionBx[nTracersB];
  1084. double TracerPositionBy[nTracersB];
  1085. double TracerPositionOx[nTracersO];
  1086. double TracerPositionOy[nTracersO];
  1087. char TracersBulkFileNames[nTracersB][100];
  1088. char TracersOutFileNames[nTracersO][100];
  1089. // Memory allocation for pressuer iteration number array
  1090. ALLOC_ARRAY(nIterStepVec, int, NumberOfSaveCycles*NumberOfDoubleStep*2+1, 1);
  1091. // Initiating something
  1092. double SimulationTime = 0;
  1093. double PrctSave=2;
  1094. // Forking for cluster or personal computer simulation setup
  1095. if(Cluster==1)
  1096. {
  1097. // Checking if the solution file already exists which signals the simulation has been terminated before
  1098. int SolutionIdx = checkSolutionFile(SaveFileName);
  1099. if (SolutionIdx < -1)
  1100. {
  1101. printf("Error: Failed to access Solutions file!\n");
  1102. if (Vx) free(Vx);
  1103. if (Vy) free(Vy);
  1104. if (Phi) free(Phi);
  1105. if (C_RD) free(C_RD);
  1106. if (C_RT) free(C_RT);
  1107. if (C_G) free(C_G);
  1108. if (C_Gc) free(C_Gc);
  1109. if (X) free(X);
  1110. if (Y) free(Y);
  1111. if (P) free(P);
  1112. if (nIterStepVec) free(nIterStepVec);
  1113. return -SolutionIdx;
  1114. }
  1115. else if (SolutionIdx == -1)
  1116. {
  1117. // Solution file does not exist
  1118. RetVAL = loadInitialConditions(InitialConditionsFile,I,J,Vx,Vy,Phi,C_RD,C_RT,C_G,C_Gc,X,Y);
  1119. if (RetVAL < 0)
  1120. {
  1121. return -RetVAL;
  1122. }
  1123. CreateSolutionsFile(SaveFileName,I,J,NumberOfSaveCycles+1);
  1124. }
  1125. else
  1126. {
  1127. // Solution file exists
  1128. loadInitialConditionsFromSolution(InitialConditionsFile,SaveFileName,SolutionIdx,
  1129. I,J,X,Y,
  1130. Vx,Vy,Phi,C_RD,C_RT,C_G,C_Gc,P,
  1131. &Time);
  1132. }
  1133. GenerateTracerPosition(Cluster,Seed,nTracersB,0,6.21,X[I-1],TracerPositionBx,TracerPositionBy);
  1134. GenerateTracerPosition(Cluster,Seed,nTracersO,6.5,11.2,X[I-1],TracerPositionOx,TracerPositionOy);
  1135. RetVAL = CrateTracerFiles(TracerPrefix,nTracersB,nTracersO,TracerPositionBx,TracerPositionBy,TracerPositionOx,TracerPositionOy,TracersBulkFileNames,TracersOutFileNames);
  1136. if (RetVAL != 0)
  1137. {
  1138. printf("Error: Failed to create Tracer files!\n");
  1139. if (Vx) free(Vx);
  1140. if (Vy) free(Vy);
  1141. if (Phi) free(Phi);
  1142. if (C_RD) free(C_RD);
  1143. if (C_RT) free(C_RT);
  1144. if (C_G) free(C_G);
  1145. if (C_Gc) free(C_Gc);
  1146. if (X) free(X);
  1147. if (Y) free(Y);
  1148. if (P) free(P);
  1149. if (nIterStepVec) free(nIterStepVec);
  1150. return RetVAL;
  1151. }
  1152. SimulationTime = simulation(SolutionIdx,Time,Cluster,SaveFileName,NumberOfSaveCycles,NumberOfDoubleStep,PrctSave,
  1153. nTracersB,nTracersO,TracerPositionBx,TracerPositionBy,TracerPositionOx,TracerPositionOy,TracersBulkFileNames,TracersOutFileNames,
  1154. I,J,dx,dy,dt,eps,nu,cG,
  1155. CO,vBB,vBS,MaxIter,PTol,nIterStepVec,
  1156. Si0,kSi,kC,Beta,Alpha,
  1157. DRD,DRT,DG,DGc,k1,k2,k3,k11,k12,R_Max,
  1158. P,Vx,Vy,Phi,C_RT,C_RD,C_G,C_Gc,X,Y);
  1159. if(SimulationTime==-1)
  1160. {
  1161. printf("Error occurred closing tracer files.\n");
  1162. if (Vx) free(Vx);
  1163. if (Vy) free(Vy);
  1164. if (Phi) free(Phi);
  1165. if (C_RD) free(C_RD);
  1166. if (C_RT) free(C_RT);
  1167. if (C_G) free(C_G);
  1168. if (C_Gc) free(C_Gc);
  1169. if (X) free(X);
  1170. if (Y) free(Y);
  1171. if (nIterStepVec) free(nIterStepVec);
  1172. return 96;
  1173. }
  1174. }
  1175. else if(Cluster==0)
  1176. {
  1177. printf("Parameters file loaded!\n\n");
  1178. printParameters(InitialConditionsFile,SaveFileName,
  1179. nTracersB,nTracersO,
  1180. NumberOfSaveCycles,NumberOfDoubleStep,MaxIter,
  1181. I, J,
  1182. dx,dy,dt,CO,vBB,vBS,
  1183. PTol,eps,nu,cG,Si0,kSi,kC,Beta, Alpha,
  1184. DRD,DRT,DG,DGc,
  1185. k1,k2,k3,k11,k12,R_Max);
  1186. printf("Checking existence of Solutions file!\n");
  1187. int SolutionIdx = checkSolutionFile(SaveFileName);
  1188. if (SolutionIdx < -1)
  1189. {
  1190. printf("Error: Failed to access Solutions file!\n");
  1191. if (Vx) free(Vx);
  1192. if (Vy) free(Vy);
  1193. if (Phi) free(Phi);
  1194. if (C_RD) free(C_RD);
  1195. if (C_RT) free(C_RT);
  1196. if (C_G) free(C_G);
  1197. if (C_Gc) free(C_Gc);
  1198. if (X) free(X);
  1199. if (Y) free(Y);
  1200. if (P) free(P);
  1201. if (nIterStepVec) free(nIterStepVec);
  1202. return -SolutionIdx;
  1203. }
  1204. else if (SolutionIdx == -1)
  1205. {
  1206. // Solution file does not exist
  1207. printf("Solution file does not exist!\n");
  1208. printf("Initial conditions will be loaded from Initial Conditions file!\n");
  1209. RetVAL = loadInitialConditions(InitialConditionsFile,I,J,Vx,Vy,Phi,C_RD,C_RT,C_G,C_Gc,X,Y);
  1210. if (RetVAL < 0)
  1211. {
  1212. return -RetVAL;
  1213. }
  1214. printf("Solutions file will be created\n");
  1215. CreateSolutionsFile(SaveFileName,I,J,NumberOfSaveCycles+1);
  1216. }
  1217. else
  1218. {
  1219. // Solution file exists
  1220. printf("Solution file exists!\n");
  1221. printf("Would you like to load initial conditions from Solutions File and continue the simulation?: (y/n) \n");
  1222. fgets(InputYN, sizeof(InputYN), stdin);
  1223. InputYN[strcspn(InputYN, "\n")] = 0;
  1224. if(strcmp(InputYN,"n")==0)
  1225. {
  1226. printf("Initial conditions will be loaded from Initial Conditions file!\n");
  1227. RetVAL = loadInitialConditions(InitialConditionsFile,I,J,Vx,Vy,Phi,C_RD,C_RT,C_G,C_Gc,X,Y);
  1228. if (RetVAL < 0)
  1229. {
  1230. return -RetVAL;
  1231. }
  1232. printf("Solutions file will be created\n");
  1233. CreateSolutionsFile(SaveFileName,I,J,NumberOfSaveCycles+1);
  1234. }
  1235. else if(strcmp(InputYN,"y")==0)
  1236. {
  1237. printf("Initial conditions will be loaded from Solutions file!\n");
  1238. printf("Solution index is = %d\n",SolutionIdx);
  1239. loadInitialConditionsFromSolution(InitialConditionsFile,SaveFileName,SolutionIdx,
  1240. I,J,X,Y,
  1241. Vx,Vy,Phi,C_RD,C_RT,C_G,C_Gc,P,
  1242. &Time);
  1243. }
  1244. else
  1245. {
  1246. printf("\nError: Please enter y for YES and n for NO.\n");
  1247. return 99;
  1248. }
  1249. }
  1250. if (nTracersB != 0) {printf("Tracer Positions in the Bulk:\n");}
  1251. GenerateTracerPosition(Cluster,Seed,nTracersB,0,6.21,X[I-1],TracerPositionBx,TracerPositionBy);
  1252. if (nTracersO != 0) {printf("Tracer Positions Outside the cell:\n");}
  1253. GenerateTracerPosition(Cluster,Seed,nTracersO,6.5,11.2,X[I-1],TracerPositionOx,TracerPositionOy);
  1254. if (nTracersO != 0 || nTracersB != 0) {printf("\nTracer Files will be created in: %s\n",TracerFolder);}
  1255. CrateTracerFiles(TracerFolder,nTracersB,nTracersO,TracerPositionBx,TracerPositionBy,TracerPositionOx,TracerPositionOy,TracersBulkFileNames,TracersOutFileNames);
  1256. printf("\n\nStart simulation? ((y/n))\n");
  1257. fgets(InputYN, sizeof(InputYN), stdin);
  1258. InputYN[strcspn(InputYN, "\n")] = 0;
  1259. if(strcmp(InputYN,"n")==0)
  1260. {
  1261. printf("Simulation aborted...\n");
  1262. if (Vx) free(Vx);
  1263. if (Vy) free(Vy);
  1264. if (Phi) free(Phi);
  1265. if (C_RD) free(C_RD);
  1266. if (C_RT) free(C_RT);
  1267. if (C_G) free(C_G);
  1268. if (C_Gc) free(C_Gc);
  1269. if (X) free(X);
  1270. if (Y) free(Y);
  1271. if (P) free(P);
  1272. if (nIterStepVec) free(nIterStepVec);
  1273. return 95;
  1274. }
  1275. else if(strcmp(InputYN,"y")==0)
  1276. {
  1277. SimulationTime = simulation(SolutionIdx,Time,Cluster,SaveFileName,NumberOfSaveCycles,NumberOfDoubleStep,PrctSave,
  1278. nTracersB,nTracersO,TracerPositionBx,TracerPositionBy,TracerPositionOx,TracerPositionOy,TracersBulkFileNames,TracersOutFileNames,
  1279. I,J,dx,dy,dt,eps,nu,cG,
  1280. CO,vBB,vBS,MaxIter,PTol,nIterStepVec,
  1281. Si0,kSi,kC,Beta,Alpha,
  1282. DRD,DRT,DG,DGc,k1,k2,k3,k11,k12,R_Max,
  1283. P,Vx,Vy,Phi,C_RT,C_RD,C_G,C_Gc,X,Y);
  1284. if(SimulationTime==-1)
  1285. {
  1286. printf("Error occurred closing tracer files.\n");
  1287. if (Vx) free(Vx);
  1288. if (Vy) free(Vy);
  1289. if (Phi) free(Phi);
  1290. if (C_RD) free(C_RD);
  1291. if (C_RT) free(C_RT);
  1292. if (C_G) free(C_G);
  1293. if (C_Gc) free(C_Gc);
  1294. if (X) free(X);
  1295. if (Y) free(Y);
  1296. if (P) free(P);
  1297. if (nIterStepVec) free(nIterStepVec);
  1298. return 96;
  1299. }
  1300. printf("Simulation took %lf s.", SimulationTime);
  1301. char ch;
  1302. scanf(" %c", &ch); // Not to close the window. When running on Windows from cmd.
  1303. }
  1304. else
  1305. {
  1306. printf("\nError: Please enter y for YES and n for NO.\n");
  1307. return 99;
  1308. }
  1309. }
  1310. else
  1311. {
  1312. printf("\nError: Value of \"Cluster\" variable neither 0 nor 1. Please check the parameters file!\n");
  1313. return 20;
  1314. }
  1315. if (Vx) free(Vx);
  1316. if (Vy) free(Vy);
  1317. if (Phi) free(Phi);
  1318. if (C_RD) free(C_RD);
  1319. if (C_RT) free(C_RT);
  1320. if (C_G) free(C_G);
  1321. if (C_Gc) free(C_Gc);
  1322. if (X) free(X);
  1323. if (Y) free(Y);
  1324. if (P) free(P);
  1325. if (nIterStepVec) free(nIterStepVec);
  1326. return 0;
  1327. }

main.c at commit c76d59f, under MIT · at the source

Overview

Authors: Blaž Ivšić1,2, Dorijan Vulić1, Igor Weber3, Piotr Nowakowski1, Ana-Sunčana Smith1,4,5
ORCID iDs: Blaž Ivšić, Dorijan Vulić, Piotr Nowakowski
  1. Division of Physical Chemistry, Ruđer Bošković Institute, Zagreb, Croatia
  2. Centre for Advanced Laser Techniques, Institute of Physics, Zagreb, Croatia
  3. Division of Molecular Biology, Ruđer Bošković Institute, Zagreb, Croatia
  4. Department of Physics, Faculty of Sciences, Friedrich-Alexander-Universität, Erlangen, Bavaria, Germany
  5. Competence Center Engineering of Advanced Materials, Friedrich-Alexander-Universität, Erlangen, Bavaria, Germany
Institutions: Ruđer Bošković Institute (Croatia); Institute of Physics (Croatia); Friedrich-Alexander-Universität Erlangen-Nürnberg (Germany)
Journal: PLoS computational biology, volume 22, issue 4, article e1014216
Dates: received 21 November 2025; accepted 9 April 2026; published online 27 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014216 · PMID 42044154 · PMCID PMC13148778 · OpenAlex W4415938307
Open access: gold, a free copy (OpenAlex)
Preprint: arxiv:2511.00558
Status: code verified
Categories: computational modeling (no new data) (modality), none (in silico) (organism), cellular / molecular (subfield)
MeSH: Cell Movement*, Cytosol*, Models, Biological*, Animals, Cell Polarity, Computational Biology, Cytoskeleton, Diffusion, Hydrodynamics, rho GTP-Binding Proteins, Signal Transduction (* major topic)
Journal subjects: Biology and Life Sciences, Cell Biology, Cell Motility, Physical Sciences, Physics, Classical Mechanics, Continuum Mechanics, Fluid Mechanics, Surface Tension, Fluid Dynamics, Fluid Flow, Cell Physiology, Cell Polarity, Biochemistry, Biochemical Simulations, Computational Biology, Cytosol, Chemistry, Physical Chemistry, Reaction Dynamics, Motion, Velocity
Topic: Cellular Mechanics and Interactions (Cell Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: Hrvatska Zaklada za Znanost (HRZZ-IP-2024-05-6331, HRZZ-IP-2025-02-3976); Deutsche Forschungsgemeinschaft (SM 289/10-1, SM 289/11-1)
Citations: not cited yet (Europe PMC); 126 references in the paper

Abstract

Cell migration is a fundamental process underlying the survival and function of both unicellular and multicellular organisms. Crawling motility in eukaryotic cells arises from cyclic protrusion and retraction driven by the cytoskeleton, whose organization is regulated by reaction–diffusion (RD) dynamics of Rho GTPases between the cytosol and the cortex. These dynamics generate spatial membrane patterning and establish front–rear polarity through the coupling of biochemical signalling and mechanical feedback. We develop a cross-scale mean-field framework that integrates RD signalling with cytosolic and cortical hydrodynamics to capture the evolution of cell shapes and emergent cellular locomotion. Our model reproduces diverse experimentally observed shape and motility phenotypes with small parameter changes, indicating that these behaviours correspond to self-organized limit cycles. Phase-space analysis reveals that coupling to both cytosolic flow and spatially varying surface tension is essential to recover the full spectrum of motility modes, providing a theoretical foundation for understanding amoeboid migration.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above.

blejzara/CellRDA-LS

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: c76d59f3ca568c1f6dedfc59c38473e5571f0987, 30 October 2025
Languages: C (3), C/C++ (2)
Size: 14 files, 5 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, license file, CITATION.cff
Not found: environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
7 files
At the source: github.com/blejzara/CellRDA-LS

Zenodo 18824137

License: CC-BY-4.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Size: 10 files
Software Heritage: not checked
Found in: “Data Availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)
  • 30 September 2026: the link answers (HTTP 200)
1 file
At the source:

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 5 scripts, each with its path and the digest of its content;
  • no match between paragraphs and code yet;
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

All data used to produce the figures presented in this manuscript are available at: DOI:10.5281/zenodo.18824137. Code used in simulations that provided presented results is available at: https://github.com/blejzara/CellRDA-LS.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 11 MeSH terms, 2 funders, 118 references.

Cite

This paper

Ivšić, B., Vulić, D., Weber, I., Nowakowski, P., & Smith, A.-S. (2026). Diversity in emergent cell locomotion from the coupling cytosolic and cortical Marangoni flows with reaction-diffusion dynamics. PLoS computational biology, 22(4), e1014216. https://doi.org/10.1371/journal.pcbi.1014216

BibTeX

@article{ivsic2026diversity,
author = {Ivšić, Blaž and Vulić, Dorijan and Weber, Igor and Nowakowski, Piotr and Smith, Ana-Sunčana},
title = {{Diversity in emergent cell locomotion from the coupling cytosolic and cortical Marangoni flows with reaction-diffusion dynamics}},
journal = {PLoS computational biology},
year = {2026},
month = apr,
volume = {22},
number = {4},
pages = {e1014216},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014216},
url = {https://doi.org/10.1371/journal.pcbi.1014216},
pmid = {42044154},
pmcid = {PMC13148778}
}

RIS

TY - JOUR
AU - Ivšić, Blaž
AU - Vulić, Dorijan
AU - Weber, Igor
AU - Nowakowski, Piotr
AU - Smith, Ana-Sunčana
TI - Diversity in emergent cell locomotion from the coupling cytosolic and cortical Marangoni flows with reaction-diffusion dynamics
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/04/27
VL - 22
IS - 4
SP - e1014216
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014216
UR - https://doi.org/10.1371/journal.pcbi.1014216
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014216",
"type": "article-journal",
"title": "Diversity in emergent cell locomotion from the coupling cytosolic and cortical Marangoni flows with reaction-diffusion dynamics",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Ivšić",
"given": "Blaž"
},
{
"family": "Vulić",
"given": "Dorijan"
},
{
"family": "Weber",
"given": "Igor"
},
{
"family": "Nowakowski",
"given": "Piotr"
},
{
"family": "Smith",
"given": "Ana-Sunčana"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "4",
"page": "e1014216",
"DOI": "10.1371/journal.pcbi.1014216",
"PMID": "42044154",
"PMCID": "PMC13148778",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014216",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
27
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41556-026-01965-1 [code]
Long-range mutual activation establishes Rho and Rac polarity during cell migration.
Journal: Nature cell biology
In common: cellular / molecular, 5 references
[2] doi:10.1038/s41467-026-72845-3 [code]
The membrane-to-cortex distance regulates mDia1 activity to control cortical mechanics.
Journal: Nature communications
In common: cellular / molecular, 4 references
[3] doi:10.1126/sciadv.aec7705 [code]
Mitotic Cdc42 waves encode PI(3,4)P&lt;sub&gt;2&lt;/sub&gt; signaling and Golgi morphological state to control spindle scaling.
Journal: Science advances
In common: cellular / molecular, 4 references
[4] doi:10.1002/advs.77115
RNA Promotes Synapsin Condensates That Organize Synaptic Vesicles in Live Neurons and Enable Localized Translation in Reconstituted System.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: cellular / molecular, 2 references
[5] doi:10.1186/s12860-026-00584-w
A role for PaxB in regulating blebbing: experimental insights and theoretical perspectives from Dictyostelium discoideum.
Journal: BMC molecular and cell biology
In common: none (in silico), 1 reference
[6] doi:10.1126/sciadv.aeb3209 [code]
The actomyosin cortex controls t-tubule remodeling in skeletal muscle.
Journal: Science advances
In common: cellular / molecular, 2 references
[7] doi:10.1085/jgp.202513926
Characterization of an open-channel structure and lateral conduction pathway in the cation-selective pentameric ligand-gated ion channel, ELIC.
Journal: The Journal of general physiology
In common: none (in silico), computational modeling (no new data), cellular / molecular
[8] doi:10.1016/j.celrep.2026.117793 [code]
Clustered inputs engage dendritic nonlinearities and calcium signaling to support efficient place-field formation in CA1 pyramidal neurons.
Journal: Cell reports
In common: none (in silico), computational modeling (no new data), cellular / molecular
[9] doi:10.1038/s41586-026-10670-w [code]
Zero-shot design of drug-binding proteins via neural iterative selection-expansion.
Journal: Nature
In common: none (in silico), computational modeling (no new data), cellular / molecular
[10] doi:10.1038/s41598-026-55218-0
Amorphous metal-oxide semiconductor thin-film neuron, synapse, and neuromorphic system.
Journal: Scientific reports
In common: none (in silico), computational modeling (no new data), cellular / molecular

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.