Skip to content

Commit 6c71986

Browse files
tradowsktradowsk
committed
Add unit tests for AbsorptionCoeff parameter derivatives
1 parent fe5aac7 commit 6c71986

File tree

1 file changed

+159
-2
lines changed

1 file changed

+159
-2
lines changed

Diff for: ‎test/unit/spectroscopic_absorption_test.C

+159-2
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@
2626

2727
#ifdef GRINS_HAVE_CPPUNIT
2828

29+
#ifdef GRINS_HAVE_ANTIOCH
30+
2931
#include <libmesh/ignore_warnings.h>
3032
#include <cppunit/extensions/HelperMacros.h>
3133
#include <cppunit/TestCase.h>
@@ -34,12 +36,18 @@
3436
#include "test_comm.h"
3537
#include "grins_test_paths.h"
3638

39+
#include "absorption_coeff_testing.h"
40+
3741
// GRINS
3842
#include "grins/math_constants.h"
3943
#include "grins/grins_enums.h"
4044
#include "grins/simulation_builder.h"
4145
#include "grins/simulation.h"
4246
#include "grins/variable_warehouse.h"
47+
#include "grins/single_variable.h"
48+
#include "grins/multicomponent_variable.h"
49+
#include "grins/chemistry_builder.h"
50+
#include "grins/antioch_chemistry.h"
4351

4452
// libMesh
4553
#include "libmesh/parsed_function.h"
@@ -53,10 +61,9 @@ namespace GRINSTesting
5361
{
5462
public:
5563
CPPUNIT_TEST_SUITE( SpectroscopicAbsorptionTest );
56-
#if GRINS_HAVE_ANTIOCH
5764
CPPUNIT_TEST( single_elem_mesh );
5865
CPPUNIT_TEST( multi_elem_mesh );
59-
#endif
66+
CPPUNIT_TEST( param_derivs );
6067
CPPUNIT_TEST_SUITE_END();
6168

6269
public:
@@ -85,6 +92,43 @@ namespace GRINSTesting
8592
this->run_test(filename,calc_answer);
8693
}
8794

95+
void param_derivs()
96+
{
97+
const std::string filename = std::string(GRINS_TEST_UNIT_INPUT_SRCDIR)+"/spectroscopic_absorption_qoi.in";
98+
this->init_sim(filename);
99+
100+
std::string material = "TestMaterial";
101+
std::string hitran_data = "./test_data/CO2_data.dat";
102+
std::string hitran_partition = "./test_data/CO2_partition_function.dat";
103+
libMesh::Real T_min = 290.0,
104+
T_max = 310.0,
105+
T_step = 0.01;
106+
GRINS::SharedPtr<GRINS::HITRAN> hitran( new GRINS::HITRAN(hitran_data,hitran_partition,T_min,T_max,T_step) );
107+
108+
std::string species = "CO2";
109+
libMesh::Real thermo_pressure = 5066.25;
110+
libMesh::Real nu_desired = 3682.7649;
111+
libMesh::Real nu_min = 3682.69;
112+
libMesh::Real nu_max = 3682.8;
113+
GRINS::ChemistryBuilder chem_builder;
114+
libMesh::UniquePtr<GRINS::AntiochChemistry> chem_ptr;
115+
chem_builder.build_chemistry(*(_input.get()),material,chem_ptr);
116+
GRINS::SharedPtr<GRINS::AntiochChemistry> chem(chem_ptr.release());
117+
GRINS::SharedPtr<AbsorptionCoeffTesting<GRINS::AntiochChemistry> > absorb = new AbsorptionCoeffTesting<GRINS::AntiochChemistry>(chem,hitran,nu_min,nu_max,nu_desired,species,thermo_pressure);
118+
119+
libMesh::Real T = 300.0;
120+
libMesh::Real P = 5066.25;
121+
122+
std::vector<libMesh::Real> Y = {0.0763662233, 0.9236337767};
123+
124+
for (unsigned int i=0; i<33; ++i)
125+
{
126+
this->T_param_derivatives(absorb,T,P,Y,i);
127+
this->P_param_derivatives(absorb,T,P,Y,i);
128+
this->Y_param_derivatives(absorb,T,P,Y,i);
129+
}
130+
}
131+
88132
private:
89133
GRINS::SharedPtr<GRINS::Simulation> _sim;
90134
GRINS::SharedPtr<GetPot> _input;
@@ -114,10 +158,123 @@ namespace GRINSTesting
114158
*TestCommWorld );
115159
}
116160

161+
void T_param_derivatives(GRINS::SharedPtr<AbsorptionCoeffTesting<GRINS::AntiochChemistry> >absorb, libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> & Y, unsigned int i)
162+
{
163+
libMesh::Real delta = 1.0e-6;
164+
165+
std::vector<libMesh::Real> T_analytic;
166+
T_analytic.push_back(absorb->dS_dT(T,P,i));
167+
T_analytic.push_back(absorb->d_nuC_dT(T,P,Y,i));
168+
T_analytic.push_back(absorb->d_nuD_dT(T,P,i));
169+
T_analytic.push_back(absorb->d_voigt_a_dT(T,P,Y,i));
170+
T_analytic.push_back(absorb->d_voigt_w_dT(T,P,i));
171+
T_analytic.push_back(absorb->d_voigt_dT(T,P,Y,i));
172+
T_analytic.push_back(absorb->d_kv_dT(T,P,Y,i));
173+
174+
std::vector<libMesh::Real> fd_plus;
175+
fd_plus.push_back(absorb->Sw(T+delta,P,i));
176+
fd_plus.push_back(absorb->nu_C(T+delta,P,Y,i));
177+
fd_plus.push_back(absorb->nu_D(T+delta,P,i));
178+
fd_plus.push_back(absorb->voigt_a(T+delta,P,Y,i));
179+
fd_plus.push_back(absorb->voigt_w(T+delta,P,i));
180+
fd_plus.push_back(absorb->voigt(T+delta,P,Y,i));
181+
fd_plus.push_back(absorb->kv(T+delta,P,Y,i));
182+
183+
std::vector<libMesh::Real> fd_minus;
184+
fd_minus.push_back(absorb->Sw(T-delta,P,i));
185+
fd_minus.push_back(absorb->nu_C(T-delta,P,Y,i));
186+
fd_minus.push_back(absorb->nu_D(T-delta,P,i));
187+
fd_minus.push_back(absorb->voigt_a(T-delta,P,Y,i));
188+
fd_minus.push_back(absorb->voigt_w(T-delta,P,i));
189+
fd_minus.push_back(absorb->voigt(T-delta,P,Y,i));
190+
fd_minus.push_back(absorb->kv(T-delta,P,Y,i));
191+
192+
check_param_derivatives(T_analytic,fd_plus,fd_minus,delta);
193+
}
194+
195+
void P_param_derivatives(GRINS::SharedPtr<AbsorptionCoeffTesting<GRINS::AntiochChemistry> >absorb, libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> & Y, unsigned int i)
196+
{
197+
libMesh::Real delta = 1.0e-3;
198+
199+
std::vector<libMesh::Real> P_analytic;
200+
P_analytic.push_back(absorb->dS_dP(T,P,i));
201+
P_analytic.push_back(absorb->d_nuC_dP(T,Y,i));
202+
P_analytic.push_back(absorb->d_nuD_dP(T,i));
203+
P_analytic.push_back(absorb->d_voigt_a_dP(T,P,Y,i));
204+
P_analytic.push_back(absorb->d_voigt_w_dP(T,P,i));
205+
P_analytic.push_back(absorb->d_voigt_dP(T,P,Y,i));
206+
P_analytic.push_back(absorb->d_kv_dP(T,P,Y,i));
207+
208+
std::vector<libMesh::Real> fd_plus;
209+
fd_plus.push_back(absorb->Sw(T,P+delta,i));
210+
fd_plus.push_back(absorb->nu_C(T,P+delta,Y,i));
211+
fd_plus.push_back(absorb->nu_D(T,P+delta,i));
212+
fd_plus.push_back(absorb->voigt_a(T,P+delta,Y,i));
213+
fd_plus.push_back(absorb->voigt_w(T,P+delta,i));
214+
fd_plus.push_back(absorb->voigt(T,P+delta,Y,i));
215+
fd_plus.push_back(absorb->kv(T,P+delta,Y,i));
216+
217+
std::vector<libMesh::Real> fd_minus;
218+
fd_minus.push_back(absorb->Sw(T,P-delta,i));
219+
fd_minus.push_back(absorb->nu_C(T,P-delta,Y,i));
220+
fd_minus.push_back(absorb->nu_D(T,P-delta,i));
221+
fd_minus.push_back(absorb->voigt_a(T,P-delta,Y,i));
222+
fd_minus.push_back(absorb->voigt_w(T,P-delta,i));
223+
fd_minus.push_back(absorb->voigt(T,P-delta,Y,i));
224+
fd_minus.push_back(absorb->kv(T,P-delta,Y,i));
225+
226+
check_param_derivatives(P_analytic,fd_plus,fd_minus,delta);
227+
}
228+
229+
void Y_param_derivatives(GRINS::SharedPtr<AbsorptionCoeffTesting<GRINS::AntiochChemistry> >absorb, libMesh::Real T, libMesh::Real P, std::vector<libMesh::Real> & Y, unsigned int i)
230+
{
231+
libMesh::Real delta = 1.0e-8;
232+
233+
libMesh::Real y_base = Y[0];
234+
std::vector<libMesh::Real> Y_analytic;
235+
Y_analytic.push_back(absorb->dX_dY(Y));
236+
Y_analytic.push_back(absorb->d_nuC_dY(T,P,Y,i));
237+
Y_analytic.push_back(absorb->d_voigt_a_dY(T,P,Y,i));
238+
Y_analytic.push_back(absorb->d_voigt_dY(T,P,Y,i));
239+
Y_analytic.push_back(absorb->d_kv_dY(T,P,Y,i));
240+
241+
std::vector<libMesh::Real> fd_plus;
242+
Y[0] = y_base + delta;
243+
fd_plus.push_back(absorb->_chemistry->X(absorb->_species_idx,absorb->_chemistry->M_mix(Y),Y[absorb->_species_idx]));
244+
fd_plus.push_back(absorb->nu_C(T,P,Y,i));
245+
fd_plus.push_back(absorb->voigt_a(T,P,Y,i));
246+
fd_plus.push_back(absorb->voigt(T,P,Y,i));
247+
fd_plus.push_back(absorb->kv(T,P,Y,i));
248+
249+
std::vector<libMesh::Real> fd_minus;
250+
Y[0] = y_base - delta;
251+
fd_minus.push_back(absorb->_chemistry->X(absorb->_species_idx,absorb->_chemistry->M_mix(Y),Y[absorb->_species_idx]));
252+
fd_minus.push_back(absorb->nu_C(T,P,Y,i));
253+
fd_minus.push_back(absorb->voigt_a(T,P,Y,i));
254+
fd_minus.push_back(absorb->voigt(T,P,Y,i));
255+
fd_minus.push_back(absorb->kv(T,P,Y,i));
256+
257+
Y[0] = y_base;
258+
259+
check_param_derivatives(Y_analytic,fd_plus,fd_minus,delta);
260+
}
261+
262+
void check_param_derivatives(std::vector<libMesh::Real> & analytic, std::vector<libMesh::Real> & fd_plus, std::vector<libMesh::Real> & fd_minus, libMesh::Real delta)
263+
{
264+
for (unsigned int d=0; d<analytic.size(); ++d)
265+
{
266+
libMesh::Real f_analytic = analytic[d];
267+
libMesh::Real fd_approx = (fd_plus[d] - fd_minus[d])/(2.0*delta);
268+
269+
CPPUNIT_ASSERT_DOUBLES_EQUAL(fd_approx,f_analytic,libMesh::TOLERANCE);
270+
}
271+
}
272+
117273
};
118274

119275
CPPUNIT_TEST_SUITE_REGISTRATION( SpectroscopicAbsorptionTest );
120276

121277
} // end namespace GRINSTesting
122278

279+
#endif // GRINS_HAVE_ANTIOCH
123280
#endif // GRINS_HAVE_CPPUNIT

0 commit comments

Comments
 (0)