Skip to content

Commit

Permalink
Fix quartic/cubic solution output and read from netcdf. (erf-model#1576)
Browse files Browse the repository at this point in the history 
Co-authored-by: Aaron Lattanzi <[email protected]>
  • Loading branch information
@AMLattanzi
AMLattanzi and Aaron Lattanzi authored Apr 11, 2024
1 parent 6acf04b commit 19de6a6
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 20 deletions.
31 changes: 15 additions & 16 deletions Source/Radiation/Modal_aero_wateruptake.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,13 +311,12 @@ class ModalAeroWateruptake {
r=rdry*(1.+p*third/(1.-slog*rdry/a));
}
else {
// AML TODO: FIX WHATEVER THIS IS TRYING TO DO...
amrex::GpuComplex<real> cx4[4];
makoh_quartic(cx4,p43,p42,p41,p40);
//find smallest real(r8) solution
r = 1000.*rdry;
auto nsol = 0;
for(auto n=0; n<4; ++n) {
int nsol = -1;
for(int n=0; n<4; ++n) {
auto xr=cx4[n].real();
auto xi=cx4[n].imag();
if(abs(xi) > abs(xr)*eps) continue;
Expand All @@ -327,7 +326,7 @@ class ModalAeroWateruptake {
r=xr;
nsol=n;
}
if(nsol != 0) {
if(nsol == -1) {
amrex::Print() << "kohlerc: no real solution found(quartic), nsol= " << nsol << "\n";
r=rdry;
}
Expand All @@ -346,7 +345,7 @@ class ModalAeroWateruptake {
makoh_cubic(cx3,p32,p31,p30);
//find smallest real(r8) solution
r=1000.*rdry;
auto nsol = 0;
int nsol = -1;
for(auto n=0; n<3; ++n) {
auto xr = cx3[n].real();
auto xi = cx3[n].imag();
Expand All @@ -357,8 +356,8 @@ class ModalAeroWateruptake {
r=xr;
nsol=n;
}
if(nsol == 0) {
amrex::Print() << "kohlerc: no real solution found (cubic)\n";
if(nsol == -1) {
amrex::Print() << "kohlerc: no real solution found (cubic), nsol= " << nsol << "\n";
r=rdry;
}
}
Expand All @@ -374,10 +373,10 @@ class ModalAeroWateruptake {
// solves x**3 + p2 x**2 + p1 x + p0 = 0
// where p0, p1, p2 are real
static
void makoh_cubic(amrex::GpuComplex<real> cx[],
const real& p2,
const real& p1,
const real& p0)
void makoh_cubic (amrex::GpuComplex<real> cx[],
const real& p2,
const real& p1,
const real& p0)
{
const real eps = 1.0e-20;
const real third=1./3.;
Expand Down Expand Up @@ -413,11 +412,11 @@ class ModalAeroWateruptake {
// solves x**4 + p3 x**3 + p2 x**2 + p1 x + p0 = 0
// where p0, p1, p2, p3 are real
static
void makoh_quartic(amrex::GpuComplex<real> cx[],
const real& p3,
const real& p2,
const real& p1,
const real& p0)
void makoh_quartic (amrex::GpuComplex<real> cx[],
const real& p3,
const real& p2,
const real& p1,
const real& p0)
{
const real third=1./3.;
auto q = -p2*p2/36.+(p3*p1-4*p0)/12.;
Expand Down
41 changes: 37 additions & 4 deletions Source/Radiation/Phys_prop.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -714,6 +714,39 @@ class PhysProp {
prop.read(asmpsw, "asmpsw" );
prop.read(absplw, "absplw" );

/*
// DEBUG NETCDF READ ORDER
{
auto lb = asmpsw.get_lbounds();
auto ub = asmpsw.get_ubounds();
amrex::Print() << "CHECK: " << ub.size() << "\n";
for (int ind(1); ind<=ub.size(); ++ind) {
amrex::Print() << "BNDS asmpsw: " << ind << ' ' << lb(ind) << ' ' << ub(ind) << ' '
<< nswbnd << ' ' << prefi << ' ' << prefr << ' ' << ncoef << "\n";
}
}
{
auto lb = abspsw.get_lbounds();
auto ub = abspsw.get_ubounds();
amrex::Print() << "CHECK: " << ub.size() << "\n";
for (int ind(1); ind<=ub.size(); ++ind) {
amrex::Print() << "BNDS abspsw: " << ind << ' ' << lb(ind) << ' ' << ub(ind) << ' '
<< nswbnd << ' ' << prefi << ' ' << prefr << ' ' << ncoef << "\n";
}
}
{
auto lb = extpsw.get_lbounds();
auto ub = extpsw.get_ubounds();
amrex::Print() << "CHECK: " << ub.size() << "\n";
for (int ind(1); ind<=ub.size(); ++ind) {
amrex::Print() << "BNDS extpsw: " << ind << ' ' << lb(ind) << ' ' << ub(ind) << ' '
<< nswbnd << ' ' << prefi << ' ' << prefr << ' ' << ncoef << "\n";
}
}
*/

phys_prop.extpsw = real4d("extpsw", nswbnd, prefi, prefr, ncoef);
phys_prop.abspsw = real4d("abspsw", nswbnd, prefi, prefr, ncoef);
phys_prop.asmpsw = real4d("asmpsw", nswbnd, prefi, prefr, ncoef);
Expand All @@ -722,15 +755,15 @@ class PhysProp {
parallel_for (SimpleBounds<4>(nswbnd, prefi, prefr, ncoef),
YAKL_LAMBDA (int i, int j, int k, int l)
{
phys_prop.extpsw(i,j,k,l) = extpsw(i,1,j,k,l);
phys_prop.abspsw(i,j,k,l) = abspsw(i,1,j,k,l);
phys_prop.asmpsw(i,j,k,l) = asmpsw(i,1,j,k,l);
phys_prop.extpsw(i,j,k,l) = extpsw(l,k,j,1,i); //extpsw(i,1,j,k,l);
phys_prop.abspsw(i,j,k,l) = abspsw(l,k,j,1,i); //abspsw(i,1,j,k,l);
phys_prop.asmpsw(i,j,k,l) = asmpsw(l,k,j,1,i); //asmpsw(i,1,j,k,l);
});

parallel_for (SimpleBounds<4>(nlwbnd, prefi, prefr, ncoef),
YAKL_LAMBDA (int i, int j, int k, int l)
{
phys_prop.absplw(i,j,k,l) = absplw(i,1,j,k,l);
phys_prop.absplw(i,j,k,l) = absplw(l,k,j,1,i); //absplw(i,1,j,k,l);
});

prop.read(phys_prop.refrtabsw, "refindex_real_sw" );
Expand Down

0 comments on commit 19de6a6

Please sign in to comment.