Skip to content

Commit

Permalink
liquid methane examples is working. got rid of that bad old full leap…
Browse files Browse the repository at this point in the history 
…r function
  • Loading branch information
@ameliajo
ameliajo committed Sep 7, 2020
1 parent b9f3121 commit 7b3ae58
Show file tree
Hide file tree
Showing 3 changed files with 103 additions and 983 deletions.
154 changes: 15 additions & 139 deletions src/leapr.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,139 +9,13 @@
#include <variant>
#include "endout/endout.h"

template<typename Range, typename Float, typename RangeZipped>
auto leaprTempLoop( int nphon, Float awr, int iel, int npr, int ncold, int lat, Range alpha, Range beta, Range temps, Float delta, Range rho, Float transWgt, Float diffusion_const, Float continWgt, RangeZipped oscEnergiesWgts, Float dka, Range kappaVals, Float cfrac, Float arat ){

Float tBar, effectiveTemp, temp, tev, sc, scaling;
std::vector<Float> sab( alpha.size()*beta.size(),0.0),
sab2(alpha.size()*beta.size(),0.0);
Range effectiveTemps(temps.size()), lambdaVals(temps.size());
std::variant<Range,bool> braggOutput;

std::vector<std::vector<Float>> sab_AllTemps(temps.size());


for ( size_t itemp = 0; itemp < temps.size(); ++itemp ){
temp = temps[itemp];
tev = kb * temp;

sc = (lat == 1) ? 0.0253/tev : 1.0;
scaling = sc/arat;
for ( auto& a : alpha ){ a *= scaling; }
for ( auto& b : beta ){ b *= sc; }

delta /= tev;

//---------------- Incoherent (Elastic and Inelastic) ----------------------
auto continOutput = contin(nphon, delta, continWgt, rho, alpha, beta, sab);

lambdaVals[itemp] = std::get<0>(continOutput);
tBar = std::get<1>(continOutput);
effectiveTemp = temp*tBar;

if (transWgt > 0){
trans( alpha, beta, transWgt, delta, diffusion_const, lambdaVals[itemp],
continWgt, effectiveTemp, temp, sab );
}

if (oscEnergiesWgts.size() > 0){
discre( lambdaVals[itemp], transWgt, continWgt, alpha, beta, temp,
oscEnergiesWgts, effectiveTemp, sab );
}

//--------------- Coherent (Inelastic Approximations) ----------------------
if (ncold > 0){
bool free = false;
coldh( tev, ncold, transWgt+continWgt, alpha, beta, dka, kappaVals, free,
sab, sab2, tBar );
}
else if (kappaVals.size() > 0){
skold( cfrac, tev, alpha, beta, kappaVals, awr, dka, sab );
}

//-------------------------- Coherent (Elastic) ----------------------------
if (iel > 0){
double emax = 5.0;
std::vector<double> bragg ( 60000, 0.0 );
auto coherOut = coher( iel, npr, bragg, emax );
int numVals = std::get<0>(coherOut);
bragg.resize(numVals);
braggOutput = bragg;
}
else { braggOutput = false; }

effectiveTemps[itemp] = effectiveTemp;

sab_AllTemps[itemp] = sab;

}
return std::make_tuple(sab,effectiveTemps,lambdaVals,braggOutput,sab_AllTemps);
}








template<typename Range, typename Float, typename RangeZipped>
auto leapr( int nphon, Float awr, int iel, int npr, int ncold, Float aws, int lat, Range alpha, Range beta, Range temps, Float delta, Range rho, Float transWgt, Float diffusion_const, Float continWgt, RangeZipped oscEnergiesWgts, Float dka, Range kappaVals, Float cfrac, std::tuple<int,int,Range> secondaryScatterInput ){

//std::variant<Range,std::tuple<Range,Range>> sab, effectiveTemps, lambdaVals;
Range sab, effectiveTemps, lambdaVals;
std::variant<Range,bool> braggOutput;

Range sab2, effectiveTemps2, lambdaVals2;
std::variant<Range,bool> braggOutput2;

//-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-
//-/-/-/-/-/-/-/-/-Principal scatterer (isecs = 0)-/-/-/-/-/-/-/-/-/-/-/-/-/-
//-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-
auto out1 = leaprTempLoop( nphon, awr, iel, npr, ncold, lat, alpha, beta,
temps, delta, rho, transWgt, diffusion_const, continWgt, oscEnergiesWgts,
dka, kappaVals, cfrac, 1.0 );

sab = std::get<0>(out1); effectiveTemps = std::get<1>(out1);
lambdaVals = std::get<2>(out1); braggOutput = std::get<3>(out1);
auto sab_AllTemps = std::get<4>(out1);

// if number of secondary scatterers is 0 or user wants free gas or diffusion
// for the secondary scatterer
int nss = std::get<0>(secondaryScatterInput),
b7 = std::get<1>(secondaryScatterInput);
if (nss == 0 or b7 > 0){
return std::make_tuple(sab, effectiveTemps, lambdaVals, braggOutput,
sab2,effectiveTemps2,lambdaVals2,braggOutput2,
sab_AllTemps);
}

//-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-
//-/-/-/-/-/-/-/-/-Secondary scatterer (isecs = 1)-/-/-/-/-/-/-/-/-/-/-/-/-/-
//-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-\-

Range secondaryRho = std::get<2>(secondaryScatterInput);
auto out2 = leaprTempLoop( nphon, awr, iel, npr, ncold, lat, alpha, beta,
temps, delta, secondaryRho, transWgt, diffusion_const, continWgt,
oscEnergiesWgts, dka, kappaVals, cfrac, aws/awr );

sab2 = std::get<0>(out2); effectiveTemps2 = std::get<1>(out2);
lambdaVals2 = std::get<2>(out2); braggOutput2 = std::get<3>(out2);

return std::make_tuple(sab, effectiveTemps, lambdaVals, braggOutput,
sab2,effectiveTemps2,lambdaVals2,braggOutput2,
sab_AllTemps);



}


template <typename Range, typename RangeOfRanges, typename Float >
auto full_LEAPR( std::vector<int> generalInfo, std::vector<int> scatterControl,
Range scatterInfo, Range temps, Range alphas, Range betas,
RangeOfRanges rhoVec, Range rho_dx_vec, RangeOfRanges transInfo,
RangeOfRanges oscE_vec, RangeOfRanges oscW_vec, Float smin, std::tuple<std::vector<double>,double> kappaInfo ){
template <typename Range, typename RangeOfRanges, typename Float>
auto leapr( std::vector<int> generalInfo, std::vector<int> scatterControl,
Range scatterInfo, Range temps, Range alphas, Range betas, RangeOfRanges rhoVec,
Range rho_dx_vec, RangeOfRanges transInfo, RangeOfRanges oscE_vec,
RangeOfRanges oscW_vec, Float smin, std::tuple<std::vector<double>,double>
kappaInfo = std::tuple<std::vector<double>,double>() ){

std::cout.precision(15);
int nphon = generalInfo[0];
Expand Down Expand Up @@ -212,6 +86,8 @@ auto full_LEAPR( std::vector<int> generalInfo, std::vector<int> scatterControl,
double dka = std::get<1>(kappaInfo);
coldh( tev, ncold, transWgt+continWgt, alphas, betas, dka, kappa, free,
sab, sab2, tempf[itemp]);
//std::cout << (sab|ranges::view::all) << std::endl;
//std::cout << (sab2|ranges::view::all) << std::endl;
}


Expand Down Expand Up @@ -254,11 +130,11 @@ auto full_LEAPR( std::vector<int> generalInfo, std::vector<int> scatterControl,
int numEdges = bragg.size();
//std::cout << numEdges << std::endl;

try {
std::get<bool>(braggOutput); // w contains int, not float: will throw
//std::cout << "No bragg!" << std::endl;
}
catch (const std::bad_variant_access&) {}
//try {
// std::get<bool>(braggOutput); // w contains int, not float: will throw
// //std::cout << "No bragg!" << std::endl;
//}
//catch (const std::bad_variant_access&) {}



Expand All @@ -285,15 +161,15 @@ auto full_LEAPR( std::vector<int> generalInfo, std::vector<int> scatterControl,
std::vector<unsigned int> numAtomsVec {npr,mss};
if (numSecondaryScatterers == 0){ numAtomsVec.resize(1); }

/*
njoy::ENDFtk::file::Type<7> MF7 = endout(sab_AllTemps,za,awrVec,spr,sps,temps,
numSecondaryScatterers, b7 , sab_AllTemps,alphas,
betas,dwpix,dwp1,iel,transWgt,bragg,nedge,tempf,tempf1,ilog,
isym,lat,numAtomsVec);

return MF7;
*/
/*
return;
*/
std::cout << numEdges << std::endl;
std::cout << ilog<< std::endl;
std::cout << za<< std::endl;
Expand Down

0 comments on commit 7b3ae58

Please sign in to comment.