grackle/src/clib/api/solve_chemistry.cpp at c762a16fdf8b486be666a8654e89f87c41eeaa65 · grackle-project/grackle · 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
/***********************************************************************
/
/ Solve the chemistry and cooling
/
/
/ Copyright (c) 2013, Enzo/Grackle Development Team.
/
/ Distributed under the terms of the Enzo Public Licence.
/
/ The full license is in the file LICENSE, distributed with this
/ software.
************************************************************************/

#include <cstdio>
#include "grackle.h"
#include "internal_units.h"
#include "solve_rate_cool_g-cpp.h"
#include "update_UVbackground_rates.hpp"
#include "utils.h"

extern "C" int local_solve_chemistry(chemistry_data *my_chemistry,
                                     chemistry_data_storage *my_rates,
                                     code_units *my_units,
                                     grackle_field_data *my_fields,
                                     double dt_value)
{

  /* Return if this doesn't concern us. */

  if (!my_chemistry->use_grackle)
    return GR_SUCCESS;

  /* Update UV background rates. */
  photo_rate_storage my_uvb_rates;

  my_uvb_rates.k24 = my_uvb_rates.k25 = my_uvb_rates.k26 =
    my_uvb_rates.k27 = my_uvb_rates.k28 = my_uvb_rates.k29 =
    my_uvb_rates.k30 = my_uvb_rates.k31 = my_uvb_rates.piHI =
    my_uvb_rates.piHeI = my_uvb_rates.piHeII = my_uvb_rates.crsHI =
    my_uvb_rates.crsHeI = my_uvb_rates.crsHeII =
    my_uvb_rates.comp_xray = my_uvb_rates.temp_xray = 0.;

  if (my_chemistry->UVbackground == 1) {
    if (grackle::impl::update_UVbackground_rates(
          my_chemistry, my_rates, &my_uvb_rates, my_units) != GR_SUCCESS) {
      std::fprintf(stderr, "Error in update_UVbackground_rates.\n");
      return GR_FAIL;
    }
  }
  else {
    my_uvb_rates.k24       = my_rates->k24;
    my_uvb_rates.k25       = my_rates->k25;
    my_uvb_rates.k26       = my_rates->k26;
    my_uvb_rates.k27       = my_rates->k27;
    my_uvb_rates.k28       = my_rates->k28;
    my_uvb_rates.k29       = my_rates->k29;
    my_uvb_rates.k30       = my_rates->k30;
    my_uvb_rates.k31       = my_rates->k31;
    my_uvb_rates.piHI      = my_rates->piHI;
    my_uvb_rates.piHeI     = my_rates->piHeI;
    my_uvb_rates.piHeII    = my_rates->piHeII;
    my_uvb_rates.crsHI     = my_rates->crsHI;
    my_uvb_rates.crsHeI    = my_rates->crsHeI;
    my_uvb_rates.crsHeII   = my_rates->crsHeII;
    my_uvb_rates.comp_xray = my_rates->comp_xray;
    my_uvb_rates.temp_xray = my_rates->temp_xray;
  }

  /* Check for a metal field. */

  int metal_field_present = TRUE;
  if (my_fields->metal_density == NULL)
    metal_field_present = FALSE;

  InternalGrUnits internalu = new_internalu_(my_units);

  /* Error checking for H2 shielding approximation */
  if (self_shielding_err_check(my_chemistry, my_fields,
                               "local_solve_chemistry") != GR_SUCCESS) {
    return GR_SUCCESS;
  }

  /* Call the routine to solve cooling equations. */

  int ierr = solve_rate_cool_g(
    metal_field_present, dt_value, internalu,
    my_chemistry, my_rates, my_fields, &my_uvb_rates
  );

  if (ierr != GR_SUCCESS) {
    std::fprintf(stderr, "Error in solve_rate_cool_g.\n");
  }

  return ierr;

}

extern "C" int solve_chemistry(code_units *my_units,
                               grackle_field_data *my_fields,
                               double dt_value)
{
  if (local_solve_chemistry(grackle_data, &grackle_rates,
                            my_units, my_fields, dt_value) != GR_SUCCESS) {
    std::fprintf(stderr, "Error in local_solve_chemistry.\n");
    return GR_FAIL;
  }
  return GR_SUCCESS;
}