diff --git a/exec/multispec/AdvanceTimestepBousq.cpp b/exec/multispec/AdvanceTimestepBousq.cpp index 4609fe14..0490df92 100644 --- a/exec/multispec/AdvanceTimestepBousq.cpp +++ b/exec/multispec/AdvanceTimestepBousq.cpp @@ -9,7 +9,6 @@ #include #include -#include // argv contains the name of the inputs file entered at the command line diff --git a/src_multispec/ComputeDivReversibleStress.cpp~ b/src_multispec/ComputeDivReversibleStress.cpp~ deleted file mode 100644 index b70aa1f9..00000000 --- a/src_multispec/ComputeDivReversibleStress.cpp~ +++ /dev/null @@ -1,493 +0,0 @@ -#include "multispec_functions.H" - -void ComputeDivFHReversibleStress(std::array& div_reversible_stress, - const MultiFab& rhotot_in, -// Array2D& fh_kappa, - MultiFab& rho_in, - const Geometry& geom) -{ - BL_PROFILE_VAR("ComputeDivFHReversibleStress()",ComputeDivReversibleStress); - - BoxArray ba = rho_in.boxArray(); - DistributionMapping dmap = rho_in.DistributionMap(); - const GpuArray dx = geom.CellSizeArray(); - - MultiFab node_grad_x_mf(convert(ba,nodal_flag), dmap, nspecies, 1); - MultiFab node_grad_y_mf(convert(ba,nodal_flag), dmap, nspecies, 1); - MultiFab node_grad_z_mf(convert(ba,nodal_flag), dmap, nspecies, 1); - - MultiFab conc(ba, dmap, nspecies, 2); - - // rho to conc - VALID REGION ONLY - ConvertRhoCToC(rho_in,rhotot_in,conc,1); - - Real scale_factor = rhobar[0]*k_B*T_init[0]/monomer_mass; - - // fill conc ghost cells - conc.FillBoundary(geom.periodicity()); - MultiFabPhysBCFH(conc,geom,0,nspecies,scale_factor); - -// JBB check what's in this -// Real scale_factor = rhobar[0]*k_B*T_init[0]/molmass[0]; - - for ( MFIter mfi(node_grad_x_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - - const Box& bx = mfi.growntilebox(1); - - const Array4& node_grad_x = node_grad_x_mf.array(mfi); - const Array4& node_grad_y = node_grad_y_mf.array(mfi); -#if (AMREX_SPACEDIM == 3) - const Array4& node_grad_z = node_grad_z_mf.array(mfi); -#endif - const Array4& c = conc.array(mfi); - - amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { -#if (AMREX_SPACEDIM == 2) - - node_grad_x(i,j,k,n) = (c(i,j,k,n)-c(i-1,j,k,n)+c(i,j-1,k,n)-c(i-1,j-1,k,n))/(2*dx[0]); - node_grad_y(i,j,k,n) = (c(i,j,k,n)-c(i,j-1,k,n)+c(i-1,j,k,n)-c(i-1,j-1,k,n))/(2*dx[1]); - -#elif (AMREX_SPACEDIM == 3) - - node_grad_x(i,j,k,n) = (c(i,j,k ,n)-c(i-1,j,k ,n)+c(i,j-1,k ,n)-c(i-1,j-1,k ,n) - +c(i,j,k-1,n)-c(i-1,j,k-1,n)+c(i,j-1,k-1,n)-c(i-1,j-1,k-1,n))/(4*dx[0]); - node_grad_y(i,j,k,n) = (c(i,j,k ,n)-c(i,j-1,k ,n)+c(i-1,j,k ,n)-c(i-1,j-1,k ,n) - +c(i,j,k-1,n)-c(i,j-1,k-1,n)+c(i-1,j,k-1,n)-c(i-1,j-1,k-1,n))/(4*dx[1]); - node_grad_z(i,j,k,n) = (c(i,j ,k,n)-c(i,j ,k-1,n)+c(i-1,j ,k,n)-c(i-1,j ,k-1,n) - +c(i,j-1,k,n)-c(i,j-1,k-1,n)+c(i-1,j-1,k,n)-c(i-1,j-1,k-1,n))/(4*dx[2]); -#endif - }); - } - -// JBB is this needed? - for (int n = 0 ; n < AMREX_SPACEDIM; n++) - div_reversible_stress[n].setVal(0.); - - for ( MFIter mfi(rhotot_in,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - - AMREX_D_TERM(const Array4& node_grad_x = node_grad_x_mf.array(mfi);, - const Array4& node_grad_y = node_grad_y_mf.array(mfi);, - const Array4& node_grad_z = node_grad_z_mf.array(mfi);); - - AMREX_D_TERM(const Array4 & forcex = div_reversible_stress[0].array(mfi);, - const Array4 & forcey = div_reversible_stress[1].array(mfi);, - const Array4 & forcez = div_reversible_stress[2].array(mfi);); - - AMREX_D_TERM(const Box & bx_x = mfi.nodaltilebox(0);, - const Box & bx_y = mfi.nodaltilebox(1);, - const Box & bx_z = mfi.nodaltilebox(2);); - -#if (AMREX_SPACEDIM == 2) - - amrex::ParallelFor(bx_x, bx_y, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for( int n=0 ; n< nspecies; n++){ - for( int m=0 ; m< nspecies; m++){ - Real cax_local_plus = 0.25*(node_grad_x(i,j,k,n)+node_grad_x(i,j+1,k,n)+node_grad_x(i+1,j+1,k,n)+node_grad_x(i+1,j,k,n)); - Real cay_local_plus = 0.25*(node_grad_y(i,j,k,n)+node_grad_y(i,j+1,k,n)+node_grad_y(i+1,j+1,k,n)+node_grad_y(i+1,j,k,n)); - Real cax_local_minus = 0.25*(node_grad_x(i,j,k,n)+node_grad_x(i,j+1,k,n)+node_grad_x(i-1,j+1,k,n)+node_grad_x(i-1,j,k,n)); - Real cay_local_minus = 0.25*(node_grad_y(i,j,k,n)+node_grad_y(i,j+1,k,n)+node_grad_y(i-1,j+1,k,n)+node_grad_y(i-1,j,k,n)); - Real cbx_local_plus = 0.25*(node_grad_x(i,j,k,m)+node_grad_x(i,j+1,k,m)+node_grad_x(i+1,j+1,k,m)+node_grad_x(i+1,j,k,m)); - Real cby_local_plus = 0.25*(node_grad_y(i,j,k,m)+node_grad_y(i,j+1,k,m)+node_grad_y(i+1,j+1,k,m)+node_grad_y(i+1,j,k,m)); - Real cbx_local_minus = 0.25*(node_grad_x(i,j,k,m)+node_grad_x(i,j+1,k,m)+node_grad_x(i-1,j+1,k,m)+node_grad_x(i-1,j,k,m)); - Real cby_local_minus = 0.25*(node_grad_y(i,j,k,m)+node_grad_y(i,j+1,k,m)+node_grad_y(i-1,j+1,k,m)+node_grad_y(i-1,j,k,m)); - - forcex(i,j,k) += scale_factor * fh_kappa(n,m)*( (node_grad_x(i,j+1,k,n)*node_grad_y(i,j+1,k,m) - node_grad_x(i,j,k,n)*node_grad_y(i,j,k,m))/dx[1] - +0.5*(cax_local_plus*cbx_local_plus-cay_local_plus*cby_local_plus -cax_local_minus*cbx_local_minus+cay_local_minus*cby_local_minus)/dx[0]); - } - } - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for( int n=0 ; n< nspecies; n++){ - for( int m=0 ; m< nspecies; m++){ - Real cax_local_plus = 0.25*(node_grad_x(i,j,k,n)+node_grad_x(i,j+1,k,n)+node_grad_x(i+1,j+1,k,n)+node_grad_x(i+1,j,k,n)); - Real cay_local_plus = 0.25*(node_grad_y(i,j,k,n)+node_grad_y(i,j+1,k,n)+node_grad_y(i+1,j+1,k,n)+node_grad_y(i+1,j,k,n)); - Real cax_local_minus = 0.25*(node_grad_x(i,j,k,n)+node_grad_x(i,j-1,k,n)+node_grad_x(i+1,j-1,k,n)+node_grad_x(i+1,j,k,n)); - Real cay_local_minus = 0.25*(node_grad_y(i,j,k,n)+node_grad_y(i,j-1,k,n)+node_grad_y(i+1,j-1,k,n)+node_grad_y(i+1,j,k,n)); - Real cbx_local_plus = 0.25*(node_grad_x(i,j,k,m)+node_grad_x(i,j+1,k,m)+node_grad_x(i+1,j+1,k,m)+node_grad_x(i+1,j,k,m)); - Real cby_local_plus = 0.25*(node_grad_y(i,j,k,m)+node_grad_y(i,j+1,k,m)+node_grad_y(i+1,j+1,k,m)+node_grad_y(i+1,j,k,m)); - Real cbx_local_minus = 0.25*(node_grad_x(i,j,k,m)+node_grad_x(i,j-1,k,m)+node_grad_x(i+1,j-1,k,m)+node_grad_x(i+1,j,k,m)); - Real cby_local_minus = 0.25*(node_grad_y(i,j,k,m)+node_grad_y(i,j-1,k,m)+node_grad_y(i+1,j-1,k,m)+node_grad_y(i+1,j,k,m)); - - forcey(i,j,k) += scale_factor*fh_kappa(n,m)*((node_grad_x(i+1,j,k,n)*node_grad_y(i+1,j,k,m) - node_grad_x(i,j,k,n)*node_grad_y(i,j,k,m))/dx[0] - +0.5*(cay_local_plus*cby_local_plus- cay_local_minus*cby_local_minus-cax_local_plus*cbx_local_plus+cax_local_minus*cbx_local_minus)/dx[1]); - } - } - }); - -#elif (AMREX_SPACEDIM == 3) - - amrex::ParallelFor(bx_x, bx_y, bx_z, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for( int n=0 ; n< nspecies; n++){ - for( int m=0 ; m< nspecies; m++){ - - Real cax_local_plus = .125*(node_grad_x(i ,j ,k ,n)+node_grad_x(i ,j+1,k ,n)+node_grad_x(i+1,j+1,k ,n)+node_grad_x(i+1,j ,k ,n) - +node_grad_x(i ,j ,k+1,n)+node_grad_x(i ,j+1,k+1,n)+node_grad_x(i+1,j+1,k+1,n)+node_grad_x(i+1,j ,k+1,n)); - - Real cax_local_minus = .125*(node_grad_x(i ,j ,k ,n)+node_grad_x(i ,j+1,k ,n)+node_grad_x(i-1,j+1,k ,n)+node_grad_x(i-1,j ,k ,n) - +node_grad_x(i ,j ,k+1,n)+node_grad_x(i ,j+1,k+1,n)+node_grad_x(i-1,j+1,k+1,n)+node_grad_x(i-1,j ,k+1,n)); - - Real cay_local_plus = .125*(node_grad_y(i ,j ,k ,n)+node_grad_y(i ,j+1,k ,n)+node_grad_y(i+1,j+1,k ,n)+node_grad_y(i+1,j ,k ,n) - +node_grad_y(i ,j ,k+1,n)+node_grad_y(i ,j+1,k+1,n)+node_grad_y(i+1,j+1,k+1,n)+node_grad_y(i+1,j ,k+1,n)); - - Real cay_local_minus = .125*(node_grad_y(i ,j ,k ,n)+node_grad_y(i ,j+1,k ,n)+node_grad_y(i-1,j+1,k ,n)+node_grad_y(i-1,j ,k ,n) - +node_grad_y(i ,j ,k+1,n)+node_grad_y(i ,j+1,k+1,n)+node_grad_y(i-1,j+1,k+1,n)+node_grad_y(i-1,j ,k+1,n)); - - Real caz_local_plus = .125*(node_grad_z(i ,j ,k ,n)+node_grad_z(i ,j+1,k ,n)+node_grad_z(i+1,j+1,k ,n)+node_grad_z(i+1,j ,k ,n) - +node_grad_z(i ,j ,k+1,n)+node_grad_z(i ,j+1,k+1,n)+node_grad_z(i+1,j+1,k+1,n)+node_grad_z(i+1,j ,k+1,n)); - - Real caz_local_minus = .125*(node_grad_z(i ,j ,k ,n)+node_grad_z(i ,j+1,k ,n)+node_grad_z(i-1,j+1,k ,n)+node_grad_z(i-1,j ,k ,n) - +node_grad_z(i ,j ,k+1,n)+node_grad_z(i ,j+1,k+1,n)+node_grad_z(i-1,j+1,k+1,n)+node_grad_z(i-1,j ,k+1,n)); - - Real cbx_local_plus = .125*(node_grad_x(i ,j ,k ,m)+node_grad_x(i ,j+1,k ,m)+node_grad_x(i+1,j+1,k ,m)+node_grad_x(i+1,j ,k ,m) - +node_grad_x(i ,j ,k+1,m)+node_grad_x(i ,j+1,k+1,m)+node_grad_x(i+1,j+1,k+1,m)+node_grad_x(i+1,j ,k+1,m)); - - Real cbx_local_minus = .125*(node_grad_x(i ,j ,k ,m)+node_grad_x(i ,j+1,k ,m)+node_grad_x(i-1,j+1,k ,m)+node_grad_x(i-1,j ,k ,m) - +node_grad_x(i ,j ,k+1,m)+node_grad_x(i ,j+1,k+1,m)+node_grad_x(i-1,j+1,k+1,m)+node_grad_x(i-1,j ,k+1,m)); - - Real cby_local_plus = .125*(node_grad_y(i ,j ,k ,m)+node_grad_y(i ,j+1,k ,m)+node_grad_y(i+1,j+1,k ,m)+node_grad_y(i+1,j ,k ,m) - +node_grad_y(i ,j ,k+1,m)+node_grad_y(i ,j+1,k+1,m)+node_grad_y(i+1,j+1,k+1,m)+node_grad_y(i+1,j ,k+1,m)); - - Real cby_local_minus = .125*(node_grad_y(i ,j ,k ,m)+node_grad_y(i ,j+1,k ,m)+node_grad_y(i-1,j+1,k ,m)+node_grad_y(i-1,j ,k ,m) - +node_grad_y(i ,j ,k+1,m)+node_grad_y(i ,j+1,k+1,m)+node_grad_y(i-1,j+1,k+1,m)+node_grad_y(i-1,j ,k+1,m)); - - Real cbz_local_plus = .125*(node_grad_z(i ,j ,k ,m)+node_grad_z(i ,j+1,k ,m)+node_grad_z(i+1,j+1,k ,m)+node_grad_z(i+1,j ,k ,m) - +node_grad_z(i ,j ,k+1,m)+node_grad_z(i ,j+1,k+1,m)+node_grad_z(i+1,j+1,k+1,m)+node_grad_z(i+1,j ,k+1,m)); - - Real cbz_local_minus = .125*(node_grad_z(i ,j ,k ,m)+node_grad_z(i ,j+1,k ,m)+node_grad_z(i-1,j+1,k ,m)+node_grad_z(i-1,j ,k ,m) - +node_grad_z(i ,j ,k+1,m)+node_grad_z(i ,j+1,k+1,m)+node_grad_z(i-1,j+1,k+1,m)+node_grad_z(i-1,j ,k+1,m)); - - forcex(i,j,k) += scale_factor * fh_kappa(n,m)*0.5* - ( ((node_grad_x(i,j+1,k,n)*node_grad_y(i,j+1,k,m) + node_grad_x(i,j+1,k+1,n)*node_grad_y(i,j+1,k+1,m)) - - (node_grad_x(i,j ,k,n)*node_grad_y(i,j ,k,m) + node_grad_x(i,j ,k+1,n)*node_grad_y(i,j ,k+1,m)))/dx[1] - +((cax_local_plus *cbx_local_plus - caz_local_plus *cbz_local_plus -cay_local_plus *cby_local_plus) - -(cax_local_minus*cbx_local_minus- caz_local_minus*cbz_local_minus-cay_local_minus*cby_local_minus))/(dx[0]) - + ((node_grad_x(i,j,k+1,n)*node_grad_z(i,j,k+1,m) + node_grad_x(i,j+1,k+1,n)*node_grad_z(i,j+1,k+1,m)) - - (node_grad_x(i,j,k ,n)*node_grad_z(i,j,k ,m) + node_grad_x(i,j+1,k ,n)*node_grad_z(i,j+1,k ,m)))/dx[2]); - - } - } - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for( int n=0 ; n< nspecies; n++){ - for( int m=0 ; m< nspecies; m++){ - - Real cax_local_plus = .125*(node_grad_x(i ,j ,k ,n)+node_grad_x(i ,j+1,k ,n)+node_grad_x(i+1,j+1,k ,n)+node_grad_x(i+1,j ,k ,n) - +node_grad_x(i ,j ,k+1,n)+node_grad_x(i ,j+1,k+1,n)+node_grad_x(i+1,j+1,k+1,n)+node_grad_x(i+1,j ,k+1,n)); - - Real cax_local_minus = .125*(node_grad_x(i ,j ,k ,n)+node_grad_x(i ,j-1,k ,n)+node_grad_x(i+1,j-1,k ,n)+node_grad_x(i+1,j ,k ,n) - +node_grad_x(i ,j ,k+1,n)+node_grad_x(i ,j-1,k+1,n)+node_grad_x(i+1,j-1,k+1,n)+node_grad_x(i+1,j ,k+1,n)); - - Real cay_local_plus = .125*(node_grad_y(i ,j ,k ,n)+node_grad_y(i ,j+1,k ,n)+node_grad_y(i+1,j+1,k ,n)+node_grad_y(i+1,j ,k ,n) - +node_grad_y(i ,j ,k+1,n)+node_grad_y(i ,j+1,k+1,n)+node_grad_y(i+1,j+1,k+1,n)+node_grad_y(i+1,j ,k+1,n)); - - Real cay_local_minus = .125*(node_grad_y(i ,j ,k ,n)+node_grad_y(i ,j-1,k ,n)+node_grad_y(i+1,j-1,k ,n)+node_grad_y(i+1,j ,k ,n) - +node_grad_y(i ,j ,k+1,n)+node_grad_y(i ,j-1,k+1,n)+node_grad_y(i+1,j-1,k+1,n)+node_grad_y(i+1,j ,k+1,n)); - - Real caz_local_plus = .125*(node_grad_z(i ,j ,k ,n)+node_grad_z(i ,j+1,k ,n)+node_grad_z(i+1,j+1,k ,n)+node_grad_z(i+1,j ,k ,n) - +node_grad_z(i ,j ,k+1,n)+node_grad_z(i ,j+1,k+1,n)+node_grad_z(i+1,j+1,k+1,n)+node_grad_z(i+1,j ,k+1,n)); - - Real caz_local_minus = .125*(node_grad_z(i ,j ,k ,n)+node_grad_z(i ,j-1,k ,n)+node_grad_z(i+1,j-1,k ,n)+node_grad_z(i+1,j ,k ,n) - +node_grad_z(i ,j ,k+1,n)+node_grad_z(i ,j-1,k+1,n)+node_grad_z(i+1,j-1,k+1,n)+node_grad_z(i+1,j ,k+1,n)); - - - Real cbx_local_plus = .125*(node_grad_x(i ,j ,k ,m)+node_grad_x(i ,j+1,k ,m)+node_grad_x(i+1,j+1,k ,m)+node_grad_x(i+1,j ,k ,m) - +node_grad_x(i ,j ,k+1,m)+node_grad_x(i ,j+1,k+1,m)+node_grad_x(i+1,j+1,k+1,m)+node_grad_x(i+1,j ,k+1,m)); - - Real cbx_local_minus = .125*(node_grad_x(i ,j ,k ,m)+node_grad_x(i ,j-1,k ,m)+node_grad_x(i+1,j-1,k ,m)+node_grad_x(i+1,j ,k ,m) - +node_grad_x(i ,j ,k+1,m)+node_grad_x(i ,j-1,k+1,m)+node_grad_x(i+1,j-1,k+1,m)+node_grad_x(i+1,j ,k+1,m)); - - Real cby_local_plus = .125*(node_grad_y(i ,j ,k ,m)+node_grad_y(i ,j+1,k ,m)+node_grad_y(i+1,j+1,k ,m)+node_grad_y(i+1,j ,k ,m) - +node_grad_y(i ,j ,k+1,m)+node_grad_y(i ,j+1,k+1,m)+node_grad_y(i+1,j+1,k+1,m)+node_grad_y(i+1,j ,k+1,m)); - - Real cby_local_minus = .125*(node_grad_y(i ,j ,k ,m)+node_grad_y(i ,j-1,k ,m)+node_grad_y(i+1,j-1,k ,m)+node_grad_y(i+1,j ,k ,m) - +node_grad_y(i ,j ,k+1,m)+node_grad_y(i ,j-1,k+1,m)+node_grad_y(i+1,j-1,k+1,m)+node_grad_y(i+1,j ,k+1,m)); - - Real cbz_local_plus = .125*(node_grad_z(i ,j ,k ,m)+node_grad_z(i ,j+1,k ,m)+node_grad_z(i+1,j+1,k ,m)+node_grad_z(i+1,j ,k ,m) - +node_grad_z(i ,j ,k+1,m)+node_grad_z(i ,j+1,k+1,m)+node_grad_z(i+1,j+1,k+1,m)+node_grad_z(i+1,j ,k+1,m)); - - Real cbz_local_minus = .125*(node_grad_z(i ,j ,k ,m)+node_grad_z(i ,j-1,k ,m)+node_grad_z(i+1,j-1,k ,m)+node_grad_z(i+1,j ,k ,m) - +node_grad_z(i ,j ,k+1,m)+node_grad_z(i ,j-1,k+1,m)+node_grad_z(i+1,j-1,k+1,m)+node_grad_z(i+1,j ,k+1,m)); - - forcey(i,j,k) += scale_factor * fh_kappa(n,m)*0.5* ( - ((node_grad_y(i+1,j,k,n)*node_grad_x(i+1,j,k,m) + node_grad_y(i+1,j,k+1,n)*node_grad_x(i+1,j,k+1,m)) - -(node_grad_y(i ,j,k,n)*node_grad_x(i ,j,k,m) + node_grad_y(i ,j,k+1,n)*node_grad_x(i ,j,k+1,m)))/dx[0] - +((cay_local_plus *cby_local_plus -caz_local_plus *cbz_local_plus -cax_local_plus *cbx_local_plus) - -(cay_local_minus*cby_local_minus-caz_local_minus*cbz_local_minus-cax_local_minus*cbx_local_minus))/dx[1] - +((node_grad_y(i,j,k+1,n)*node_grad_z(i,j,k+1,m) + node_grad_y(i+1,j,k+1,n)*node_grad_z(i+1,j,k+1,m)) - - (node_grad_y(i,j,k ,n)*node_grad_z(i,j,k ,m) + node_grad_y(i+1,j,k ,n)*node_grad_z(i+1,j,k ,m)))/dx[2]); - - } - } - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for( int n=0 ; n< nspecies; n++){ - for( int m=0 ; m< nspecies; m++){ - - Real cax_local_plus = .125*(node_grad_x(i ,j ,k ,n)+node_grad_x(i ,j+1,k ,n)+node_grad_x(i+1,j+1,k ,n)+node_grad_x(i+1,j ,k ,n) - +node_grad_x(i ,j ,k+1,n)+node_grad_x(i ,j+1,k+1,n)+node_grad_x(i+1,j+1,k+1,n)+node_grad_x(i+1,j ,k+1,n)); - - Real cax_local_minus = .125*(node_grad_x(i ,j ,k ,n)+node_grad_x(i ,j+1,k ,n)+node_grad_x(i+1,j+1,k ,n)+node_grad_x(i+1,j ,k ,n) - +node_grad_x(i ,j ,k-1,n)+node_grad_x(i ,j+1,k-1,n)+node_grad_x(i+1,j+1,k-1,n)+node_grad_x(i+1,j ,k-1,n)); - - Real cay_local_plus = .125*(node_grad_y(i ,j ,k ,n)+node_grad_y(i ,j+1,k ,n)+node_grad_y(i+1,j+1,k ,n)+node_grad_y(i+1,j ,k ,n) - +node_grad_y(i ,j ,k+1,n)+node_grad_y(i ,j+1,k+1,n)+node_grad_y(i+1,j+1,k+1,n)+node_grad_y(i+1,j ,k+1,n)); - - Real cay_local_minus = .125*(node_grad_y(i ,j ,k ,n)+node_grad_y(i ,j+1,k ,n)+node_grad_y(i+1,j+1,k ,n)+node_grad_y(i+1,j ,k ,n) - +node_grad_y(i ,j ,k-1,n)+node_grad_y(i ,j+1,k-1,n)+node_grad_y(i+1,j+1,k-1,n)+node_grad_y(i+1,j ,k-1,n)); - - Real caz_local_plus = .125*(node_grad_z(i ,j ,k ,n)+node_grad_z(i ,j+1,k ,n)+node_grad_z(i+1,j+1,k ,n)+node_grad_z(i+1,j ,k ,n) - +node_grad_z(i ,j ,k+1,n)+node_grad_z(i ,j+1,k+1,n)+node_grad_z(i+1,j+1,k+1,n)+node_grad_z(i+1,j ,k+1,n)); - - Real caz_local_minus = .125*(node_grad_z(i ,j ,k ,n)+node_grad_z(i ,j+1,k ,n)+node_grad_z(i+1,j+1,k ,n)+node_grad_z(i+1,j ,k ,n) - +node_grad_z(i ,j ,k-1,n)+node_grad_z(i ,j+1,k-1,n)+node_grad_z(i+1,j+1,k-1,n)+node_grad_z(i+1,j ,k-1,n)); - - Real cbx_local_plus = .125*(node_grad_x(i ,j ,k ,m)+node_grad_x(i ,j+1,k ,m)+node_grad_x(i+1,j+1,k ,m)+node_grad_x(i+1,j ,k ,m) - +node_grad_x(i ,j ,k+1,m)+node_grad_x(i ,j+1,k+1,m)+node_grad_x(i+1,j+1,k+1,m)+node_grad_x(i+1,j ,k+1,m)); - - Real cbx_local_minus = .125*(node_grad_x(i ,j ,k ,m)+node_grad_x(i ,j+1,k ,m)+node_grad_x(i+1,j+1,k ,m)+node_grad_x(i+1,j ,k ,m) - +node_grad_x(i ,j ,k-1,m)+node_grad_x(i ,j+1,k-1,m)+node_grad_x(i+1,j+1,k-1,m)+node_grad_x(i+1,j ,k-1,m)); - - Real cby_local_plus = .125*(node_grad_y(i ,j ,k ,m)+node_grad_y(i ,j+1,k ,m)+node_grad_y(i+1,j+1,k ,m)+node_grad_y(i+1,j ,k ,m) - +node_grad_y(i ,j ,k+1,m)+node_grad_y(i ,j+1,k+1,m)+node_grad_y(i+1,j+1,k+1,m)+node_grad_y(i+1,j ,k+1,m)); - - Real cby_local_minus = .125*(node_grad_y(i ,j ,k ,m)+node_grad_y(i ,j+1,k ,m)+node_grad_y(i+1,j+1,k ,m)+node_grad_y(i+1,j ,k ,m) - +node_grad_y(i ,j ,k-1,m)+node_grad_y(i ,j+1,k-1,m)+node_grad_y(i+1,j+1,k-1,m)+node_grad_y(i+1,j ,k-1,m)); - - Real cbz_local_plus = .125*(node_grad_z(i ,j ,k ,m)+node_grad_z(i ,j+1,k ,m)+node_grad_z(i+1,j+1,k ,m)+node_grad_z(i+1,j ,k ,m) - +node_grad_z(i ,j ,k+1,m)+node_grad_z(i ,j+1,k+1,m)+node_grad_z(i+1,j+1,k+1,m)+node_grad_z(i+1,j ,k+1,m)); - - Real cbz_local_minus = .125*(node_grad_z(i ,j ,k ,m)+node_grad_z(i ,j+1,k ,m)+node_grad_z(i+1,j+1,k ,m)+node_grad_z(i+1,j ,k ,m) - +node_grad_z(i ,j ,k-1,m)+node_grad_z(i ,j+1,k-1,m)+node_grad_z(i+1,j+1,k-1,m)+node_grad_z(i+1,j ,k-1,m)); - - - - forcez(i,j,k) += scale_factor * fh_kappa(n,m) * 0.5 * ( - ((node_grad_z(i+1,j,k,n)*node_grad_x(i+1,j,k,m) + node_grad_z(i+1,j+1,k,n)*node_grad_x(i+1,j+1,k,m)) - -(node_grad_z(i ,j,k,n)*node_grad_x(i ,j,k,m) + node_grad_z(i ,j+1,k,n)*node_grad_x(i ,j+1,k,m)))/dx[0] - +((caz_local_plus *cbz_local_plus -cay_local_plus *cby_local_plus -cax_local_plus *cbx_local_plus) - -(caz_local_minus*cbz_local_minus-cay_local_minus*cby_local_minus-cax_local_minus*cbx_local_minus))/dx[2] - +((node_grad_z(i,j+1,k,n)*node_grad_y(i,j+1,k,m) + node_grad_z(i+1,j+1,k,n)*node_grad_y(i+1,j+1,k,m)) - - (node_grad_z(i,j ,k,n)*node_grad_y(i,j ,k,m) + node_grad_z(i+1,j ,k,n)*node_grad_y(i+1,j ,k,m)))/dx[1]); - - } - } - }); - -#endif - } - - // set force on walls to be zero since normal velocity is zero - ZeroEdgevalWalls(div_reversible_stress, geom, 0, 1); - -} - -void ComputeDivReversibleStress(std::array& div_reversible_stress, - const MultiFab& rhotot_in, - MultiFab& rho_in, - const Geometry& geom) -{ - BL_PROFILE_VAR("ComputeDivReversibleStress()",ComputeDivReversibleStress); - - BoxArray ba = rho_in.boxArray(); - DistributionMapping dmap = rho_in.DistributionMap(); - const GpuArray dx = geom.CellSizeArray(); - - MultiFab node_grad_c_mf(convert(ba,nodal_flag), dmap, AMREX_SPACEDIM, 1); - - MultiFab conc(ba, dmap, nspecies, 2); - - // rho to conc - VALID REGION ONLY - ConvertRhoCToC(rho_in,rhotot_in,conc,1); - - // fill conc ghost cells - conc.FillBoundary(geom.periodicity()); - MultiFabPhysBC(conc,geom,0,nspecies,SPEC_BC_COMP); - - Real scale_factor = rhobar[0]*k_B*T_init[0]/molmass[0]; - - for ( MFIter mfi(node_grad_c_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - - const Box& bx = mfi.growntilebox(1); - - const Array4& node_grad_c = node_grad_c_mf.array(mfi); - const Array4& c = conc.array(mfi); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { -#if (AMREX_SPACEDIM == 2) - - node_grad_c(i,j,k,0) = (c(i,j,k)-c(i-1,j,k)+c(i,j-1,k)-c(i-1,j-1,k))/(2*dx[0]); - node_grad_c(i,j,k,1) = (c(i,j,k)-c(i,j-1,k)+c(i-1,j,k)-c(i-1,j-1,k))/(2*dx[1]); - -#elif (AMREX_SPACEDIM == 3) - - node_grad_c(i,j,k,0) = (c(i,j,k)-c(i-1,j,k)+c(i,j-1,k)-c(i-1,j-1,k) - +c(i,j,k-1)-c(i-1,j,k-1)+c(i,j-1,k-1)-c(i-1,j-1,k-1))/(4*dx[0]); - node_grad_c(i,j,k,1) = (c(i,j,k)-c(i,j-1,k)+c(i-1,j,k)-c(i-1,j-1,k) - +c(i,j,k-1)-c(i,j-1,k-1)+c(i-1,j,k-1)-c(i-1,j-1,k-1))/(4*dx[1]); - node_grad_c(i,j,k,2) = (c(i,j,k)-c(i,j,k-1)+c(i-1,j,k)-c(i-1,j,k-1) - +c(i,j-1,k)-c(i,j-1,k-1)+c(i-1,j-1,k)-c(i-1,j-1,k-1))/(4*dx[2]); -#endif - }); - } - - for ( MFIter mfi(rhotot_in,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - - const Array4& node_grad_c = node_grad_c_mf.array(mfi); - - AMREX_D_TERM(const Array4 & forcex = div_reversible_stress[0].array(mfi);, - const Array4 & forcey = div_reversible_stress[1].array(mfi);, - const Array4 & forcez = div_reversible_stress[2].array(mfi);); - - AMREX_D_TERM(const Box & bx_x = mfi.nodaltilebox(0);, - const Box & bx_y = mfi.nodaltilebox(1);, - const Box & bx_z = mfi.nodaltilebox(2);); - -#if (AMREX_SPACEDIM == 2) - - amrex::ParallelFor(bx_x, bx_y, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real cx_local_plus = 0.25*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+node_grad_c(i+1,j+1,k,0)+node_grad_c(i+1,j,k,0)); - Real cy_local_plus = 0.25*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i+1,j+1,k,1)+node_grad_c(i+1,j,k,1)); - Real cx_local_minus = 0.25*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+node_grad_c(i-1,j+1,k,0)+node_grad_c(i-1,j,k,0)); - Real cy_local_minus = 0.25*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i-1,j+1,k,1)+node_grad_c(i-1,j,k,1)); - - forcex(i,j,k) = -(node_grad_c(i,j+1,k,0)*node_grad_c(i,j+1,k,1) - node_grad_c(i,j,k,0)*node_grad_c(i,j,k,1))/dx[1] - +(0.5*(cy_local_plus*cy_local_plus-cx_local_plus*cx_local_plus) - -0.5*(cy_local_minus*cy_local_minus-cx_local_minus*cx_local_minus))/(dx[0]); - forcex(i,j,k) = scale_factor * kc_tension*forcex(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real cx_local_plus = 0.25*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+node_grad_c(i+1,j+1,k,0)+node_grad_c(i+1,j,k,0)); - Real cy_local_plus = 0.25*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i+1,j+1,k,1)+node_grad_c(i+1,j,k,1)); - Real cx_local_minus = 0.25*(node_grad_c(i,j,k,0)+node_grad_c(i,j-1,k,0)+node_grad_c(i+1,j-1,k,0)+node_grad_c(i+1,j,k,0)); - Real cy_local_minus = 0.25*(node_grad_c(i,j,k,1)+node_grad_c(i,j-1,k,1)+node_grad_c(i+1,j-1,k,1)+node_grad_c(i+1,j,k,1)); - - forcey(i,j,k) = -(node_grad_c(i+1,j,k,0)*node_grad_c(i+1,j,k,1) - node_grad_c(i,j,k,0)*node_grad_c(i,j,k,1))/dx[0] - +(0.5*(cx_local_plus*cx_local_plus-cy_local_plus*cy_local_plus) - -0.5*(cx_local_minus*cx_local_minus-cy_local_minus*cy_local_minus))/(dx[1]); - forcey(i,j,k) = scale_factor * kc_tension*forcey(i,j,k); - }); - -#elif (AMREX_SPACEDIM == 3) - - amrex::ParallelFor(bx_x, bx_y, bx_z, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real cx_local_plus = (.125)*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+node_grad_c(i+1,j+1,k,0) - +node_grad_c(i+1,j,k,0)+node_grad_c(i+1,j+1,k+1,0)+node_grad_c(i,j,k+1,0) - +node_grad_c(i,j+1,k+1,0)+node_grad_c(i+1,j,k+1,0)); - Real cy_local_plus = (.125)*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i+1,j+1,k,1) - +node_grad_c(i+1,j,k,1)+node_grad_c(i+1,j+1,k+1,1)+node_grad_c(i,j,k+1,1) - +node_grad_c(i,j+1,k+1,1)+node_grad_c(i+1,j,k+1,1)); - Real cz_local_plus = (.125)*(node_grad_c(i,j,k,2)+node_grad_c(i,j+1,k,2)+node_grad_c(i+1,j+1,k,2) - +node_grad_c(i+1,j,k,2)+node_grad_c(i+1,j+1,k+1,2)+node_grad_c(i,j,k+1,2)+ - node_grad_c(i,j+1,k+1,2)+node_grad_c(i+1,j,k+1,2)); - Real cx_local_minus = (.125)*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+node_grad_c(i-1,j+1,k,0)+ - node_grad_c(i-1,j,k,0)+node_grad_c(i,j,k+1,0)+node_grad_c(i,j+1,k+1,0)+ - node_grad_c(i-1,j+1,k+1,0)+node_grad_c(i-1,j,k+1,0)); - Real cy_local_minus = (.125)*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i-1,j+1,k,1)+ - node_grad_c(i-1,j,k,1)+node_grad_c(i,j,k+1,1)+node_grad_c(i,j+1,k+1,1)+ - node_grad_c(i-1,j+1,k+1,1)+node_grad_c(i-1,j,k+1,1)); - Real cz_local_minus = (.125)*(node_grad_c(i,j,k,2)+node_grad_c(i,j+1,k,2)+node_grad_c(i-1,j+1,k,2) - +node_grad_c(i-1,j,k,2)+node_grad_c(i,j,k+1,2)+node_grad_c(i,j+1,k+1,2)+ - node_grad_c(i-1,j+1,k+1,2)+node_grad_c(i-1,j,k+1,2)); - - forcex(i,j,k) = -(.5*(node_grad_c(i,j+1,k,0)*node_grad_c(i,j+1,k,1) - +node_grad_c(i,j+1,k+1,0)*node_grad_c(i,j+1,k+1,1)) - - .5*(node_grad_c(i,j,k,0)*node_grad_c(i,j,k,1)+ - node_grad_c(i,j,k+1,0)*node_grad_c(i,j,k+1,1)))/dx[1] - +(0.5*(cy_local_plus*cy_local_plus+cz_local_plus*cz_local_plus-cx_local_plus*cx_local_plus) - -0.5*(cy_local_minus*cy_local_minus+cz_local_minus*cz_local_minus-cx_local_minus*cx_local_minus))/(dx[0]) - -(.5*(node_grad_c(i,j,k+1,0)*node_grad_c(i,j,k+1,2) - +node_grad_c(i,j+1,k+1,0)*node_grad_c(i,j+1,k+1,2)) - - .5*(node_grad_c(i,j,k,0)*node_grad_c(i,j,k,2)+node_grad_c(i,j+1,k,0)*node_grad_c(i,j+1,k,2)))/dx[2]; - - forcex(i,j,k) = scale_factor * kc_tension*forcex(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real cx_local_plus = (.125)*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+ - node_grad_c(i+1,j+1,k,0)+node_grad_c(i+1,j,k,0)+node_grad_c(i+1,j+1,k+1,0)+ - node_grad_c(i,j,k+1,0)+node_grad_c(i,j+1,k+1,0)+node_grad_c(i+1,j,k+1,0)); - Real cy_local_plus = (.125)*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i+1,j+1,k,1)+ - node_grad_c(i+1,j,k,1)+node_grad_c(i+1,j+1,k+1,1)+ - node_grad_c(i,j,k+1,1)+node_grad_c(i,j+1,k+1,1)+node_grad_c(i+1,j,k+1,1)); - Real cz_local_plus = (.125)*(node_grad_c(i,j,k,2)+node_grad_c(i,j+1,k,2)+ - node_grad_c(i+1,j+1,k,2)+node_grad_c(i+1,j,k,2)+node_grad_c(i+1,j+1,k+1,2)+ - node_grad_c(i,j,k+1,2)+node_grad_c(i,j+1,k+1,2)+node_grad_c(i+1,j,k+1,2)); - - Real cx_local_minus = (.125)*(node_grad_c(i,j,k,0)+node_grad_c(i,j-1,k,0)+ - node_grad_c(i+1,j-1,k,0)+node_grad_c(i+1,j,k,0)+node_grad_c(i,j,k+1,0)+ - node_grad_c(i,j-1,k+1,0)+node_grad_c(i+1,j-1,k+1,0)+node_grad_c(i+1,j,k+1,0)); - Real cy_local_minus = (.125)*(node_grad_c(i,j,k,1)+node_grad_c(i,j-1,k,1)+ - node_grad_c(i+1,j-1,k,1)+node_grad_c(i+1,j,k,1)+node_grad_c(i,j,k+1,1)+ - node_grad_c(i,j-1,k+1,1)+node_grad_c(i+1,j-1,k+1,1)+node_grad_c(i+1,j,k+1,1)); - Real cz_local_minus = (.125)*(node_grad_c(i,j,k,2)+node_grad_c(i,j-1,k,2)+ - node_grad_c(i+1,j-1,k,2)+node_grad_c(i+1,j,k,2)+node_grad_c(i,j,k+1,2)+ - node_grad_c(i,j-1,k+1,2)+node_grad_c(i+1,j-1,k+1,2)+node_grad_c(i+1,j,k+1,2)); - - forcey(i,j,k) = -(.5*(node_grad_c(i+1,j,k,0)*node_grad_c(i+1,j,k,1) - +node_grad_c(i+1,j,k+1,0)*node_grad_c(i+1,j,k+1,1)) - - .5*(node_grad_c(i,j,k,0)*node_grad_c(i,j,k,1) - +node_grad_c(i,j,k+1,0)*node_grad_c(i,j,k+1,1)))/dx[0] - +(0.5*(cx_local_plus*cx_local_plus+cz_local_plus*cz_local_plus-cy_local_plus*cy_local_plus) - -0.5*(cx_local_minus*cx_local_minus+cz_local_minus*cz_local_minus-cy_local_minus*cy_local_minus))/(dx[1]) - -(.5*(node_grad_c(i,j,k+1,1)*node_grad_c(i,j,k+1,2) - +node_grad_c(i+1,j,k+1,1)*node_grad_c(i+1,j,k+1,2)) - - .5*(node_grad_c(i,j,k,1)*node_grad_c(i,j,k,2)+node_grad_c(i+1,j,k,1)*node_grad_c(i+1,j,k,2)))/dx[2]; - - forcey(i,j,k) = scale_factor * kc_tension*forcey(i,j,k); - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real cx_local_plus = (.125)*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+ - node_grad_c(i+1,j+1,k,0)+node_grad_c(i+1,j,k,0)+node_grad_c(i+1,j+1,k+1,0)+ - node_grad_c(i,j,k+1,0)+node_grad_c(i,j+1,k+1,0)+node_grad_c(i+1,j,k+1,0)); - Real cy_local_plus = (.125)*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i+1,j+1,k,1)+ - node_grad_c(i+1,j,k,1)+node_grad_c(i+1,j+1,k+1,1)+node_grad_c(i,j,k+1,1)+ - node_grad_c(i,j+1,k+1,1)+node_grad_c(i+1,j,k+1,1)); - Real cz_local_plus = (.125)*(node_grad_c(i,j,k,2)+node_grad_c(i,j+1,k,2)+node_grad_c(i+1,j+1,k,2)+ - node_grad_c(i+1,j,k,2)+node_grad_c(i+1,j+1,k+1,2)+node_grad_c(i,j,k+1,2)+ - node_grad_c(i,j+1,k+1,2)+node_grad_c(i+1,j,k+1,2)); - - Real cx_local_minus = (.125)*(node_grad_c(i,j,k,0)+node_grad_c(i,j+1,k,0)+node_grad_c(i+1,j+1,k,0)+ - node_grad_c(i+1,j,k,0)+node_grad_c(i,j,k-1,0)+node_grad_c(i,j+1,k-1,0) - +node_grad_c(i+1,j+1,k-1,0)+node_grad_c(i+1,j,k-1,0)); - Real cy_local_minus = (.125)*(node_grad_c(i,j,k,1)+node_grad_c(i,j+1,k,1)+node_grad_c(i+1,j+1,k,1) - +node_grad_c(i+1,j,k,1)+node_grad_c(i,j,k-1,1)+node_grad_c(i,j+1,k-1,1) - +node_grad_c(i+1,j+1,k-1,1)+node_grad_c(i+1,j,k-1,1)); - Real cz_local_minus = (.125)*(node_grad_c(i,j,k,2)+node_grad_c(i,j+1,k,2)+node_grad_c(i+1,j+1,k,2) - +node_grad_c(i+1,j,k,2)+node_grad_c(i,j,k-1,2)+node_grad_c(i,j+1,k-1,2) - +node_grad_c(i+1,j+1,k-1,2)+node_grad_c(i+1,j,k-1,2)); - - forcez(i,j,k) = -(.5*(node_grad_c(i+1,j,k,0)*node_grad_c(i+1,j,k,2) - +(node_grad_c(i+1,j+1,k,0)*node_grad_c(i+1,j+1,k,2))) - - .5*(node_grad_c(i,j,k,0)*node_grad_c(i,j,k,2) - +node_grad_c(i,j+1,k,0)*node_grad_c(i,j+1,k,2)))/dx[0] - +(0.5*(cx_local_plus*cx_local_plus+cy_local_plus*cy_local_plus-cz_local_plus*cz_local_plus) - -0.5*(cx_local_minus*cx_local_minus+cy_local_minus*cy_local_minus-cz_local_minus*cz_local_minus))/(dx[2]) - -(.5*(node_grad_c(i,j+1,k,1)*node_grad_c(i,j+1,k,2) - +node_grad_c(i+1,j+1,k,1)*node_grad_c(i+1,j+1,k,2)) - - .5*(node_grad_c(i,j,k,1)*node_grad_c(i,j,k,2)+node_grad_c(i+1,j,k,1)*node_grad_c(i+1,j,k,2)))/dx[1]; - - forcez(i,j,k) = scale_factor * kc_tension*forcez(i,j,k); - }); - -#endif - } - - // set force on walls to be zero since normal velocity is zero - ZeroEdgevalWalls(div_reversible_stress, geom, 0, 1); - -} diff --git a/src_multispec/multispec_functions.H~ b/src_multispec/multispec_functions.H~ deleted file mode 100644 index 6d92f2e4..00000000 --- a/src_multispec/multispec_functions.H~ +++ /dev/null @@ -1,1605 +0,0 @@ -#ifndef _multispec_functions_H_ -#define _multispec_functions_H_ - -#include - -#include "multispec_namespace.H" - -#include "common_functions.H" - -#include "StochMassFlux.H" -#include "StochMomFlux.H" - -#include - -using namespace multispec; -using namespace amrex; - -///////////////////////////////////////////////////////////////////////////////// -// in multispec_functions.cpp - -void InitializeMultispecNamespace(); - -///////////////////////////////////////////////////////////////////////////////// -// in ComputeDivReversibleStress.cpp -void ComputeDivFHReversibleStress(std::array& div_reversible_stress, - const MultiFab& rhotot_in, - MultiFab& rho_in, - const Geometry& geom); - - -///////////////////////////////////////////////////////////////////////////////// -// in ComputeDivReversibleStress.cpp -void ComputeDivReversibleStress(std::array& div_reversible_stress, - const MultiFab& rhotot_in, - MultiFab& rho_in, - const Geometry& geom); - -///////////////////////////////////////////////////////////////////////////////// -// in ComputeMassFluxdiv.cpp - -void ComputeMassFluxdiv(MultiFab& rho, - MultiFab& rhotot, - const MultiFab& Temp, - MultiFab& diff_mass_fluxdiv, - MultiFab& stoch_mass_fluxdiv, - std::array& diff_mass_flux, - std::array& stoch_mass_flux, - StochMassFlux& sMassFlux, - const Real& dt, const Real& stage_time, const Geometry& geom, - Vector& weights, - MultiFab& charge, - std::array& grad_Epot, - MultiFab& Epot, - MultiFab& permittivity, - const int& zero_initial_Epot=1); - -void ComputeHigherOrderTerm(MultiFab& molarconc, - std::array& diff_mass_flux, - const Geometry& geom); - -void ComputeFHHigherOrderTerm(MultiFab& molarconc, - std::array& diff_mass_flux, - const Geometry& geom); - -void MultiFabPhysBCFH(MultiFab& phi, const Geometry& geom, int scomp, int ncomp, - const Real& scale, const Real& time = 0.); - -///////////////////////////////////////////////////////////////////////////////// -// in DiffusiveMassFlux.cpp - -void DiffusiveMassFluxdiv(const MultiFab& rho, - const MultiFab& rhotot, - MultiFab& molarconc, - const MultiFab& rhoWchi, - const MultiFab& Gamma, - MultiFab& diff_mass_fluxdiv, - std::array< MultiFab, AMREX_SPACEDIM >& diff_mass_flux, - const Geometry& geom); - -void DiffusiveMassFlux(const MultiFab& rho, - const MultiFab& rhotot, - MultiFab& molarconc, - const MultiFab& rhoWchi, - const MultiFab& Gamma, - std::array< MultiFab, AMREX_SPACEDIM >& diff_mass_flux, - const Geometry& geom); - -///////////////////////////////////////////////////////////////////////////////// -// in ComputeMixtureProperties.cpp -void ComputeMixtureProperties(const MultiFab& rho, - const MultiFab& rhotot, - MultiFab& D_bar, - MultiFab& D_therm, - MultiFab& Hessian); - -void ComputeEta(const MultiFab& rho_in, - const MultiFab& rhotot_in, - MultiFab& eta_in); - -///////////////////////////////////////////////////////////////////////////////// -// in CorrectionFlux.cpp - -void CorrectionFlux(const MultiFab& rho, const MultiFab& rhotot, - std::array< MultiFab, AMREX_SPACEDIM >& flux); - -///////////////////////////////////////////////////////////////////////////////// -// in CorrectionFlux.cpp - -void CorrectionFlux(const MultiFab& rho, const MultiFab& rhotot, - std::array< MultiFab, AMREX_SPACEDIM >& flux); - -///////////////////////////////////////////////////////////////////////////////// -// in ElectroDiffusiveMassFluxdiv.cpp - -void ElectroDiffusiveMassFluxdiv(const MultiFab& rho, - const MultiFab& Temp, - const MultiFab& rhoWchi, - std::array< MultiFab, AMREX_SPACEDIM >& diff_mass_flux, - MultiFab& diff_mass_fluxdiv, - std::array< MultiFab, AMREX_SPACEDIM >& stoch_mass_flux, - MultiFab& charge, - std::array< MultiFab, AMREX_SPACEDIM >& grad_Epot, - MultiFab& Epot, - const MultiFab& permittivity, - Real dt, - int zero_initial_Epot, - const Geometry& geom); - -void ElectroDiffusiveMassFlux(const MultiFab& rho, - const MultiFab& Temp, - const MultiFab& rhoWchi, - std::array< MultiFab, AMREX_SPACEDIM >& electro_mass_flux, - std::array< MultiFab, AMREX_SPACEDIM >& diff_mass_flux, - std::array< MultiFab, AMREX_SPACEDIM >& stoch_mass_flux, - MultiFab& charge, - std::array< MultiFab, AMREX_SPACEDIM >& grad_Epot, - MultiFab& Epot, - const MultiFab& permittivity, - Real dt, - int zero_initial_Epot, - const Geometry& geom); - -void LimitEMF(const MultiFab& rho_in, - std::array< MultiFab, AMREX_SPACEDIM >& electro_mass_flux); - -///////////////////////////////////////////////////////////////////////////////// -// in FluicCharge.cpp - -void DotWithZ(const MultiFab& mf, - MultiFab& mfdotz, - int abs_z=0); - -void DotWithZFace(std::array< const MultiFab, AMREX_SPACEDIM >& mf, - std::array< MultiFab, AMREX_SPACEDIM >& mfdotz, - int abs_0); - -void ComputeChargeCoef(const MultiFab& rho_in, - const MultiFab& Temp_in, - MultiFab& charge_coef_in); - -void EnforceChargeNeutrality(); - -void ImplicitPotentialCoef(); - -void ModifyS(); - -void ComputePermittivity(); - -void ComputeLorentzForce(std::array< MultiFab, AMREX_SPACEDIM >& Lorentz_force, - std::array< MultiFab, AMREX_SPACEDIM >& grad_Epot, - const MultiFab& permittivity, - const MultiFab& charge, - const Geometry& geom); - -void ComputeE_ext(std::array< MultiFab, AMREX_SPACEDIM >& E_ext); - -void ZeroEpsOnWall(std::array< MultiFab, AMREX_SPACEDIM >& beta); - -///////////////////////////////////////////////////////////////////////////////// -// in InitialProjection.cpp - -void InitialProjection(std::array< MultiFab, AMREX_SPACEDIM >& umac, - MultiFab& rho, MultiFab& rhotot, - MultiFab& diff_mass_fluxdiv, - MultiFab& stoch_mass_fluxdiv, - std::array< MultiFab, AMREX_SPACEDIM >& stoch_mass_flux, - StochMassFlux& sMassFlux, - const MultiFab& Temp, const MultiFab& eta, - const std::array< MultiFab, NUM_EDGE >& eta_ed, - const Real& dt, const Real& time, const Geometry& geom, - MultiFab& charge_old, - std::array& grad_Epot_old, - MultiFab& Epot, - MultiFab& permittivity); - -///////////////////////////////////////////////////////////////////////////////// -// in MassFluxUtil.cpp - -void ComputeMolconcMolmtot(const MultiFab& rho, - const MultiFab& rhotot, - MultiFab& molarconc, - MultiFab& molmtot); - -void ComputeMassfrac(const MultiFab& rho, - const MultiFab& rhotot, - MultiFab& massfrac); - -void ComputeGamma(const MultiFab& molarconc, - const MultiFab& Hessian, - MultiFab& Gamma); - -void ComputeFHGamma(const MultiFab& massfrac, - MultiFab& Gamma); - -void ComputeRhoWChi(const MultiFab& rho, - const MultiFab& rhotot, - const MultiFab& molarconc, - MultiFab& rhoWchi, - const MultiFab& D_bar); - -void ComputeZetaByTemp(const MultiFab& molarconc, - const MultiFab& D_bar, - const MultiFab& Temp, - MultiFab& zeta_by_Temp, - const MultiFab& D_therm); - -void ComputeSqrtLonsagerFC(const MultiFab& rho, - const MultiFab& rhotot, - std::array< MultiFab, AMREX_SPACEDIM >& sqrtLonsager_fc, - const Geometry& geom); - -///////////////////////////////////////////////////////////////////////////////// -// in MatvecMul.cpp - -void MatvecMul(MultiFab& x, - const MultiFab& A); - -///////////////////////////////////////////////////////////////////////////////// -// in MkDiffusiveMFluxdiv.cpp - -void MkDiffusiveMFluxdiv(std::array & m_update, - const std::array & umac, - const MultiFab& eta, - const std::array & eta_ed, - const MultiFab& kappa, - const Geometry& geom, - const Real* dx, - const int& increment); - -///////////////////////////////////////////////////////////////////////////////// -// in ProjectOntoEOS.cpp - -void ProjectOntoEOS(MultiFab& rho_in); - -///////////////////////////////////////////////////////////////////////////////// -// in RhoUtil.cpp - -void RhototBCInit(); - -void ComputeRhotot(const MultiFab& rho, MultiFab& rhotot, int include_ghost=0); - -void ConvertRhoCToC(MultiFab& rho, const MultiFab& rhotot, MultiFab& conc, int rho_to_c); - -void FillRhoRhototGhost(MultiFab& rho, MultiFab& rhotot, const Geometry& geom); - -void FillRhototGhost(MultiFab& rhotot_in, const MultiFab& conc_in, const Geometry& geom); - -///////////////////////////////////////////////////////////////////////////////// -// Device calls from MassFluxUtil.cpp - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeChi(int nspecies_in, - GpuArray& molmass_in, - GpuArray& rhoN, - Real rhotot, - GpuArray& MolarConcN, - Array2D& chi, - Array2D& D_barN); - -/** - * \brief Compute total molar mass - * - * \param[in] W Mass fractions. - * \param[in] molmass_in - * \param[in] nspecies_in Number of species - * \param[out] molmtot Total molar mass - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void GetMolTot ( GpuArray& W, - GpuArray& molmass_in, - Real& molmtot_in, - int nspecies_in) -{ - Real sumWOverN = 0.0; - - for (int n=0; n& W, - GpuArray& monomer_in, - Real& monomertot, - int nspecies_in) -{ - Real sumWOverN = 0.0; - - for (int n=0; n& molmass_in, - GpuArray& RhoN, - Real rhotot_in, - GpuArray& MolarConcN, - Real& molmtot_in){ - - GpuArray W; /**< Mass fraction w_i = rho_i/rho */ - - // calculate mass fraction and total molar mass (1/m=Sum(w_i/m_i)) - for (int n=0; n nspecies_in - W[n] = RhoN[n] / rhotot_in; - } - - GetMolTot(W, molmass_in, molmtot_in, nspecies_in); - - // calculate molar concentrations in each cell (x_i=m*w_i/m_i) - for (int n=0; n& matrixIn) { -//HACK:required - - for (int n=1; n<=size; ++n ){ - for (int m=1; m<=size; ++m ){ - matrixIn(n,m) = 0.0; - } - } -} - -/** - * \brief Populate X_xxT - * - * \param[out] X_xxT - * \param[in] MolarConcN - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void PopulateX_xxT ( - Array2D& X_xxT, - GpuArray& MolarConcN) { - - for (int row=1; row<=nspecies;++row ){ - // diagonal entries - X_xxT(row,row) = MolarConcN[row-1] - MolarConcN[row-1]*MolarConcN[row-1]; - for (int column=1; column<=row-1; ++column ){ - // off-diagonal entries - // form x*traspose(x) off diagonals - X_xxT(row,column) = -MolarConcN[row-1]*MolarConcN[column-1]; - //symmetric - X_xxT(column,row) = X_xxT(row,column); - } - } -} - -/** - * \brief Compute Gamma = I + matmul(X_xxT, Hessian). - * - * Adds the identity matrix during initialization. Then adds product of X_xxT and Hessian matrices. - * - * \param[out] GammaN - * \param[in] X_xxT - * \param[in] HessianN - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void GammaIPlusMatmul ( - Array2D& GammaN, - Array2D& X_xxT, - Array2D& HessianN) { - - //Compute Gamma - for (int row=1; row<=nspecies; ++row){ - for (int column=1; column<=nspecies; ++column){ - if (row == column) { - GammaN(row,column) = 1.0; // add the identity matrix - } else { - GammaN(row,column) = 0.0; // initialize off-diagonals to 0 - } - for (int n=1; n<=nspecies; ++n){ - GammaN(row, column) += X_xxT(row,n) * HessianN(n,column); - } - } - } -} - -/** - * \brief Compute Gamma - * - * \param[in] MolarConcN molar concentration - * \param[in] HessianN - * \param[in,out] GammaN - * - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeGammaLocal ( - GpuArray& MolarConcN, - Array2D& HessianN, - Array2D& GammaN){ - - //Local variables - Array2D X_xxT; - - if ((use_multiphase == 1) && (nspecies == 2)){ - - // local to define Gamma matrix - Real w1 = MolarConcN[0]; - Real w2 = MolarConcN[1]; - - if (std::abs(w1+w2-1.0) > 1e-14){ - printf("Mole fractions do not add up in gamma computation\n"); - } - if (w1 < 0){ - w1 = 0.0; - w2 = 1.0; - } - if (w2 < 0){ - w1 = 1.0; - w2 = 0.0; - } - - if( n_gex == 1 ){ - - GammaN(1,2) = w1 * n_gex * n_gex * alpha_gex * std::pow(w1,n_gex-1) * std::pow(w2,n_gex-1); - GammaN(2,1) = w2 * n_gex * n_gex * alpha_gex * std::pow(w2,n_gex-1) * std::pow(w1,n_gex-1); - GammaN(1,1) = 1.0 ; - GammaN(2,2) = 1.0 ; - - } else { - - GammaN(1,2) = w1 * n_gex * n_gex * alpha_gex * std::pow(w1,n_gex-1) * std::pow(w2,n_gex-1); - GammaN(2,1) = w2 * n_gex * n_gex * alpha_gex * std::pow(w2,n_gex-1) * std::pow(w1,n_gex-1); - GammaN(1,1) = 1.0 + w1 * n_gex * (n_gex-1) * alpha_gex * std::pow(w1,n_gex-2) * std::pow(w2,n_gex); - GammaN(2,2) = 1.0 + w2 * n_gex * (n_gex-1) * alpha_gex * std::pow(w2,n_gex-2) * std::pow(w1,n_gex); - - } - - } else { - - //populate X_xxT - if (is_ideal_mixture == 1){ - MatrixToZeros(nspecies, X_xxT); - } else { - PopulateX_xxT(X_xxT, MolarConcN); - } - - //Compute Gamma - GammaIPlusMatmul(GammaN, X_xxT, HessianN); - - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeFHGammaLocal ( - GpuArray& massFrac, - Array2D& GammaN){ - - Real Nbar; - - GetMonomerTot ( massFrac, fh_monomers, Nbar, nspecies); - - for (int n=0; n dest, - GpuArray& molmass_in, - GpuArray& rhoN, - Array2D& D_barN, - Array2D& rhoWchiN){ - - GpuArray molmass_sub, rhoN_sub, MolarConcN_sub; - Array1D W_sub; - Array2D D_barN_sub, chi_sub; - Real rhotot_sub, molmtot_sub; - Real Deff, tmp; - - // create a vector of non-trace densities and molmass for the subsystem - for (int row=1; row<=nspecies; ++row){ - if ( dest(row) != 0 ){ - molmass_sub[dest(row)-1] = molmass_in[row-1]; - rhoN_sub[dest(row)-1] = rhoN[row-1]; - } - } - - // renormalize total density and mass fractions - rhotot_sub = 0.0; - for (int row=1; row<=nspecies_sub; ++row){ - rhotot_sub += rhoN_sub[row-1]; - } - - for (int row=1; row<=nspecies_sub; ++row){ - W_sub(row) = rhoN_sub[row-1]/rhotot_sub; - } - - // construct D_bar_sub by mapping the full D_bar into D_bar_sub - // you could read in only the lower diagonals, - // reflect, and set the diagnals to zero if you want - - for (int row=1; row<=nspecies; ++row){ - if (dest(row) == 0){ - continue; - } - for (int column=1; column<=nspecies; ++column){ - if (dest(column) != 0){ - D_barN_sub(dest(row),dest(column)) = D_barN(row,column); - } - } - } - - // compute molarconc_sub and molmtot_sub - ComputeMolconcMolmtotLocal(nspecies_sub, molmass_sub, rhoN_sub, rhotot_sub, MolarConcN_sub, molmtot_sub); - - // compute chi_sub -// amrex::Print() << " entering reduced system " << nspecies_sub << std::endl; -// for (int n=1; n<=nspecies_sub; ++n){ -// amrex::Print() << "n , dest,molmass rho massf " << n << " " << dest(n) << " " << molmass_sub[n-1] << " " << -// rhoN_sub[n-1] << " " << MolarConcN_sub[n-1] << std::endl; -// for (int m=1; m<=nspecies_sub; ++m){ -// amrex::Print() << "m d chi " << m << " " << D_barN_sub(n,m) << " " << chi_sub(n,m) << std::endl; -// } -// } - ComputeChi(nspecies_sub, molmass_sub, rhoN_sub, rhotot_sub, MolarConcN_sub, chi_sub, D_barN_sub); - - - // compute full rho*W*chi - MatrixToZeros(nspecies, rhoWchiN); - for (int column=1; column<=nspecies; ++column){ - if (dest(column) == 0){ //column of trace species - // compute Deff - Deff = 0.0; - for (int k=1; k<=nspecies; ++k){ - if (dest(k) != 0){ - Deff = Deff + MolarConcN_sub[dest(k)-1]/D_barN(k,column); - } - } - Deff = 1.0/Deff; - - // assign rhoWchi - for (int row=1; row<=nspecies; ++row){ - if ( row == column ){ - rhoWchiN(row,column) = rhotot_sub*Deff*molmass_in[row-1]/molmtot_sub; - } else if ( dest(row) == 0 ){ - rhoWchiN(row,column) = 0.0; - } else { - tmp = 0.0; - for (int k=1; k<=nspecies; ++k){ - if ( dest(k)!=0 ){ - tmp = tmp + chi_sub(dest(row),dest(k))*MolarConcN_sub[dest(k)-1]/D_barN(k,column); - } - } - - - rhoWchiN(row,column) = Deff*rhoN_sub[dest(row)-1]*(tmp - molmass_in[column-1]/molmtot_sub); - } - } - } else { // column of non-trace species - // assign rhoWchi - for (int row=1; row<=nspecies; ++row){ - if ( dest(row) == 0 ){ - rhoWchiN(row,column) = 0.0; - } else { - rhoWchiN(row,column) = rhoN_sub[dest(row)-1]*chi_sub(dest(row),dest(column)); - } - } - } - } -} - - -/** - * \brief Set 2-D square Matrix = 2-D squareMatrix - * - * \param[out] Matrix_out - * \param[in] Matrix_in - * \param[in] size - * Takes two Array2Ds of size `size` and sets Matrix_out = Matrix_in - * - */ - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void Matrix2DSetEqual ( - Array2D& Matrix_out, - Array2D& Matrix_in, - int size){ - - for (int i=1; i<=size; ++i){ - for (int j=1; j<=size; ++j){ - Matrix_out(i,j) = Matrix_in(i,j); - } - } - -} - - -/** - * \brief D_bar to chi - iterative - * - * \param[in] D_bar matrix of Maxwell-Stefan binary diffusion coefficient - * \param[in] Xk mole fractions --- MUST NOT BE ZERO - * \param[in] molmass_local - * \param[out] chi multispecies diffusion matrix - * - * The number of terms in the sum is determined by the global variable `chi_iterations`. - * 3-5 is a reasonable range for the value. - * - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void Dbar2chiIterative ( int nspecies_in, - Array2D& D_bar, - GpuArray& Xk, - GpuArray& molmass_local, - Array2D& chi){ - - Array1D Xkp, Ykp, Di, Mmat, Minv; - Array2D Pmat, Deltamat, Zmat, Jmat, PJ, matrix1, matrix2; - - // mole fractions correction - Real sum = 0.; - for (int ii=1; ii<=nspecies_in; ++ii){ - Xkp(ii) = Xk[ii-1]; - sum += Xkp(ii); - } -// if(amrex::abs(sum-1.) > 1.e-12){ -// amrex::Print() << " bad x in dbarier " << sum << std::endl; -// } - - // molecular weight of mixture - EGLIB - Real MWmix = 0.0; - for (int ii=1; ii<=nspecies_in; ++ii){ - MWmix = MWmix + Xkp(ii)*molmass_local[ii-1]; - - } - - // mass fractions correction - EGLIB - for (int ii=1; ii<=nspecies_in; ++ii){ - Ykp(ii) = molmass_local[ii-1] / MWmix * Xkp(ii); - } - - // Find Di matrix - Real term2; - for (int i=1; i<=nspecies_in; ++i){ - term2 = 0.0; - for (int j=1; j<=nspecies_in; ++j){ - if (j!=i){ - term2 = term2 + Xkp(j)/D_bar(i,j); - } - } - Di(i) = (1.0 - Ykp(i))/term2; - } - - // Compute Mmat and Minv - for (int i=1; i<=nspecies_in; ++i){ - Mmat(i) = Xkp(i)/Di(i); - Minv(i) = Di(i)/Xkp(i); - } - - // Compute P matrix - MatrixToZeros(nspecies_in, Pmat); - for (int i=1; i<=nspecies_in; ++i){ - for (int j=1; j<=nspecies_in; ++j){ - Pmat(i,j) = - Ykp(j); - if (i==j){ - Pmat(i,j) = Pmat(i,j) + 1.0; - } - } - } - - // Compute Deltamat - MatrixToZeros(nspecies_in, Deltamat); - Real term1; - - for (int i=1; i<=nspecies_in; ++i){ - for (int j=1; j<=nspecies_in; ++j){ - if (i==j){ - term1 = 0.0; - for (int k=1; k<=nspecies_in; ++k){ - if (k!=i){ - term1 = term1 + Xkp(i)*Xkp(k)/D_bar(i,k); - } - } - Deltamat(i,i) = term1; - } else { - Deltamat(i,j) = -Xkp(i)*Xkp(j)/D_bar(i,j); - } - Zmat(i,j) = -Deltamat(i,j); - } - } - - // Compute Zmat - for (int i=1; i<=nspecies_in; ++i){ - Zmat(i,i) = Zmat(i,i) + Mmat(i); - } - - // Compute Jmat - for (int i=1; i<=nspecies_in; ++i){ - for (int j=1; j<=nspecies_in; ++j){ - Jmat(i,j) = Minv(i)*Zmat(i,j); - } - } - - // Compute PJ - MatrixToZeros(nspecies_in, PJ); - for (int i=1; i<=nspecies_in; ++i){ - for (int j=1; j<=nspecies_in; ++j){ - for (int k=1; k<=nspecies_in; ++k){ - PJ(i,j) = PJ(i,j) + Pmat(i,k)*Jmat(k,j); - } - } - } - - // Compute P M^-1 Pt; store it in matrix2 - Real scr; - for (int i=1; i<=nspecies_in; ++i){ - for (int j=1; j<=nspecies_in; ++j){ - scr = 0.0; - for (int k=1; k<=nspecies_in; ++k){ - // notice the change in indices for Pmat to represent Pmat^t - scr = scr + Pmat(i,k)*Minv(k)*Pmat(j,k); - } - matrix2(i,j) = scr; - chi(i,j) = scr; - } - } - - - for (int jj=1; jj<=chi_iterations; ++jj){ - for (int i=1; i<=nspecies_in; ++i){ - for (int j=1; j<=nspecies_in; ++j){ - scr = 0.0; - for (int k=1; k<=nspecies_in; ++k){ - scr = scr + PJ(i,k)*chi(k,j); - } - matrix1(i,j) = scr + matrix2(i,j); - } - } - Matrix2DSetEqual(chi, matrix1, MAX_SPECIES); - } -} - - - - - - -/** - * \brief Compute Chi - * - * \param[in] nspecies_in Number of species - * \param[in] molmass_in - * \param[in] rhoN - * \param[in] rhotot - * \param[in] MolarConcN - * \param[in,out] chi multispcies diffusion matrix. - * \param[in] D_barN matrix of Maxwell-Stefan binary diffusion coefficient - * - * - */ - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeChi (int nspecies_in, - GpuArray& molmass_in, - GpuArray& rhoN, - Real rhotot, - GpuArray& MolarConcN, - Array2D& chi, - Array2D& D_barN){ - - //local variables - Real eepsilon=1.e-16; - - Array1D W; - - /* - * if nspecies_in = 2, use analytic formulas - * note: nspecies_in = 1 (that is, ntrace = nspecies-1) is treated separately (chi=0) - * before this routine is called - */ - - if ( nspecies_in == 2){ - - W(1) = rhoN[0]/rhotot; - W(2) = rhoN[1]/rhotot; - - if ( use_multiphase == 1 ){ - W(1) = std::max(std::min(W(1),1.0),eepsilon); - W(2) = std::max(std::min(W(2),1.0),eepsilon); - } - - Real tmp = molmass_in[0]*W(2) + molmass_in[1]*W(1); - tmp = D_barN(1,2)*tmp*tmp/molmass_in[0]/molmass_in[1]; - - chi(1,1) = tmp*W(2)/W(1); - chi(1,2) = -tmp; - chi(2,1) = -tmp; - chi(2,2) = tmp*W(1)/W(2); - - return; - } - - // compute chi either selecting inverse/pseudoinverse or iterative methods - if (use_lapack == 1){ - printf("Compute Chi: use_lapack not supported\n"); - //call amrex_error('compute_chi: use_lapack not supported') - } - - Dbar2chiIterative(nspecies_in, D_barN,MolarConcN,molmass_in,chi); - -} - - - -/** - * \brief Compute rhoWchi - * - * \param[in] rhoN densities - * \param[in] rhotot total density - * \param[in] MolarConcN molar concentration - * \param[in,out] rhoWchiN - * \param[in] D_barN MS diff-coefs - * \param[in] molmass_in - * - */ - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeRhoWChiLocal ( - GpuArray& rhoN, - Real rhotot, - GpuArray& MolarConcN, - Array2D& rhoWchiN, - Array2D& D_barN, - GpuArray& molmass_in){ - - //local variables - Array1D W; - Array2D chi; - int ntrace; /**< number of trace species with w_k < fractional_tolerance */ - int nspecies_sub; /**< dim of subsystem = nspecies - ntrace */ - Real w1,w2; - - /* - * this is a mapping used to eliminate the elements in D_bar that we don't need (for D_bar_sub) - * and for expanding sqrtLonsager_sub into sqrtLonsager - * it will contain the numbers 1, 2, ..., (nspecies-ntrace) - * with zeros in elements corresponding to trace elements - * (example) for a 5-species system having trace species 2 and 5: - * species 1 2 3 4 5 - * dest(species) 1 0 2 3 0 - */ - Array1D dest; - - if ((use_multiphase == 1) && (nspecies == 2)) { - - w1 = MolarConcN[0]; - w2 = MolarConcN[1]; - if (w1<0){ - w1=0.0; - w2=1.0; - } - if (w2<0){ - w1=1.0; - w2=0.0; - } - rhoWchiN(1,1) = w2*rho0*D_barN(1,2); - rhoWchiN(1,2) = -w1*rho0*D_barN(1,2); - rhoWchiN(2,1) = -w2*rho0*D_barN(1,2); - rhoWchiN(2,2) = w1*rho0*D_barN(1,2); - - } else { - - // compute the number of trace species - // build the mapping for expanding/contracting arrays - ntrace = 0; - for (int row=1; row<=nspecies; ++row){ - W(row) = rhoN[row-1] / rhotot; - if (W(row) < fraction_tolerance){ - ntrace = ntrace + 1; - dest(row) = 0; - } else { - dest(row) = row - ntrace; - } - } - - if (ntrace == nspecies -1){ - // this is all trace species except for 1 (essentially pure solvent); - MatrixToZeros(nspecies, rhoWchiN); //set rhoWchi to zero - - } else if (ntrace == 0){ - // there are no trace species, hence chi = chi_sub - - ComputeChi(nspecies, molmass_in, rhoN, rhotot, MolarConcN, chi, D_barN); - - for (int row=1; row<=nspecies; ++row){ - for (int column=1; column<=nspecies; ++column){ - rhoWchiN(row, column) = rhoN[row-1] * chi(row,column); - } - } - } - else { - // if there are trace species, we consider a subsystem - // consisting of non-trace species - nspecies_sub = nspecies - ntrace; - -// amrex::Print() << " entering subspecies " << ntrace << " " << nspecies_sub << std::endl; -// for (int n=1; n<=nspecies; ++n){ -// amrex::Print() << "n , dest,molmass rho " << n << " " << dest(n) << " " << molmass_in[n-1] << " " << rhoN[n-1] << std::endl; -// for (int m=1; m<=nspecies; ++m){ -// amrex::Print() << "m d rhowchi " << m << " " << D_barN(n,m) << " " << rhoWchiN(n,m) << std::endl; -// } -// } - ComputeChiSub(nspecies_sub, - dest, - molmass_in, - rhoN, - D_barN, - rhoWchiN); - - } - } -} - -/** - * \brief Compute nonnegative rho face-centered average - * - * \param[in] Rho1 density from a neighboring cell - * \param[in] Rho2 density from a neighboring cell - * \param[in] dx - * \param[out] RhoAv face-centered average - * - */ - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeNonnegativeRhoAv ( - GpuArray& Rho1, - GpuArray& Rho2, - const GpuArray& dx, - GpuArray& molmass_in, - GpuArray& RhoAv){ - - Real dv, value1, value2, tmp1, tmp2; - -#if (AMREX_SPACEDIM == 2) - dv = 1.0; - for (int dim=0; dim<2; ++dim){ - dv *= dx[dim]; - } - dv *= cell_depth; -#elif (AMREX_SPACEDIM == 3) - dv = 1.0; - for (int dim=0; dim<3; ++dim){ - dv *= dx[dim]; - } -#endif - - if ((use_multiphase == 1) && (nspecies == 2)){ - value1 = Rho1[0]; - value2 = Rho2[0]; - if ((value1 <= 0.0) || (value2 <= 0.0)){ - RhoAv[0] = 0.0; - } else { - RhoAv[0] = std::min( 0.5*(value1 + value2), rho0); - } - RhoAv[1] = rho0 - RhoAv[0]; - - } else { - - for (int comp=1; comp<=nspecies; ++comp){ - value1 = Rho1[comp-1]/molmass_in[comp-1]; //Convert to number density - value2 = Rho2[comp-1]/molmass_in[comp-1]; - - switch (avg_type) { - - case 1 : { // Arithmetic with a C0-smoothed Heaviside - - if ( (value1 <= 0.0) || (value2 <= 0.0) ) { - RhoAv[comp-1] = 0.0; - } else { - tmp1 = std::min(dv*value1,1.0); - tmp2 = std::min(dv*value2,1.0); - RhoAv[comp-1] = molmass_in[comp-1]*(value1+value2)/2.0*tmp1*tmp2; - } - - } break; - - case 2 : // Geometric - RhoAv[comp-1] = molmass_in[comp-1]*std::sqrt(std::max(value1,0.0)*std::max(value2,0.0)); - break; - - case 3 : { // Harmonic - // What we want here is the harmonic mean of max(value1,0) and max(value2,0) - // Where we define the result to be zero if either one is zero - // But numerically we want to avoid here division by zero - Real tiny = std::numeric_limits::min(); //replacement for fortran tiny(1.d0) function call - if ( (value1 <= 10.0*tiny) || (value2 <= 10.0*tiny) ) { - RhoAv[comp-1] = 0.0; - } else { - RhoAv[comp-1] = molmass_in[comp-1]*2.0/(1.0/value1+1.0/value2); - } - } break; - - case 10 : { // Arithmetic with (discontinuous) Heaviside - if ( (value1 <= 0.0) || (value2 <= 0.0) ) { - RhoAv[comp-1] = 0.0; - } else { - RhoAv[comp-1] = molmass_in[comp-1]*(value1+value2)/2.0; - } - } break; - - case 11 : { // Arithmetic with C1-smoothed Heaviside - if ( (value1 <= 0.0) || (value2 <= 0.0) ) { - RhoAv[comp-1] = 0.0; - } else { - tmp1 = dv*value1; - if (tmp1 < 1.0) { - tmp1 = (3.0-2.0*tmp1)*tmp1*tmp1; - } else { - tmp1 = 1.0; - } - tmp2 = dv*value2; - if (tmp2 < 1.0) { - tmp2 = (3.0 - 2.0*tmp2)*tmp2*tmp2; - } else { - tmp2 = 1.0; - } - RhoAv[comp-1] = molmass_in[comp-1]*(value1 + value2)/2.0*tmp1*tmp2; - } - } break; - - case 12 : { // Arithmetic with C2-smoothed Heaviside - if ( (value1 <= 0.0) || (value2 <= 0.0) ) { - RhoAv[comp-1] = 0.0; - } else { - tmp1 = dv*value1; - if (tmp1 < 1.0) { - tmp1=(10.0 - 15.0*tmp1 + 6.0*tmp1*tmp1)*std::pow(tmp1,3); - } else { - tmp1=1.0; - } - tmp2 = dv*value2; - if (tmp2 < 1.0) { - tmp2 = (10.0-15.0*tmp2+6.0*tmp2*tmp2)*std::pow(tmp2,3); - } else { - tmp2 = 1.0; - } - RhoAv[comp-1] = molmass_in[comp-1]*(value1+value2)/2.0*tmp1*tmp2; - } - } break; - - default : - printf("compute_nonnegative_rho_av: invalid avg_type\n"); - Abort(); - } - } - } -} - - - -/** - * \brief Water-glycerol mixtures near room temperature - * - * \param[in,out] chi - * \param[in] RhoN - * \param[in] rhotot - */ - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ChiWaterGlycerol ( - Real& chi, - GpuArray& RhoN, - Real rhotot_in){ - - // local - Real c_loc; - - // mass fraction of glycerol - c_loc = RhoN[0]/rhotot_in; - - // chi = chi0 * rational function - chi = Dbar[0]*(1.024-1.001692692*c_loc)/(1.0+0.6632641981*c_loc); - -} - - -/** - * \brief compute D_bar - * - * \param[in] RhoN - * \param[in] rhotot - * \param[in,out] D_bar - */ - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeDBarLocal ( - GpuArray& RhoN, - Real rhotot_in, - Array2D& D_bar){ - - // off-diagonal components of symmetric matrices - Array1D DBar_loc; - - switch (std::abs(mixture_type)) { - - case 1 : { // water-glycerol - - //Note: This case could not be verified. - printf("mixture_type=1 is not currently supported.\n"); - Abort(); - - // we require nspecies=2 - // Dbar(1) = chi0 in the binary notation - if (nspecies != 2){ - printf("mixture_properties_mass_local assumes nspecies=2 if mixture_type=3 (water-glycerol)\n"); - Abort(); - } - - ChiWaterGlycerol(DBar_loc(1), RhoN, rhotot_in); - - } break; - - case 2 : { // Electrolyte mixture - - //Note: The case could not be verified. - printf("mixture_type=2 is not currently supported.\n"); - Abort(); - - if (nspecies != 3){ - printf("mixture_properties_mass_local assumes nspecies=3 if mixture_type=2 (water-glycerol)\n"); - Abort(); - } - - // This is the leading-order correction for dilute solutions - // In particular, the cross-diffusion coefficient ~sqrt(concentration) as per renormalization theory - // The ordering of the values is D_12; D_13, D_23 - - DBar_loc(1) = Dbar[0]*std::sqrt(RhoN[0]/rhotot_in); // counter-ion cross coefficient D_12~sqrt(w) - DBar_loc(2) = Dbar[1]; // D_13 = self diffusion of first ion - DBar_loc(3) = Dbar[2]; // D_23 = self diffusion of second ion - - } break; - - default : //Keep it constant - - for ( int n=1; n <= nspecies*(nspecies-1)/2; ++n){ - DBar_loc(n) = Dbar[n-1]; - } - - } - - // Complete the process by filling the matrices using generic formulae -- this part should not change - // populate D_bar and Hessian matrix - int n=0; - for (int row=1; row<=nspecies; ++row){ - for (int column=1; column<=row-1; ++column){ - n = n+1; - D_bar(row, column) = DBar_loc(n); // SM-diffcoeff's read from input - D_bar(column, row) = D_bar(row, column); // symmetric - } - // populate diagonals - D_bar(row, row) = 0.0; // as self-diffusion is zero - } -} - -/** - * \brief lower triangle and diagonal are overwritten by cholesky factor - * - * \param[in] A input matrix - * \param[in] np - * - * A is ths input matrix. Upon return the lower triangle and diagonal - * are overwritten by the cholesky factor. - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void CholeskyDC (Array2D& A, - int np){ - - int ising; - Real sum1; - - Array1D P; - - Real small_number = 0.0; // Some tolerance - - for (int i=1; i<=np; ++i){ - ising = 0; - - for (int j=i; j<=np; ++j){ - sum1 = A(i,j); - - for (int k=i-1; k>=1; --k){ - sum1 = sum1 - A(i,k)*A(j,k); - } - if(i==j){ - - if(sum1<=small_number){ - P(i) = 0.0; - ising = 1; - - } else { - P(i) = std::sqrt(sum1); - } - } else { - - if( ising==0 ){ - A(j,i) = sum1/P(i); - } else { - A(j,i) = 0.0; - } - } - } - } - - for (int i = 1; i<=np; ++i){ - for (int j = i+1; j<=np; ++j){ - A(i,j) = 0.0; // Zero upper triangle - } - A(i,i) = P(i); - } - -} - -/** - * \brief Compute Sqrt L-Onsager on subset - * - * \param[in] dest map from full system to subset - * \param[in] Molmass_sub - * \param[in] Rho_sub - * \param[in] rhotot_sub - * \param[in] nspecies_sub - * \param[in,out] SqrtLOnsager_sub - * - * dest is a mapping used to eliminate elements in D_bar we don't need (for D_bar_sub) - * and for expanding sqrtLonsager_sub into sqrtLonsager - * it will contain the numbers 1, 2, ..., (nspecies-ntrace) - * with zeros in elements corresponding to trace elements - * (example) for a 5-species system having trace species 2 and 5: - * species 1 2 3 4 5 - * dest(species) 1 0 2 3 0 - * - */ - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeSqrtLOnsagerSub ( - Array1D& dest, - GpuArray& Molmass_in, - GpuArray& Rho_in, - Array2D& DBar_in, - Real rhotot_in, - int nspecies_sub, - Array2D& SqrtLOnsager){ - - // local variables - GpuArray Molmass_sub; - GpuArray Rho_sub; - GpuArray Molarconc_sub; - - Real molmtot_sub; - Array1D W_sub; - Array2D Chi_sub; - Array2D DBar_sub; - Array2D SqrtLOnsager_sub; - - // create a vector of non-trace densities and molmass for the subsystem - for (int row=1; row<=nspecies; ++row){ - if (dest(row) != 0){ - Molmass_sub[dest(row)-1] = Molmass_in[row-1]; - Rho_sub[dest(row)-1] = Rho_in[row-1]; - } - } - - // renormalize total density and mass fractions - Real rhotot_sub = 0.0; - for (int n=0; n& molmass_in, - GpuArray& RhoN, - Real rhotot_in, - Array2D& SqrtLOnsager){ - - // local variables - Array1D W; - - GpuArray Molarconc_loc; - Real molmtot_loc; - Array2D Chi_loc; - Array2D DBar_loc; - - int ntrace; /**< number of trace species with w_k < fractional_tolerance */ - int nspecies_sub; /**< dim of subsystem = nspecies-ntrace */ - - /* - * this is a mapping used to eliminate elements in D_bar we don't need (for D_bar_sub) - * and for expanding sqrtLonsager_sub into sqrtLonsager - * it will contain the numbers 1, 2, ..., (nspecies-ntrace) - * with zeros in elements corresponding to trace elements - * (example) for a 5-species system having trace species 2 and 5: - * species 1 2 3 4 5 - * dest(species) 1 0 2 3 0 - */ - Array1D dest; - - // compute the number of trace species - // build the mapping for expanding/contracting arrays - ntrace = 0; - for (int row=1; row <= nspecies; ++row){ - W(row) = RhoN[row-1]/rhotot_in; - if (W(row) < fraction_tolerance){ - ntrace = ntrace + 1; - dest(row) = 0; - } else { - dest(row) = row - ntrace; - } - } - - if (ntrace == nspecies-1){ - - // this is all trace species except for 1 (essentially pure solvent); - // set sqrtLonsager to zero - MatrixToZeros(MAX_SPECIES, SqrtLOnsager); - - } else if (ntrace == 0) { // there are no trace species - - // compute molarconc and molmtot - ComputeMolconcMolmtotLocal(nspecies, molmass_in, RhoN, rhotot_in, Molarconc_loc, molmtot_loc); - - // compute D_bar - ComputeDBarLocal(RhoN, rhotot_in, DBar_loc); - - // compute chi - ComputeChi(nspecies, molmass_in, RhoN, rhotot_in, Molarconc_loc, Chi_loc, DBar_loc); - - // compute Onsager matrix L (store in sqrtLonsager) - for (int column=1; column<= nspecies; ++column){ - for (int row=1; row<=nspecies; ++row){ - SqrtLOnsager(row, column) = molmtot_loc*rhotot_in*W(row)*Chi_loc(row,column)*W(column)/k_B; - } - } - - // compute cell-centered Cholesky factor, sqrtLonsager - CholeskyDC(SqrtLOnsager,nspecies); - - } else { - - // if there are trace species, we consider a subsystem - // consisting of non-trace species - - nspecies_sub = nspecies - ntrace; - ComputeSqrtLOnsagerSub( dest, - molmass_in, - RhoN, - DBar_loc, - rhotot_in, - nspecies_sub, - SqrtLOnsager); - - } - -} - -/** - * \brief Compute Zeta by temperature - * - * \param[in] MolarConc molar concentration - * \param[in] DBar MS diff-coeffs - * \param[in] Temp Temperature - * \param[in,out] ZetaByTemp zeta/T - * \param[in] DTerm Thermo diff-coeffs - * - * - */ - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void ComputeZetaByTempLocal ( - GpuArray& MolarConc, - Array2D& DBar, - Real Temp, - GpuArray& ZetaByTemp, - GpuArray& DTherm){ - - // local variables - Real Sum_knoti; - Array2D Lambda; - - // compute zeta_by_Temp for thermodiffusion - if (is_nonisothermal == 1){ - - // compute Lambda_ij matrix; molarconc is - // expressed in terms of molmtot,mi,rhotot etc. - for (int row=1; row<=nspecies; ++row){ - for (int column=1; column<=row-1; ++column){ - Lambda(row, column) = -MolarConc[row-1]*MolarConc[column-1]/DBar(row,column); - Lambda(column, row) = Lambda(row, column); - } - } - for (int row=1; row<=nspecies; ++row){ - Sum_knoti = 0.0; - for (int column=1; column<=nspecies; ++column){ - if (column != row){ - Sum_knoti = Sum_knoti + Lambda(row,column)*(DTherm[row-1] - DTherm[column-1]); - } - ZetaByTemp[row-1] = Sum_knoti/Temp; - } - } - } - -} - -/** - * \brief Mixture Mass Properties - * - * \params[in] Rho - * \params[in] rhotot - * \params[out] DBar - * \params[out] DTherm - * \params[out] Hessian - * - * The default case should be to simply set Dbar, Dtherm and Hessian to constants, and read from the input file. - * - * - */ -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void MixturePropsMassLocal (GpuArray& Rho, - Real rhotot, - Array2D& DBar, - GpuArray& DThermLoc, - Array2D& Hessian){ - - // Local values of transport and thermodynamic coefficients (functions of composition!): - GpuArray HOffdiagLoc; //off-diagonal components of symmetric matrices - GpuArray HDiagLoc; // Diagonal component - - // populate D_bar - ComputeDBarLocal(Rho, rhotot, DBar); - - // For now we only encode constant Hessian matrices since we do not have any thermodynamic - // models coded up (Wilson, NTLR, UNIQUAC, etc.) - - for (int n=0; n multispec::fh_kappa; AMREX_GPU_MANAGED amrex::Array2D multispec::fh_chi; AMREX_GPU_MANAGED amrex::GpuArray multispec::fh_monomers; -AMREX_GPU_MANAGED amrex::Real multispec::fh_tension; -AMREX_GPU_MANAGED amrex::GpuArray multispec::contact_angle_lo; -AMREX_GPU_MANAGED amrex::GpuArray multispec::contact_angle_hi; AMREX_GPU_MANAGED amrex::Real multispec::monomer_mass; AMREX_GPU_MANAGED int multispec::n_gex; AMREX_GPU_MANAGED amrex::GpuArray multispec::c_init_1; @@ -90,7 +86,6 @@ void InitializeMultispecNamespace() { use_multiphase = 0; // for RTIL use_flory_huggins = 0; // for flory huggins use_ice_nucleation = 0; // for ice nucleation simulations - use_log_potential = 0; // use quartic potential by default EnScale = 2.79e09; GradEnCoef = 2.73e-6; PotWellDepr = -2.687e-05*T_init[0]*T_init[0]+0.009943*T_init[0]-0.7128; // For polynomial potential @@ -99,12 +94,6 @@ void InitializeMultispecNamespace() { alpha_gex = 0; // for RTIL n_gex = 1; // for RTIL chi_iterations = 10; // number of iterations used in Dbar2chi_iterative - fh_tension = 0.; - - for (int i=0; i fh_kappa; extern AMREX_GPU_MANAGED amrex::Array2D fh_chi; extern AMREX_GPU_MANAGED amrex::GpuArray fh_monomers; - extern AMREX_GPU_MANAGED amrex::Real fh_tension; - extern AMREX_GPU_MANAGED amrex::GpuArray contact_angle_lo; - extern AMREX_GPU_MANAGED amrex::GpuArray contact_angle_hi; extern AMREX_GPU_MANAGED amrex::Real monomer_mass; extern AMREX_GPU_MANAGED amrex::Real kc_tension; extern AMREX_GPU_MANAGED amrex::Real alpha_gex; @@ -61,9 +58,8 @@ namespace multispec { extern amrex::Real L_trans; extern amrex::Real L_zero; - // ice-nucleation + extern AMREX_GPU_MANAGED int use_ice_nucleation; - extern AMREX_GPU_MANAGED int use_log_potential; extern AMREX_GPU_MANAGED amrex::Real EnScale; extern AMREX_GPU_MANAGED amrex::Real GradEnCoef; extern AMREX_GPU_MANAGED amrex::Real PotWellDepr;