partmc/src/camp_debug.h at ea78eecd9e1c8bed5c1aa3755d2e738e77de534e · compdyn/partmc · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
/* Copyright (C) 2015-2018 Matthew Dawson
 * Licensed under the GNU General Public License version 2 or (at your
 * option) any later version. See the file COPYING for details.
 *
 * Header file with some debugging functions for use with camp_solver.c
 *
 */
#ifndef CAMP_DEBUG_H
#define CAMP_DEBUG_H
#ifdef PMC_DEBUG
#define PMC_DEBUG_SPEC_ 0
#define PMC_DEBUG_PRINT(x) \
  pmc_debug_print(sd->cvode_mem, x, false, 0, __LINE__, __func__)
#define PMC_DEBUG_PRINT_INT(x, y) \
  pmc_debug_print(sd->cvode_mem, x, false, y, __LINE__, __func__)
#define PMC_DEBUG_PRINT_FULL(x) \
  pmc_debug_print(sd->cvode_mem, x, true, 0, __LINE__, __func__)
#define PMC_DEBUG_JAC_STRUCT(J, x) pmc_debug_print_jac_struct((void *)sd, J, x)
#define PMC_DEBUG_JAC(J, x) pmc_debug_print_jac((void *)sd, J, x)
void pmc_debug_print(void *cvode_mem, const char *message, bool do_full,
                     const int int_val, const int line, const char *func) {
#ifdef PMC_USE_SUNDIALS
  CVodeMem cv_mem = (CVodeMem)cvode_mem;
  if (!(cv_mem->cv_debug_out)) return;
  printf(
      "\n[DEBUG] line %4d in %-20s(): %-25s %-4.0d t_n = %le h = %le q = %d "
      "hin = %le species %d(zn[0] = %le zn[1] = %le tempv = %le tempv1 = %le "
      "tempv2 = %le acor_init = %le last_yn = %le",
      line, func, message, int_val, cv_mem->cv_tn, cv_mem->cv_h, cv_mem->cv_q,
      cv_mem->cv_hin, PMC_DEBUG_SPEC_,
      NV_DATA_S(cv_mem->cv_zn[0])[PMC_DEBUG_SPEC_],
      NV_DATA_S(cv_mem->cv_zn[1])[PMC_DEBUG_SPEC_],
      NV_DATA_S(cv_mem->cv_tempv)[PMC_DEBUG_SPEC_],
      NV_DATA_S(cv_mem->cv_tempv1)[PMC_DEBUG_SPEC_],
      // NV_DATA_S(cv_mem->cv_tempv2)[PMC_DEBUG_SPEC_],
      NV_DATA_S(cv_mem->cv_acor_init)[PMC_DEBUG_SPEC_],
      NV_DATA_S(cv_mem->cv_last_yn)[PMC_DEBUG_SPEC_]);
  if (do_full) {
    for (int i = 0; i < NV_LENGTH_S(cv_mem->cv_y); i++) {
      printf(
          "\n  zn[0][%3d] = % -le zn[1][%3d] = % -le tempv[%3d] = % -le "
          "tempv1[%3d] = % -le tempv2[%3d] = % -le acor_init[%3d] = % -le "
          "last_yn[%3d] = % -le",
          i, NV_DATA_S(cv_mem->cv_zn[0])[i], i, NV_DATA_S(cv_mem->cv_zn[1])[i],
          i, NV_DATA_S(cv_mem->cv_tempv)[i], i, NV_DATA_S(cv_mem->cv_tempv1)[i],
          // i, NV_DATA_S(cv_mem->cv_tempv2)[i], i,
          NV_DATA_S(cv_mem->cv_acor_init)[i], i,
          NV_DATA_S(cv_mem->cv_last_yn)[i]);
    }
  }
#endif
}
void pmc_debug_print_jac_struct(void *solver_data, SUNMatrix J,
                                const char *message) {
#ifdef PMC_USE_SUNDIALS
  SolverData *sd = (SolverData *)solver_data;

  if (!(sd->debug_out)) return;
  int n_state_var = SM_COLUMNS_S(J);
  int i_elem = 0;
  int next_col = 0;
  printf("\n\n   Jacobian structure (↓ind →dep) - %s\n     ", message);
  for (int i_dep = 0; i_dep < n_state_var; i_dep++) printf("[%3d]", i_dep);
  for (int i_ind = 0; i_ind < n_state_var; i_ind++) {
    printf("\n[%3d]", i_ind);
    next_col = SM_INDEXPTRS_S(J)[i_ind + 1];
    for (int i_dep = 0; i_dep < n_state_var; i_dep++) {
      if (i_dep == SM_INDEXVALS_S(J)[i_elem] && i_elem < next_col) {
        printf(" %3d ", i_elem++);
      } else {
        printf("  -  ");
      }
    }
  }
#endif
}
void pmc_debug_print_jac(void *solver_data, SUNMatrix J, const char *message) {
#ifdef PMC_USE_SUNDIALS
  SolverData *sd = (SolverData *)solver_data;

  if (!(sd->debug_out)) return;
  int n_state_var = SM_COLUMNS_S(J);
  int i_elem = 0;
  int next_col = 0;
  printf("\n\n   Jacobian (↓ind →dep) - %s\n     ", message);
  for (int i_dep = 0; i_dep < n_state_var; i_dep++)
    printf("      [%3d]", i_dep);
  for (int i_ind = 0; i_ind < n_state_var; i_ind++) {
    printf("\n[%3d]   ", i_ind);
    next_col = SM_INDEXPTRS_S(J)[i_ind + 1];
    for (int i_dep = 0; i_dep < n_state_var; i_dep++) {
      if (i_dep == SM_INDEXVALS_S(J)[i_elem] && i_elem < next_col) {
        printf(" % -1.2le ", SM_DATA_S(J)[i_elem++]);
      } else {
        printf("     -     ");
      }
    }
  }
#endif
}

#else
#define PMC_DEBUG_PRINT(x)
#define PMC_DEBUG_PRINT_INT(x, y)
#define PMC_DEBUG_PRINT_FULL(x)
#define PMC_DEBUG_JAC_STRUCT(J, x)
#define PMC_DEBUG_JAC(J, x)
#endif

/** \brief Print some camp-chem data sizes
 *
 * \param md Pointer to the model data
 */
static void print_data_sizes(ModelData *md) {
  int *ptr = md->rxn_int_data;
  int n_rxn = ptr[0];

  printf("n_rxn: %d ", n_rxn);
  printf("n_state_var: %d", md->n_per_cell_state_var * md->n_cells);
  printf("n_dep_var: %d", md->n_per_cell_dep_var * md->n_cells);
}

/** \brief Print Jacobian matrix in format KLU SPARSE
 *
 * \param M Jacobian matrix
 */
static void print_jacobian(SUNMatrix M) {
  printf("\n NNZ JAC: %lld \n", SM_NNZ_S(M));
  printf("DATA | INDEXVALS:\n");
  for (int i = 0; i < SM_NNZ_S(M); i++) {
    printf("% -le \n", (SM_DATA_S(M))[i]);
    printf("%lld \n", (SM_INDEXVALS_S(M))[i]);
  }
  printf("PTRS:\n");
  for (int i = 0; i <= SM_NP_S(M); i++) {
    printf("%lld \n", (SM_INDEXPTRS_S(M))[i]);
  }
}

/** \brief Print derivative array
 *
 * \param deriv Derivative array
 */
static void print_derivative(N_Vector deriv) {
  // printf(" deriv length: %d\n", NV_LENGTH_S(deriv));
  for (int i = 0; i < NV_LENGTH_S(deriv); i++) {  // NV_LENGTH_S(deriv)
    printf(" deriv: % -le", NV_DATA_S(deriv)[i]);
    printf(" index: %d \n", i);
  }
}

#endif