diff --git a/Exec/GNUmakefile b/Exec/GNUmakefile index c8e3091..fa46309 100644 --- a/Exec/GNUmakefile +++ b/Exec/GNUmakefile @@ -1,5 +1,5 @@ # AMREX_HOME defines the directory in which we will find all the AMReX code. -AMREX_HOME ?= ../../amrex +AMREX_HOME = ../../amrex CODE_HOME := .. DEBUG = FALSE diff --git a/Exec/inputs_mfim_eb b/Exec/inputs_mfim_eb new file mode 100644 index 0000000..d5dd15a --- /dev/null +++ b/Exec/inputs_mfim_eb @@ -0,0 +1,135 @@ +################################# +###### PROBLEM DOMAIN ###### +################################# +amrex.the_arena_is_managed = 1 + +domain.prob_lo = -64.e-8 -4.e-8 0.0e-8 +domain.prob_hi = 64.e-8 4.e-8 6.0e-8 + +domain.n_cell = 1024 64 64 + +# domain.prob_lo = -4.e-8 -4.e-8 0.0e-8 +# domain.prob_hi = 4.e-8 4.e-8 6.0e-8 + +# domain.n_cell = 64 64 64 + +domain.max_grid_size = 64 64 64 +domain.blocking_factor = 64 64 64 + +prob_type = 1 + +TimeIntegratorOrder = 1 + +nsteps = 3 +plot_int = 1 +# plot_int = 10 +Remnant_P = 0.0 0.0 0.26 +noise_amplitude = 0.001 + +hardswitch_flag = 1 +nucleation_flag = 1 +hardswitch_ratio = 0.2 +nucleation_ratio = 0.3 +hardswitch_alpha_ratio = 13.0 +plot_mat_alpha = 1 +dt = 10.0e-8 + +############################################ +###### POLARIZATION BOUNDARY CONDITIONS #### +############################################ + +P_BC_flag_lo = 3 3 1 # dP/dz= P_lo/lambda; +P_BC_flag_hi = 3 3 1 # dP/dz= P_hi/lambda; +lambda = 5.e-9 + +############################################ +###### ELECTRICAL BOUNDARY CONDITIONS ###### +############################################ + +domain.is_periodic = 1 1 0 + +boundary.hi = per per dir(0.3) +boundary.lo = per per dir(0.0) + +voltage_sweep = 1 +Phi_Bc_lo = 0.0 +Phi_Bc_hi = 0.3 + +Phi_Bc_inc = -0.03 +Phi_Bc_hi_max = 0.3 +phi_tolerance = 1.e-6 +num_Vapp_max = 41 #121 +theta_dep = 0.80 + +use_Fermi_Dirac = 1 +use_work_function = 0 #1 +################################# +###### STACK GEOMETRY ########### +################################# + +#unused +SC_lo = -1.0 -1.0 -1.0 +SC_hi = -1.0 -1.0 -1.0 + +DE_lo = -64.0e-8 -4.0e-8 0.0 +DE_hi = 64.0e-8 4.0e-8 5.0e-9 + +DE1_lo = -64.0e-8 -4.0e-8 55.0e-9 +DE1_hi = 64.0e-8 4.0e-8 60.0e-9 + +FE_lo = -64.e-8 -4.e-8 5.0e-9 +FE_hi = 64.e-8 4.e-8 55.0e-9 + +#device_geom.device_geom_function(x,y,z) = "0.0" +#nucleation_geom.nucleation_geom_function(x,y,z) = "0.0" +#device_geom.device_geom_function(x,y,z) = "0.0 + 1.0*(z < 5.0e-9)" +#device_geom.device_geom_function(x,y,z) = "0.0 + 1.0*(z < 2.5e-9) + 1.0*(z > 52.5e-9)*(z < 55.0e-9)" +#nucleation_geom.nucleation_geom_function(x,y,z) = "1.0 - 1.0*(z < 2.5e-9) - 1.0*(z > 52.5e-9)*(z < 55.0e-9)" + +#nucleation_geom.nucleation_geom_function(x,y,z) = "0.0" #"-1.*(x > 4.0e-8) * (x < 4.25e-8) * (z > 4.84375e-8) * (z < 5.0e-8) + # +1.*(x > -4.255e-8) * (x < -4.0e-8) * (z > 0.) * (z < 1.5625e-9)" + +################################# +###### MATERIAL PROPERTIES ###### +################################# + +epsilon_0 = 8.85e-12 +# epsilonX_fe = 600.0 +# epsilonZ_fe = 600.0 +# epsilon_de = 300.0 +# epsilon_si = 80.0 +# # -13860999.999999998 -1618875000.0 32016000000.0 0 +# #-7423200.0 -209700000.0 7764000000.0 0 --> Ec=15 kV/cm +# #-18558000.0 -209700000.0 10352000000.0 0 --> Ec=30kV/cm +# alpha = -18558000.0 #-3.7116e7 #-1.386e7 #-7.588e7 # 3.34*(T-381)*1e5*2; T=298; ref: Rabe book +# beta = -209700000.0 #-2.097e8 #-1.618875e9 #-8.388e8 #-2.5902e9 # (4.69*(T-393)-202)*1e6*4 +# gamma = 10352000000.0 #1.0352e10 #3.2e10 #7.764e9 # 4.8024e10 # (-5.52*(298-393) + 276)*1e7*6 +# BigGamma = 1.e-3 +# g11 = 20.e-11 # zz +# g44 = 1.e-11 +# g44_p = 0.0 +# g12 = 0.0 +# alpha_12 = 0.0 +# alpha_112 = 0.0 +# alpha_123 = 0.0 + +material_properties.BigGamma_function(x,y,z) = 1.e-3 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.landau_alpha_function(x,y,z) = -18558000.0 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.landau_beta_function(x,y,z) = -209700000.0 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.landau_gamma_function(x,y,z) = 10352000000.0 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.g11_function(x,y,z) = 20.e-11 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.g44_function(x,y,z) = 1.e-11 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.g44_p_function(x,y,z) = 0.0 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.g12_function(x,y,z) = 0.0 +material_properties.alpha_12_function(x,y,z) = 0.0 +material_properties.alpha_112_function(x,y,z) = 0.0 +material_properties.alpha_123_function(x,y,z) = 0.0 + +material_properties.epsilonX_fe_function(x,y,z) = 600.0 * (z<=55.e-9 ) * ( z >= 5.e-9) + 1.0 * (z>55.e-9 ) + 1.0*(z<5.e-9 ) +material_properties.epsilonZ_fe_function(x,y,z) = 600.0 * (z<=55.e-9 ) * ( z >= 5.e-9) + 1.0 * (z>55.e-9 ) + 1.0*(z<5.e-9 ) +material_properties.epsilon_de_function(x,y,z) = 300.0 * (z>55.e-9) + 300.0 * (z<5.e-9) + 1.0 * (z<=55.e-9 ) * ( z >= 5.e-9) +material_properties.epsilon_si_function(x,y,z) = 80.0 + + +donor_doping = 0.25e28 +acceptor_doping = 0.1e28 \ No newline at end of file diff --git a/Source/FerroX.H b/Source/FerroX.H index 8541df5..6a560c9 100644 --- a/Source/FerroX.H +++ b/Source/FerroX.H @@ -134,6 +134,21 @@ void WritePlotfile(c_FerroX& rFerroX, MultiFab& beta_cc, MultiFab& MaterialMask, MultiFab& tphaseMask, + MultiFab& BigGamma, + MultiFab& alpha, + MultiFab& beta, + MultiFab& gamma, + MultiFab& epsilonX_fe, + MultiFab& epsilonZ_fe, + MultiFab& epsilon_de, + MultiFab& epsilon_si, + MultiFab& g11, + MultiFab& g44, + MultiFab& g44_p, + MultiFab& g12, + MultiFab& alpha_12, + MultiFab& alpha_112, + MultiFab& alpha_123, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 20ef757..8495941 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -1,4 +1,4 @@ -/* Contributors: Prabhat Kumar, Saurabh Sawant +/* Contributors: Prabhat Kumar, Saurabh Sawant, Zhi Jackie Yao * */ #include "FerroX.H" @@ -164,6 +164,19 @@ c_FerroX::InitData () AMREX_GPU_MANAGED int FerroX::nsteps; AMREX_GPU_MANAGED int FerroX::plot_int; + +//hardtoswitch & nucleation mask +AMREX_GPU_MANAGED int FerroX::hardswitch_flag; +AMREX_GPU_MANAGED int FerroX::nucleation_flag; + +//hardtoswitch & nucleation ratio +AMREX_GPU_MANAGED amrex::Real FerroX::hardswitch_ratio; +AMREX_GPU_MANAGED amrex::Real FerroX::nucleation_ratio; + +//hardtoswitch & nucleation alpha +AMREX_GPU_MANAGED amrex::Real FerroX::hardswitch_alpha_ratio; + + // time step AMREX_GPU_MANAGED amrex::Real FerroX::dt; @@ -176,17 +189,34 @@ int FerroX::plot_charge; int FerroX::plot_epsilon; int FerroX::plot_mask; int FerroX::plot_tphase; -int FerroX::plot_alpha; -int FerroX::plot_beta; -int FerroX::plot_theta; +int FerroX::plot_mat_BigGamma; +int FerroX::plot_mat_alpha; +int FerroX::plot_mat_beta; +int FerroX::plot_mat_gamma; +int FerroX::plot_mat_epsilonX_fe; +int FerroX::plot_mat_epsilonZ_fe; +int FerroX::plot_mat_epsilon_de; +int FerroX::plot_mat_epsilon_si; +int FerroX::plot_mat_g11; +int FerroX::plot_mat_g44; +int FerroX::plot_mat_g44_p; +int FerroX::plot_mat_g12; +int FerroX::plot_mat_alpha_12; +int FerroX::plot_mat_alpha_112; +int FerroX::plot_mat_alpha_123; +int FerroX::plot_angle_alpha; +int FerroX::plot_angle_beta; +int FerroX::plot_angle_theta; int FerroX::plot_PhiDiff; // multimaterial stack geometry AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE_lo; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE1_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::FE_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::SC_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE_hi; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE1_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::FE_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::SC_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Channel_hi; @@ -194,22 +224,22 @@ AMREX_GPU_MANAGED amrex::GpuArray FerroX::Channel_l // material parameters AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_0; -AMREX_GPU_MANAGED amrex::Real FerroX::epsilonX_fe; +// AMREX_GPU_MANAGED amrex::Real FerroX::epsilonX_fe; AMREX_GPU_MANAGED amrex::Real FerroX::epsilonX_fe_tphase; -AMREX_GPU_MANAGED amrex::Real FerroX::epsilonZ_fe; -AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_de; -AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_si; -AMREX_GPU_MANAGED amrex::Real FerroX::alpha; // alpha = 2*alpha_1 -AMREX_GPU_MANAGED amrex::Real FerroX::beta; // beta = 4*alpha_11 -AMREX_GPU_MANAGED amrex::Real FerroX::gamma; // gamma = 6*alpha_111 -AMREX_GPU_MANAGED amrex::Real FerroX::BigGamma; -AMREX_GPU_MANAGED amrex::Real FerroX::g11; -AMREX_GPU_MANAGED amrex::Real FerroX::g44; -AMREX_GPU_MANAGED amrex::Real FerroX::g44_p; -AMREX_GPU_MANAGED amrex::Real FerroX::g12; -AMREX_GPU_MANAGED amrex::Real FerroX::alpha_12; -AMREX_GPU_MANAGED amrex::Real FerroX::alpha_112; -AMREX_GPU_MANAGED amrex::Real FerroX::alpha_123; +// AMREX_GPU_MANAGED amrex::Real FerroX::epsilonZ_fe; +// AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_de; +// AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_si; +// AMREX_GPU_MANAGED amrex::Real FerroX::alpha; // alpha = 2*alpha_1 +// AMREX_GPU_MANAGED amrex::Real FerroX::beta; // beta = 4*alpha_11 +// AMREX_GPU_MANAGED amrex::Real FerroX::gamma; // gamma = 6*alpha_111 +// AMREX_GPU_MANAGED amrex::Real FerroX::BigGamma; +// AMREX_GPU_MANAGED amrex::Real FerroX::g11; +// AMREX_GPU_MANAGED amrex::Real FerroX::g44; +// AMREX_GPU_MANAGED amrex::Real FerroX::g44_p; +// AMREX_GPU_MANAGED amrex::Real FerroX::g12; +// AMREX_GPU_MANAGED amrex::Real FerroX::alpha_12; +// AMREX_GPU_MANAGED amrex::Real FerroX::alpha_112; +// AMREX_GPU_MANAGED amrex::Real FerroX::alpha_123; // Constants for SC layer calculations AMREX_GPU_MANAGED amrex::Real FerroX::Nc; @@ -233,6 +263,8 @@ AMREX_GPU_MANAGED amrex::Real FerroX::lambda; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Remnant_P; +AMREX_GPU_MANAGED amrex::Real FerroX::noise_amplitude; +AMREX_GPU_MANAGED amrex::Real FerroX::theta_dep; //problem type : initialization of P for 2D/3D/convergence problems AMREX_GPU_MANAGED int FerroX::prob_type; @@ -293,24 +325,56 @@ void InitializeFerroXNamespace(const amrex::GpuArray DE_lo; + extern AMREX_GPU_MANAGED amrex::GpuArray DE1_lo; extern AMREX_GPU_MANAGED amrex::GpuArray FE_lo; extern AMREX_GPU_MANAGED amrex::GpuArray SC_lo; extern AMREX_GPU_MANAGED amrex::GpuArray DE_hi; + extern AMREX_GPU_MANAGED amrex::GpuArray DE1_hi; extern AMREX_GPU_MANAGED amrex::GpuArray FE_hi; extern AMREX_GPU_MANAGED amrex::GpuArray SC_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Channel_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Channel_lo; - // material parameters + // // material parameters extern AMREX_GPU_MANAGED amrex::Real epsilon_0; - extern AMREX_GPU_MANAGED amrex::Real epsilonX_fe; + // extern AMREX_GPU_MANAGED amrex::Real epsilonX_fe; extern AMREX_GPU_MANAGED amrex::Real epsilonX_fe_tphase; - extern AMREX_GPU_MANAGED amrex::Real epsilonZ_fe; - extern AMREX_GPU_MANAGED amrex::Real epsilon_de; - extern AMREX_GPU_MANAGED amrex::Real epsilon_si; - extern AMREX_GPU_MANAGED amrex::Real alpha; // alpha = 2*alpha_1 - extern AMREX_GPU_MANAGED amrex::Real beta; // beta = 4*alpha_11 - extern AMREX_GPU_MANAGED amrex::Real gamma; // gamma = 6*alpha_111 - extern AMREX_GPU_MANAGED amrex::Real BigGamma; - extern AMREX_GPU_MANAGED amrex::Real g11; - extern AMREX_GPU_MANAGED amrex::Real g44; - extern AMREX_GPU_MANAGED amrex::Real g44_p; - extern AMREX_GPU_MANAGED amrex::Real g12; - extern AMREX_GPU_MANAGED amrex::Real alpha_12; - extern AMREX_GPU_MANAGED amrex::Real alpha_112; - extern AMREX_GPU_MANAGED amrex::Real alpha_123; + // extern AMREX_GPU_MANAGED amrex::Real epsilonZ_fe; + // extern AMREX_GPU_MANAGED amrex::Real epsilon_de; + // extern AMREX_GPU_MANAGED amrex::Real epsilon_si; + // extern AMREX_GPU_MANAGED amrex::Real alpha; // alpha = 2*alpha_1 + // extern AMREX_GPU_MANAGED amrex::Real beta; // beta = 4*alpha_11 + // extern AMREX_GPU_MANAGED amrex::Real gamma; // gamma = 6*alpha_111 + // extern AMREX_GPU_MANAGED amrex::Real BigGamma; + // extern AMREX_GPU_MANAGED amrex::Real g11; + // extern AMREX_GPU_MANAGED amrex::Real g44; + // extern AMREX_GPU_MANAGED amrex::Real g44_p; + // extern AMREX_GPU_MANAGED amrex::Real g12; + // extern AMREX_GPU_MANAGED amrex::Real alpha_12; + // extern AMREX_GPU_MANAGED amrex::Real alpha_112; + // extern AMREX_GPU_MANAGED amrex::Real alpha_123; // Constants for SC layer calculations extern AMREX_GPU_MANAGED amrex::Real Nc; @@ -73,7 +99,8 @@ namespace FerroX { extern AMREX_GPU_MANAGED amrex::GpuArray P_BC_flag_lo; extern AMREX_GPU_MANAGED amrex::GpuArray P_BC_flag_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Remnant_P; - + extern AMREX_GPU_MANAGED amrex::Real noise_amplitude; + extern AMREX_GPU_MANAGED amrex::Real theta_dep; //problem type : initialization of P for 2D/3D/convergence problems extern AMREX_GPU_MANAGED int prob_type; diff --git a/Source/Plotfile.cpp b/Source/Plotfile.cpp index 0887f47..2aada24 100644 --- a/Source/Plotfile.cpp +++ b/Source/Plotfile.cpp @@ -13,6 +13,21 @@ void WritePlotfile(c_FerroX& rFerroX, MultiFab& beta_cc, MultiFab& MaterialMask, MultiFab& tphaseMask, + MultiFab& BigGamma, + MultiFab& alpha, + MultiFab& beta, + MultiFab& gamma, + MultiFab& epsilonX_fe, + MultiFab& epsilonZ_fe, + MultiFab& epsilon_de, + MultiFab& epsilon_si, + MultiFab& g11, + MultiFab& g44, + MultiFab& g44_p, + MultiFab& g12, + MultiFab& alpha_12, + MultiFab& alpha_112, + MultiFab& alpha_123, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, @@ -85,19 +100,80 @@ void WritePlotfile(c_FerroX& rFerroX, var_names.push_back("tphase"); } - if (plot_alpha) { + if (plot_mat_BigGamma) { + ++nvar; + var_names.push_back("BigGamma"); + } + if (plot_mat_alpha) { ++nvar; var_names.push_back("alpha"); } - - if (plot_beta) { + if (plot_mat_beta) { ++nvar; var_names.push_back("beta"); } + if (plot_mat_gamma) { + ++nvar; + var_names.push_back("gamma"); + } + if (plot_mat_epsilonX_fe) { + ++nvar; + var_names.push_back("epsilonX_fe"); + } + if (plot_mat_epsilonZ_fe) { + ++nvar; + var_names.push_back("epsilonZ_fe"); + } + if (plot_mat_epsilon_de) { + ++nvar; + var_names.push_back("epsilon_de"); + } + if (plot_mat_epsilon_si) { + ++nvar; + var_names.push_back("epsilon_si"); + } + if (plot_mat_g11) { + ++nvar; + var_names.push_back("g11"); + } + if (plot_mat_g44) { + ++nvar; + var_names.push_back("g44"); + } + if (plot_mat_g44_p) { + ++nvar; + var_names.push_back("g44_p"); + } + if (plot_mat_g12) { + ++nvar; + var_names.push_back("g12"); + } + if (plot_mat_alpha_12) { + ++nvar; + var_names.push_back("alpha_12"); + } + if (plot_mat_alpha_112) { + ++nvar; + var_names.push_back("alpha_112"); + } + if (plot_mat_alpha_123) { + ++nvar; + var_names.push_back("alpha_123"); + } + + if (plot_angle_alpha) { + ++nvar; + var_names.push_back("angle_alpha"); + } + + if (plot_angle_beta) { + ++nvar; + var_names.push_back("angle_beta"); + } - if (plot_theta) { + if (plot_angle_theta) { ++nvar; - var_names.push_back("theta"); + var_names.push_back("angle_theta"); } if (plot_PhiDiff) { @@ -158,15 +234,61 @@ void WritePlotfile(c_FerroX& rFerroX, MultiFab::Copy(Plt, tphaseMask, 0, counter++, 1, 0); } - if (plot_alpha) { + if (plot_mat_BigGamma) { + MultiFab::Copy(Plt, BigGamma, 0, counter++, 1, 0); + } + if (plot_mat_alpha) { + MultiFab::Copy(Plt, alpha, 0, counter++, 1, 0); + } + if (plot_mat_beta) { + MultiFab::Copy(Plt, beta, 0, counter++, 1, 0); + } + if (plot_mat_gamma) { + MultiFab::Copy(Plt, gamma, 0, counter++, 1, 0); + } + if (plot_mat_epsilonX_fe) { + MultiFab::Copy(Plt, epsilonX_fe, 0, counter++, 1, 0); + } + if (plot_mat_epsilonZ_fe) { + MultiFab::Copy(Plt, epsilonZ_fe, 0, counter++, 1, 0); + } + if (plot_mat_epsilon_de) { + MultiFab::Copy(Plt, epsilon_de, 0, counter++, 1, 0); + } + if (plot_mat_epsilon_si) { + MultiFab::Copy(Plt, epsilon_si, 0, counter++, 1, 0); + } + if (plot_mat_g11) { + MultiFab::Copy(Plt, g11, 0, counter++, 1, 0); + } + if (plot_mat_g44) { + MultiFab::Copy(Plt, g44, 0, counter++, 1, 0); + } + if (plot_mat_g44_p) { + MultiFab::Copy(Plt, g44_p, 0, counter++, 1, 0); + } + if (plot_mat_g12) { + MultiFab::Copy(Plt, g12, 0, counter++, 1, 0); + } + if (plot_mat_alpha_12) { + MultiFab::Copy(Plt, alpha_12, 0, counter++, 1, 0); + } + if (plot_mat_alpha_112) { + MultiFab::Copy(Plt, alpha_12, 0, counter++, 1, 0); + } + if (plot_mat_alpha_123) { + MultiFab::Copy(Plt, alpha_123, 0, counter++, 1, 0); + } + + if (plot_angle_alpha) { MultiFab::Copy(Plt, angle_alpha, 0, counter++, 1, 0); } - if (plot_beta) { + if (plot_angle_beta) { MultiFab::Copy(Plt, angle_beta, 0, counter++, 1, 0); } - if (plot_theta) { + if (plot_angle_theta) { MultiFab::Copy(Plt, angle_theta, 0, counter++, 1, 0); } diff --git a/Source/Solver/DerivativeAlgorithm.H b/Source/Solver/DerivativeAlgorithm.H index 80ee376..4f918d0 100644 --- a/Source/Solver/DerivativeAlgorithm.H +++ b/Source/Solver/DerivativeAlgorithm.H @@ -193,8 +193,9 @@ using namespace FerroX; } else if (P_BC_flag_lo[2] == 1){ - Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); - return (dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]); // dP/dz = P_lo/lambda; + //Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); + Real dF_lo = F(i,j,k)/lambda; + return (dx[2]*dF_lo - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]); // dP/dz = P_lo/lambda; // Real F_lo = (9. * F(i,j,k) - F(i,j,k+1)) / (3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return -(dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]);// dP/dz = P_lo/lambda; @@ -222,8 +223,9 @@ using namespace FerroX; } else if (P_BC_flag_hi[2] == 1){ - Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); - return (dx[2]*F_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; + //Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); + Real dF_hi = -F(i,j,k)/lambda; + return (dx[2]*dF_hi + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; // Real F_hi = (9. * F(i,j,k) - F(i,j,k-1)) / ( - 3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return -(dx[2]*F_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; @@ -407,8 +409,9 @@ using namespace FerroX; } else if (P_BC_flag_lo[2] == 1){ - Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); - return (-dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];//dPdz = P_lo/lambda; + //Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda); + Real dF_lo = F(i,j,k)/lambda; + return (-dx[2]*dF_lo - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];//dPdz = P_lo/lambda; // Real F_lo = (9. * F(i,j,k) - F(i,j,k+1)) / (3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return (-dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];// dPdz = P_lo/lambda; @@ -436,8 +439,9 @@ using namespace FerroX; } else if (P_BC_flag_hi[2] == 1){ - Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); - return (dx[2]*F_hi/lambda - F(i,j,k) + F(i,j,k-1))/dx[2]/dx[2];//dPdz = P_hi/lambda; + //Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); + Real dF_hi = -F(i,j,k)/lambda; + return (dx[2]*dF_hi - F(i,j,k) + F(i,j,k-1))/dx[2]/dx[2];//dPdz = P_hi/lambda; // Real F_hi = (9. * F(i,j,k) - F(i,j,k-1)) / ( - 3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative // return (dx[2]*F_hi/lambda - F(i,j,k) + F(i,j,k-1))/dx[2]/dx[2]; // dPdz = P_hi/lambda; diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index 128e6cb..d581b5c 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -19,7 +19,9 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, Array &P_old, MultiFab& rho, MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, const Geometry& geom); //void ComputeEfromPhi(MultiFab& PoissonPhi, @@ -30,14 +32,19 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, // void ComputeEfromPhi(MultiFab& PoissonPhi, Array &E, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, const Geometry& geom, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); //mask based permittivity void InitializePermittivity(std::array,2>& LinOpBCType_2d, - MultiFab& beta_cc, + MultiFab& beta_cc, + MultiFab& epsilonX_fe, + MultiFab& epsilon_de, + MultiFab& epsilon_si, const MultiFab& MaterialMask, const MultiFab& tphaseMask, const amrex::GpuArray& n_cell, @@ -53,7 +60,9 @@ void dF_dPhi(MultiFab& alpha_cc, MultiFab& e_den, MultiFab& p_den, MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, const Geometry& geom, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); @@ -63,6 +72,7 @@ void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& alpha_cc); void SetPhiBC_z(MultiFab& PossonPhi, const amrex::GpuArray& n_cell, const Geometry& geom); +void SetPhiBC_z(MultiFab& PoissonPhi, Array &P_old, MultiFab& MaterialMask, const amrex::GpuArray& n_cell, const Geometry& geom); void SetPoissonBC(c_FerroX& rFerroX, std::array,2>& LinOpBCType_2d, bool& all_homogeneous_boundaries, bool& some_functionbased_inhomogeneous_boundaries, bool& some_constant_inhomogeneous_boundaries); @@ -70,11 +80,14 @@ void Fill_Constant_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& Poisson void Fill_FunctionBased_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time); void CheckSteadyState(MultiFab& PoissonPhi, MultiFab& PoissonPhi_Old, MultiFab& Phidiff, Real phi_tolerance, int step, int& steady_state_step, int& inc_step); +void SetNucleation(Array &P_old, MultiFab& NucleationMask, const amrex::GpuArray& n_cell, amrex::Real hardswitch_ratio, amrex::Real nucleation_ratio); void SetupMLMG(std::unique_ptr& pMLMG, std::unique_ptr& p_mlabec, std::array,2>& LinOpBCType_2d, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + Array& P_old, + MultiFab& MaterialMask, c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info); #ifdef AMREX_USE_EB @@ -94,14 +107,17 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, MultiFab& PoissonPhi, MultiFab& PoissonPhi_Prev, MultiFab& PhiErr, - Array& P_old, + Array& P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + MultiFab& MaterialMask, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, const Geometry& geom, - const amrex::GpuArray& prob_lo, + const amrex::GpuArray& n_cell, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); #ifdef AMREX_USE_EB @@ -112,13 +128,15 @@ void ComputePhi_Rho_EB(std::unique_ptr& pMLMG, MultiFab& PoissonPhi, MultiFab& PoissonPhi_Prev, MultiFab& PhiErr, - Array& P_old, + Array& P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, - const Geometry& geom, - const amrex::GpuArray& prob_lo, + MultiFab& MaterialMask, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); #endif diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 9286cf5..fb15f0c 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -3,15 +3,29 @@ #include "ChargeDensity.H" #include "Utils/eXstaticUtils/eXstaticUtil.H" #include "Utils/FerroXUtils/FerroXUtil.H" +#include "TotalEnergyDensity.H" - -void ComputePoissonRHS(MultiFab& PoissonRHS, +void ComputePoissonRHS(MultiFab& PoissonRHS, Array &P_old, MultiFab& rho, - MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + MultiFab& MaterialMask, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, const Geometry& geom) { + //Real average_P_r = 0.; + //Real total_P_r = 0.; + //Real FE_index_counter = 0.; + + //Compute_P_av(P_old, total_P_r, MaterialMask, FE_index_counter, average_P_r); + + //Real FE_thickness = FE_hi[2] - FE_lo[2]; + //Real one_m_theta = 1.0 - theta_dep; + //Real E_dep = average_P_r/(epsilonZ_fe*epsilon_0)*one_m_theta; + + //amrex::Print() << " E_dep = " << E_dep << "\n"; + //amrex::Print() << " FE_thickness = " << FE_thickness << "\n"; for ( MFIter mfi(PoissonRHS); mfi.isValid(); ++mfi ) { const Box& bx = mfi.validbox(); @@ -73,7 +87,7 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, } else { //mask(i,j,k) == 0.0 FE region RHS(i,j,k) = - (R_11*DPDx(pOld_p, mask, i, j, k, dx) + R_12*DPDy(pOld_p, mask, i, j, k, dx) + R_13*DPDz(pOld_p, mask, i, j, k, dx)) - (R_21*DPDx(pOld_q, mask, i, j, k, dx) + R_22*DPDy(pOld_q, mask, i, j, k, dx) + R_23*DPDz(pOld_q, mask, i, j, k, dx)) - - (R_31*DPDx(pOld_r, mask, i, j, k, dx) + R_32*DPDy(pOld_r, mask, i, j, k, dx) + R_33*DPDz(pOld_r, mask, i, j, k, dx)); + - (R_31*DPDx(pOld_r, mask, i, j, k, dx) + R_32*DPDy(pOld_r, mask, i, j, k, dx) + R_33*DPDz(pOld_r, mask, i, j, k, dx));// - E_dep*FE_thickness; } @@ -85,14 +99,16 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, void dF_dPhi(MultiFab& alpha_cc, MultiFab& PoissonRHS, MultiFab& PoissonPhi, - Array& P_old, + Array& P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, - const Geometry& geom, - const amrex::GpuArray& prob_lo, + MultiFab& MaterialMask, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { @@ -133,10 +149,12 @@ void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, void ComputeEfromPhi(MultiFab& PoissonPhi, Array& E, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, - const Geometry& geom, - const amrex::GpuArray& prob_lo, - const amrex::GpuArray& prob_hi) + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi) { // Calculate E from Phi @@ -204,12 +222,15 @@ void ComputeEfromPhi(MultiFab& PoissonPhi, void InitializePermittivity(std::array,2>& LinOpBCType_2d, MultiFab& beta_cc, - const MultiFab& MaterialMask, - const MultiFab& tphaseMask, - const amrex::GpuArray& n_cell, - const Geometry& geom, + MultiFab& epsilonX_fe, + MultiFab& epsilon_de, + MultiFab& epsilon_si, + const MultiFab& MaterialMask, + const MultiFab& tphaseMask, + const amrex::GpuArray& n_cell, + const Geometry& geom, const amrex::GpuArray& prob_lo, - const amrex::GpuArray& prob_hi) + const amrex::GpuArray& prob_hi) { beta_cc.setVal(0.0); @@ -222,6 +243,9 @@ void InitializePermittivity(std::array& beta = beta_cc.array(mfi); + const Array4& mat_epsilonX_fe = epsilonX_fe.array(mfi); + const Array4& mat_epsilon_de = epsilon_de.array(mfi); + const Array4& mat_epsilon_si = epsilon_si.array(mfi); const Array4& mask = MaterialMask.array(mfi); const Array4& tphase = tphaseMask.array(mfi); @@ -236,18 +260,21 @@ void InitializePermittivity(std::array= t_phase_lo[0] && y <= t_phase_hi[1] && y >= t_phase_lo[1] && z <= t_phase_hi[2] && z >= t_phase_lo[2]){ if(tphase(i,j,k) == 1.0){ beta(i,j,k) = epsilonX_fe_tphase * epsilon_0; } } else if(mask(i,j,k) == 1.0) { - beta(i,j,k) = epsilon_de * epsilon_0; //DE layer + beta(i,j,k) = mat_epsilon_de(i,j,k) * epsilon_0; //DE layer + // printf("epsilon in DE is %g \n", beta(i,j,k)); } else if (mask(i,j,k) >= 2.0){ - beta(i,j,k) = epsilon_si * epsilon_0; //SC layer + beta(i,j,k) = mat_epsilon_si(i,j,k) * epsilon_0; //SC layer } else { - beta(i,j,k) = epsilon_de * epsilon_0; //Spacer is same as DE + beta(i,j,k) = mat_epsilon_de(i,j,k) * epsilon_0; //Spacer is same as DE } }); @@ -499,8 +526,47 @@ void Fill_FunctionBased_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& Po } } +void SetPhiBC_z(MultiFab& PoissonPhi, Array &P_old, MultiFab& MaterialMask, const amrex::GpuArray& n_cell, const Geometry& geom) +{ +// +// Real average_P_r = 0.; +// Real total_P_r = 0.; +// Real FE_index_counter = 0.; +// +// Compute_P_av(P_old, total_P_r, MaterialMask, FE_index_counter, average_P_r); +// +// Real FE_thickness = FE_hi[2] - FE_lo[2]; +// Real one_m_theta = 1.0 - theta_dep; +// Real E_dep = -1.0*average_P_r/(epsilonZ_fe*epsilon_0)*one_m_theta; +// +// amrex::Print() << "average_P_r = " << average_P_r << ", "<< ", E_dep = " << E_dep << ", theta = " << theta_dep << ", phi_Bc_hi" << Phi_Bc_hi << "\n"; +// + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.growntilebox(1); + + const Array4& Phi = PoissonPhi.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + if(k < 0) { + Phi(i,j,k) = Phi_Bc_lo; + } else if(k >= n_cell[2]){ + amrex::Real Eg = bandgap; + amrex::Real Chi = affinity; + amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc/Nv)/q; + amrex::Real phi_m = use_work_function ? metal_work_function : phi_ref; //in eV When not used, applied voltgae is set as the potential on the metal interface + Phi(i,j,k) = (Phi_Bc_hi - (phi_m - phi_ref))*theta_dep;// - E_dep*FE_thickness; //multiplying by theta_dep + //if(i == 5 && j == 5) amrex::Print() << "Phi(i,j,k) = "<< Phi(i,j,k) << std::endl; + } + }); + } + PoissonPhi.FillBoundary(geom.periodicity()); +} + void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray& n_cell, const Geometry& geom) { + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(1); @@ -516,13 +582,86 @@ void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray amrex::Real Chi = affinity; amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc/Nv)/q; amrex::Real phi_m = use_work_function ? metal_work_function : phi_ref; //in eV When not used, applied voltgae is set as the potential on the metal interface - Phi(i,j,k) = Phi_Bc_hi - (phi_m - phi_ref); + Phi(i,j,k) = (Phi_Bc_hi - (phi_m - phi_ref)); + //if(i == 5 && j == 5) amrex::Print() << "Phi(i,j,k) = "<< Phi(i,j,k) << std::endl; } }); } PoissonPhi.FillBoundary(geom.periodicity()); } +void SetNucleation(Array &P_old, MultiFab& NucleationMask, const amrex::GpuArray& n_cell, amrex::Real hardswitch_ratio, amrex::Real nucleation_ratio) +{ + int seed = random_seed; + + //process values (e.g hardswitch_ratio 0.01, nucleation_ratio 0.02) + amrex::Real lower_bd = hardswitch_ratio; // 0.01 + amrex::Real mid_bd = hardswitch_ratio + nucleation_ratio / 2; //0.01 + 0.01 = 0.02 + amrex::Real higher_bd = hardswitch_ratio + nucleation_ratio; //0.01 + 0.02 = 0.03 + + int nprocs = ParallelDescriptor::NProcs(); + + if (prob_type == 1) { + amrex::InitRandom(seed , nprocs, seed ); // give all MPI ranks the same seed + } else { + amrex::InitRandom(seed+ParallelDescriptor::MyProc(), nprocs, seed+ParallelDescriptor::MyProc()); // give all MPI ranks a different seed + } + + int nrand = n_cell[0]*n_cell[2]; + amrex::Gpu::ManagedVector rngs(nrand, 0.0); + + // generate random numbers on the host + for (int i=0; i &pOld_p = P_old[0].array(mfi); + const Array4 &pOld_q = P_old[1].array(mfi); + const Array4 &pOld_r = P_old[2].array(mfi); + const Array4& mask = NucleationMask.array(mfi); + + Real* rng = rngs.data(); + + amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::RandomEngine const& engine) noexcept + { + if (mask(i,j,k) == 0.) { + if (prob_type == 1) { //2D + if (rng[i + k*n_cell[2]] > lower_bd && rng[i + k*n_cell[2]] <= mid_bd){ + pOld_p(i,j,k) = Remnant_P[0]; + pOld_q(i,j,k) = Remnant_P[1]; + pOld_r(i,j,k) = Remnant_P[2]; + } else if (rng[i + k*n_cell[2]] > mid_bd && rng[i + k*n_cell[2]] <= higher_bd){ + pOld_p(i,j,k) = -Remnant_P[0]; + pOld_q(i,j,k) = -Remnant_P[1]; + pOld_r(i,j,k) = -Remnant_P[2]; + } else { + pOld_p(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[0]*noise_amplitude; + pOld_q(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[1]*noise_amplitude; + pOld_r(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[2]*noise_amplitude; + } + } else if (prob_type == 2) { //3D + Real rand = Random(engine); + if (rand <= 0.04) { + pOld_p(i,j,k) = Remnant_P[0]; + pOld_q(i,j,k) = Remnant_P[1]; + pOld_r(i,j,k) = Remnant_P[2]; + } else { + pOld_p(i,j,k) += (-1.0 + 2.0*rand)*Remnant_P[0]*noise_amplitude; + pOld_q(i,j,k) += (-1.0 + 2.0*rand)*Remnant_P[1]*noise_amplitude; + pOld_r(i,j,k) += (-1.0 + 2.0*rand)*Remnant_P[2]*noise_amplitude; + } + } + } + }); + } + +} + void CheckSteadyState(MultiFab& PoissonPhi, MultiFab& PoissonPhi_Old, MultiFab& Phidiff, Real phi_tolerance, int step, int& steady_state_step, int& inc_step) { @@ -559,12 +698,17 @@ void CheckSteadyState(MultiFab& PoissonPhi, MultiFab& PoissonPhi_Old, MultiFab& } -void SetupMLMG(std::unique_ptr& pMLMG, - std::unique_ptr& p_mlabec, +void SetupMLMG(std::unique_ptr& pMLMG, + std::unique_ptr& p_mlabec, std::array,2>& LinOpBCType_2d, - const amrex::GpuArray& n_cell, - std::array< MultiFab, AMREX_SPACEDIM >& beta_face, - c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info) + const amrex::GpuArray& n_cell, + std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + Array& P_old, + MultiFab& MaterialMask, + c_FerroX& rFerroX, + MultiFab& PoissonPhi, + amrex::Real& time, + amrex::LPInfo& info) { auto& rGprop = rFerroX.get_GeometryProperties(); auto& geom = rGprop.geom; @@ -597,7 +741,8 @@ void SetupMLMG(std::unique_ptr& pMLMG, PoissonPhi.FillBoundary(geom.periodicity()); // set Dirichlet BC by reading in the ghost cell values - SetPhiBC_z(PoissonPhi, n_cell, geom); + SetPhiBC_z(PoissonPhi, n_cell, geom); + //SetPhiBC_z(PoissonPhi, P_old, MaterialMask, n_cell, geom); p_mlabec->setLevelBC(amrlev, &PoissonPhi); // (A*alpha_cc - B * div beta grad) phi = rhs @@ -617,7 +762,10 @@ void SetupMLMG(std::unique_ptr& pMLMG, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, MultiFab& beta_cc, - c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info) + c_FerroX& rFerroX, + MultiFab& PoissonPhi, + amrex::Real& time, + amrex::LPInfo& info) { auto& rGprop = rFerroX.get_GeometryProperties(); auto& geom = rGprop.geom; @@ -687,8 +835,11 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, MultiFab& e_den, MultiFab& p_den, MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, - const Geometry& geom, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, + const Geometry& geom, + const amrex::GpuArray& n_cell, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) @@ -714,6 +865,8 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, //Initial guess for phi PoissonPhi.setVal(0.); + SetPhiBC_z(PoissonPhi, n_cell, geom); + //SetPhiBC_z(PoissonPhi, P_old, MaterialMask, n_cell, geom); //Poisson Solve pMLMG->solve({&PoissonPhi}, {&PoissonRHS}, 1.e-10, -1); @@ -759,8 +912,10 @@ void ComputePhi_Rho_EB(std::unique_ptr& pMLMG, MultiFab& e_den, MultiFab& p_den, MultiFab& MaterialMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, - const Geometry& geom, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, + const Geometry& geom, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) diff --git a/Source/Solver/Initialization.H b/Source/Solver/Initialization.H index 8a929d5..7b00cfa 100644 --- a/Source/Solver/Initialization.H +++ b/Source/Solver/Initialization.H @@ -8,15 +8,15 @@ using namespace amrex; using namespace FerroX; void InitializePandRho(Array &P_old, - MultiFab& Gamma, + MultiFab& BigGamma, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - const MultiFab& MaterialMask, - const MultiFab& tphaseMask, + const MultiFab& MaterialMask, + const MultiFab& tphaseMask, const amrex::GpuArray& n_cell, const Geometry& geom, - const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); void InitializeMaterialMask(MultiFab& MaterialMask, @@ -28,3 +28,25 @@ void InitializeMaterialMask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& M void Initialize_tphase_Mask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& tphaseMask); void Initialize_Euler_angles(c_FerroX& rFerroX, const Geometry& geom, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta); +void Initialize_MaterialProperties(c_FerroX& rFerroX, const Geometry& geom, + MultiFab& BigGamma, + MultiFab& alpha, + MultiFab& beta, + MultiFab& gamma, + MultiFab& epsilonX_fe, + MultiFab& epsilonZ_fe, + MultiFab& epsilon_de, + MultiFab& epsilon_si, + MultiFab& g11, + MultiFab& g44, + MultiFab& g44_p, + MultiFab& g12, + MultiFab& alpha_12, + MultiFab& alpha_112, + MultiFab& alpha_123); + +void SetHardToSwitchNucleation(MultiFab& alpha, + MultiFab& NucleationMask, + const amrex::GpuArray& n_cell, + amrex::Real hardswitch_ratio, + amrex::Real hardswitch_alpha_ratio); \ No newline at end of file diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 447adf4..6946f1e 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -4,15 +4,15 @@ // INITIALIZE rho in SC region void InitializePandRho(Array &P_old, - MultiFab& Gamma, + MultiFab& BigGamma, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - const MultiFab& MaterialMask, - const MultiFab& tphaseMask, + const MultiFab& MaterialMask, + const MultiFab& tphaseMask, const amrex::GpuArray& n_cell, const Geometry& geom, - const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { @@ -73,7 +73,7 @@ void InitializePandRho(Array &P_old, const Array4 &pOld_p = P_old[0].array(mfi); const Array4 &pOld_q = P_old[1].array(mfi); const Array4 &pOld_r = P_old[2].array(mfi); - const Array4& Gam = Gamma.array(mfi); + const Array4& mat_BigGamma = BigGamma.array(mfi); const Array4& mask = MaterialMask.array(mfi); const Array4& tphase = tphaseMask.array(mfi); @@ -110,8 +110,6 @@ void InitializePandRho(Array &P_old, } - Gam(i,j,k) = BigGamma; - //set t_phase Pz to zero //if(x <= t_phase_hi[0] && x >= t_phase_lo[0] && y <= t_phase_hi[1] && y >= t_phase_lo[1] && z <= t_phase_hi[2] && z >= t_phase_lo[2]){ if(tphase(i,j,k) == 1.0){ @@ -122,7 +120,7 @@ void InitializePandRho(Array &P_old, pOld_p(i,j,k) = 0.0; pOld_q(i,j,k) = 0.0; pOld_r(i,j,k) = 0.0; - Gam(i,j,k) = 0.0; + mat_BigGamma(i,j,k) = 0.0; // Note this is overwriting the initialized Gamma, therefore this function must be called after InitializeMaterialProperties } if (is_polarization_scalar == 1){ @@ -167,8 +165,8 @@ void InitializePandRho(Array &P_old, // create a mask filled with integers to idetify different material types void InitializeMaterialMask(MultiFab& MaterialMask, - const Geometry& geom, - const amrex::GpuArray& prob_lo, + const Geometry& geom, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { // loop over boxes @@ -192,6 +190,8 @@ void InitializeMaterialMask(MultiFab& MaterialMask, mask(i,j,k) = 0.; } else if (x <= DE_hi[0] && x >= DE_lo[0] && y <= DE_hi[1] && y >= DE_lo[1] && z <= DE_hi[2] && z >= DE_lo[2]) { mask(i,j,k) = 1.; + } else if (x <= DE1_hi[0] && x >= DE1_lo[0] && y <= DE1_hi[1] && y >= DE1_lo[1] && z <= DE1_hi[2] && z >= DE1_lo[2]) { + mask(i,j,k) = 1.; } else if (x <= SC_hi[0] && x >= SC_lo[0] && y <= SC_hi[1] && y >= SC_lo[1] && z <= SC_hi[2] && z >= SC_lo[2]) { mask(i,j,k) = 2.; if (x <= Channel_hi[0] && x >= Channel_lo[0] && y <= Channel_hi[1] && y >= Channel_lo[1] && z <= Channel_hi[2] && z >= Channel_lo[2]){ @@ -218,32 +218,33 @@ void InitializeMaterialMask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& M for (MFIter mfi(MaterialMask, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const auto& mask_arr = MaterialMask.array(mfi); - const auto& bx = mfi.tilebox(); + // const auto& bx = mfi.tilebox(); + const Box& bx = mfi.growntilebox(MaterialMask.nGrow()); - std::string m_mask_s; - std::unique_ptr m_mask_parser; + std::string m_mask_s; + std::unique_ptr m_mask_parser; std::string m_str_device_geom_function; - ParmParse pp_mask("device_geom"); + ParmParse pp_mask("device_geom"); - if (pp_mask.query("device_geom_function(x,y,z)", m_str_device_geom_function) ) { - m_mask_s = "parse_device_geom_function"; - } + if (pp_mask.query("device_geom_function(x,y,z)", m_str_device_geom_function) ) { + m_mask_s = "parse_device_geom_function"; + } - if (m_mask_s == "parse_device_geom_function") { - Store_parserString(pp_mask, "device_geom_function(x,y,z)", m_str_device_geom_function); - m_mask_parser = std::make_unique( - makeParser(m_str_device_geom_function,{"x","y","z"})); - } + if (m_mask_s == "parse_device_geom_function") { + Store_parserString(pp_mask, "device_geom_function(x,y,z)", m_str_device_geom_function); + m_mask_parser = std::make_unique( + makeParser(m_str_device_geom_function,{"x","y","z"})); + } - const auto& macro_parser = m_mask_parser->compile<3>(); + const auto& macro_parser = m_mask_parser->compile<3>(); - amrex::ParallelFor(bx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv,macro_parser,mask_arr); - }); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv,macro_parser,mask_arr); + }); } MaterialMask.FillBoundary(geom.periodicity()); @@ -302,78 +303,73 @@ void Initialize_Euler_angles(c_FerroX& rFerroX, const Geometry& geom, MultiFab& const auto dx = rGprop.geom.CellSizeArray(); const auto& real_box = rGprop.geom.ProbDomain(); - const auto iv_alpha = angle_alpha.ixType().toIntVect(); - const auto iv_beta = angle_beta.ixType().toIntVect(); - const auto iv_theta = angle_theta.ixType().toIntVect(); + const auto iv_angle_alpha = angle_alpha.ixType().toIntVect(); + const auto iv_angle_beta = angle_beta.ixType().toIntVect(); + const auto iv_angle_theta = angle_theta.ixType().toIntVect(); for (MFIter mfi(angle_alpha, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const auto& alpha_arr = angle_alpha.array(mfi); - const auto& beta_arr = angle_beta.array(mfi); - const auto& theta_arr = angle_theta.array(mfi); + const auto& angle_alpha_arr = angle_alpha.array(mfi); + const auto& angle_beta_arr = angle_beta.array(mfi); + const auto& angle_theta_arr = angle_theta.array(mfi); const auto& bx = mfi.tilebox(); - std::string alpha_s; - std::unique_ptr alpha_parser; - std::string m_str_alpha_function; + std::string angle_alpha_s; + std::unique_ptr angle_alpha_parser; + std::string m_str_angle_alpha_function; - std::string beta_s; - std::unique_ptr beta_parser; - std::string m_str_beta_function; + std::string angle_beta_s; + std::unique_ptr angle_beta_parser; + std::string m_str_angle_beta_function; - std::string theta_s; - std::unique_ptr theta_parser; - std::string m_str_theta_function; + std::string angle_theta_s; + std::unique_ptr angle_theta_parser; + std::string m_str_angle_theta_function; - ParmParse pp_alpha("angle_alpha"); + ParmParse pp_mask("o_phase_angle"); - if (pp_alpha.query("alpha_function(x,y,z)", m_str_alpha_function) ) { - alpha_s = "parse_alpha_function"; + if (pp_mask.query("angle_alpha_function(x,y,z)", m_str_angle_alpha_function) ) { + angle_alpha_s = "parse_angle_alpha_function"; } - if (alpha_s == "parse_alpha_function") { - Store_parserString(pp_alpha, "alpha_function(x,y,z)", m_str_alpha_function); - alpha_parser = std::make_unique( - makeParser(m_str_alpha_function,{"x","y","z"})); + if (angle_alpha_s == "parse_angle_alpha_function") { + Store_parserString(pp_mask, "angle_alpha_function(x,y,z)", m_str_angle_alpha_function); + angle_alpha_parser = std::make_unique( + makeParser(m_str_angle_alpha_function,{"x","y","z"})); } - ParmParse pp_beta("angle_beta"); - - if (pp_beta.query("beta_function(x,y,z)", m_str_beta_function) ) { - beta_s = "parse_beta_function"; + if (pp_mask.query("angle_beta_function(x,y,z)", m_str_angle_beta_function) ) { + angle_beta_s = "parse_angle_beta_function"; } - if (beta_s == "parse_beta_function") { - Store_parserString(pp_beta, "beta_function(x,y,z)", m_str_beta_function); - beta_parser = std::make_unique( - makeParser(m_str_beta_function,{"x","y","z"})); + if (angle_beta_s == "parse_angle_beta_function") { + Store_parserString(pp_mask, "angle_beta_function(x,y,z)", m_str_angle_beta_function); + angle_beta_parser = std::make_unique( + makeParser(m_str_angle_beta_function,{"x","y","z"})); } - ParmParse pp_theta("angle_theta"); - - - if (pp_theta.query("theta_function(x,y,z)", m_str_theta_function) ) { - theta_s = "parse_theta_function"; + if (pp_mask.query("angle_theta_function(x,y,z)", m_str_angle_theta_function) ) { + angle_theta_s = "parse_angle_theta_function"; } - if (theta_s == "parse_theta_function") { - Store_parserString(pp_theta, "theta_function(x,y,z)", m_str_theta_function); - theta_parser = std::make_unique( - makeParser(m_str_theta_function,{"x","y","z"})); + if (angle_theta_s == "parse_angle_theta_function") { + Store_parserString(pp_mask, "angle_theta_function(x,y,z)", m_str_angle_theta_function); + angle_theta_parser = std::make_unique( + makeParser(m_str_angle_theta_function,{"x","y","z"})); } - const auto& macro_parser_alpha = alpha_parser->compile<3>(); - const auto& macro_parser_beta = beta_parser->compile<3>(); - const auto& macro_parser_theta = theta_parser->compile<3>(); + const auto& macro_parser_angle_alpha = angle_alpha_parser->compile<3>(); + const auto& macro_parser_angle_beta = angle_beta_parser->compile<3>(); + const auto& macro_parser_angle_theta = angle_theta_parser->compile<3>(); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_alpha,macro_parser_alpha,alpha_arr); - eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_beta, macro_parser_beta, beta_arr ); - eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_theta,macro_parser_theta,theta_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_angle_alpha,macro_parser_angle_alpha,angle_alpha_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_angle_beta, macro_parser_angle_beta, angle_beta_arr ); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_angle_theta,macro_parser_angle_theta,angle_theta_arr); }); } @@ -382,3 +378,364 @@ void Initialize_Euler_angles(c_FerroX& rFerroX, const Geometry& geom, MultiFab& angle_theta.FillBoundary(geom.periodicity()); } +// initialization of Material Properties +void Initialize_MaterialProperties(c_FerroX& rFerroX, const Geometry& geom, + MultiFab& BigGamma, + MultiFab& alpha, + MultiFab& beta, + MultiFab& gamma, + MultiFab& epsilonX_fe, + MultiFab& epsilonZ_fe, + MultiFab& epsilon_de, + MultiFab& epsilon_si, + MultiFab& g11, + MultiFab& g44, + MultiFab& g44_p, + MultiFab& g12, + MultiFab& alpha_12, + MultiFab& alpha_112, + MultiFab& alpha_123) +{ + auto& rGprop = rFerroX.get_GeometryProperties(); + Box const& domain = rGprop.geom.Domain(); + + const auto dx = rGprop.geom.CellSizeArray(); + const auto& real_box = rGprop.geom.ProbDomain(); + const auto iv_mat_BigGamma = BigGamma.ixType().toIntVect(); + const auto iv_mat_alpha = alpha.ixType().toIntVect(); + const auto iv_mat_beta = beta.ixType().toIntVect(); + const auto iv_mat_gamma = gamma.ixType().toIntVect(); + const auto iv_mat_epsilonX_fe = epsilonX_fe.ixType().toIntVect(); + const auto iv_mat_epsilonZ_fe = epsilonZ_fe.ixType().toIntVect(); + const auto iv_mat_epsilon_de = epsilon_de.ixType().toIntVect(); + const auto iv_mat_epsilon_si = epsilon_si.ixType().toIntVect(); + const auto iv_mat_g11 = g11.ixType().toIntVect(); + const auto iv_mat_g44 = g44.ixType().toIntVect(); + const auto iv_mat_g44_p = g44_p.ixType().toIntVect(); + const auto iv_mat_g12 = g12.ixType().toIntVect(); + const auto iv_mat_alpha_12 = alpha_12.ixType().toIntVect(); + const auto iv_mat_alpha_112 = alpha_112.ixType().toIntVect(); + const auto iv_mat_alpha_123 = alpha_123.ixType().toIntVect(); + + for (MFIter mfi(alpha, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto& mat_BigGamma_arr = BigGamma.array(mfi); + const auto& mat_alpha_arr = alpha.array(mfi); + const auto& mat_beta_arr = beta.array(mfi); + const auto& mat_gamma_arr = gamma.array(mfi); + const auto& mat_epsilonX_fe_arr = epsilonX_fe.array(mfi); + const auto& mat_epsilonZ_fe_arr = epsilonZ_fe.array(mfi); + const auto& mat_epsilon_de_arr = epsilon_de.array(mfi); + const auto& mat_epsilon_si_arr = epsilon_si.array(mfi); + const auto& mat_g11_arr = g11.array(mfi); + const auto& mat_g44_arr = g44.array(mfi); + const auto& mat_g44_p_arr = g44_p.array(mfi); + const auto& mat_g12_arr = g12.array(mfi); + const auto& mat_alpha_12_arr = alpha_12.array(mfi); + const auto& mat_alpha_112_arr = alpha_112.array(mfi); + const auto& mat_alpha_123_arr = alpha_123.array(mfi); + const auto& bx = mfi.tilebox(); + + std::string mat_BigGamma_s; + std::unique_ptr mat_BigGamma_parser; + std::string m_str_mat_BigGamma_function; + + std::string mat_alpha_s; + std::unique_ptr mat_alpha_parser; + std::string m_str_mat_alpha_function; + + std::string mat_beta_s; + std::unique_ptr mat_beta_parser; + std::string m_str_mat_beta_function; + + std::string mat_gamma_s; + std::unique_ptr mat_gamma_parser; + std::string m_str_mat_gamma_function; + + std::string mat_epsilonX_fe_s; + std::unique_ptr mat_epsilonX_fe_parser; + std::string m_str_mat_epsilonX_fe_function; + + std::string mat_epsilonZ_fe_s; + std::unique_ptr mat_epsilonZ_fe_parser; + std::string m_str_mat_epsilonZ_fe_function; + + std::string mat_epsilon_de_s; + std::unique_ptr mat_epsilon_de_parser; + std::string m_str_mat_epsilon_de_function; + + std::string mat_epsilon_si_s; + std::unique_ptr mat_epsilon_si_parser; + std::string m_str_mat_epsilon_si_function; + + std::string mat_g11_s; + std::unique_ptr mat_g11_parser; + std::string m_str_mat_g11_function; + + std::string mat_g44_s; + std::unique_ptr mat_g44_parser; + std::string m_str_mat_g44_function; + + std::string mat_g44_p_s; + std::unique_ptr mat_g44_p_parser; + std::string m_str_mat_g44_p_function; + + std::string mat_g12_s; + std::unique_ptr mat_g12_parser; + std::string m_str_mat_g12_function; + + std::string mat_alpha_12_s; + std::unique_ptr mat_alpha_12_parser; + std::string m_str_mat_alpha_12_function; + + std::string mat_alpha_112_s; + std::unique_ptr mat_alpha_112_parser; + std::string m_str_mat_alpha_112_function; + + std::string mat_alpha_123_s; + std::unique_ptr mat_alpha_123_parser; + std::string m_str_mat_alpha_123_function; + + ParmParse pp_mask("material_properties"); + + if (pp_mask.query("BigGamma_function(x,y,z)", m_str_mat_BigGamma_function) ) { + mat_BigGamma_s = "parse_BigGamma_function"; + } + if (mat_BigGamma_s == "parse_BigGamma_function") { + Store_parserString(pp_mask, "BigGamma_function(x,y,z)", m_str_mat_BigGamma_function); + mat_BigGamma_parser = std::make_unique( + makeParser(m_str_mat_BigGamma_function,{"x","y","z"})); + } + + if (pp_mask.query("landau_alpha_function(x,y,z)", m_str_mat_alpha_function) ) { + mat_alpha_s = "parse_landau_alpha_function"; + } + if (mat_alpha_s == "parse_landau_alpha_function") { + Store_parserString(pp_mask, "landau_alpha_function(x,y,z)", m_str_mat_alpha_function); + mat_alpha_parser = std::make_unique( + makeParser(m_str_mat_alpha_function,{"x","y","z"})); + } + + if (pp_mask.query("landau_beta_function(x,y,z)", m_str_mat_beta_function) ) { + mat_beta_s = "parse_landau_beta_function"; + } + if (mat_beta_s == "parse_landau_beta_function") { + Store_parserString(pp_mask, "landau_beta_function(x,y,z)", m_str_mat_beta_function); + mat_beta_parser = std::make_unique( + makeParser(m_str_mat_beta_function,{"x","y","z"})); + } + + if (pp_mask.query("landau_gamma_function(x,y,z)", m_str_mat_gamma_function) ) { + mat_gamma_s = "parse_landau_gamma_function"; + } + if (mat_gamma_s == "parse_landau_gamma_function") { + Store_parserString(pp_mask, "landau_gamma_function(x,y,z)", m_str_mat_gamma_function); + mat_gamma_parser = std::make_unique( + makeParser(m_str_mat_gamma_function,{"x","y","z"})); + } + + if (pp_mask.query("epsilonX_fe_function(x,y,z)", m_str_mat_epsilonX_fe_function) ) { + mat_epsilonX_fe_s = "parse_epsilonX_fe_function"; + } + if (mat_epsilonX_fe_s == "parse_epsilonX_fe_function") { + Store_parserString(pp_mask, "epsilonX_fe_function(x,y,z)", m_str_mat_epsilonX_fe_function); + mat_epsilonX_fe_parser = std::make_unique( + makeParser(m_str_mat_epsilonX_fe_function,{"x","y","z"})); + } + + if (pp_mask.query("epsilonZ_fe_function(x,y,z)", m_str_mat_epsilonZ_fe_function) ) { + mat_epsilonZ_fe_s = "parse_epsilonZ_fe_function"; + } + if (mat_epsilonZ_fe_s == "parse_epsilonZ_fe_function") { + Store_parserString(pp_mask, "epsilonZ_fe_function(x,y,z)", m_str_mat_epsilonZ_fe_function); + mat_epsilonZ_fe_parser = std::make_unique( + makeParser(m_str_mat_epsilonZ_fe_function,{"x","y","z"})); + } + + if (pp_mask.query("epsilon_de_function(x,y,z)", m_str_mat_epsilon_de_function) ) { + mat_epsilon_de_s = "parse_epsilon_de_function"; + } + if (mat_epsilon_de_s == "parse_epsilon_de_function") { + Store_parserString(pp_mask, "epsilon_de_function(x,y,z)", m_str_mat_epsilon_de_function); + mat_epsilon_de_parser = std::make_unique( + makeParser(m_str_mat_epsilon_de_function,{"x","y","z"})); + } + + if (pp_mask.query("epsilon_si_function(x,y,z)", m_str_mat_epsilon_si_function) ) { + mat_epsilon_si_s = "parse_epsilon_si_function"; + } + if (mat_epsilon_si_s == "parse_epsilon_si_function") { + Store_parserString(pp_mask, "epsilon_si_function(x,y,z)", m_str_mat_epsilon_si_function); + mat_epsilon_si_parser = std::make_unique( + makeParser(m_str_mat_epsilon_si_function,{"x","y","z"})); + } + + if (pp_mask.query("g11_function(x,y,z)", m_str_mat_g11_function) ) { + mat_g11_s = "parse_g11_function"; + } + if (mat_g11_s == "parse_g11_function") { + Store_parserString(pp_mask, "g11_function(x,y,z)", m_str_mat_g11_function); + mat_g11_parser = std::make_unique( + makeParser(m_str_mat_g11_function,{"x","y","z"})); + } + + if (pp_mask.query("g44_function(x,y,z)", m_str_mat_g44_function) ) { + mat_g44_s = "parse_g44_function"; + } + if (mat_g44_s == "parse_g44_function") { + Store_parserString(pp_mask, "g44_function(x,y,z)", m_str_mat_g44_function); + mat_g44_parser = std::make_unique( + makeParser(m_str_mat_g44_function,{"x","y","z"})); + } + + if (pp_mask.query("g44_p_function(x,y,z)", m_str_mat_g44_p_function) ) { + mat_g44_p_s = "parse_g44_p_function"; + } + if (mat_g44_p_s == "parse_g44_p_function") { + Store_parserString(pp_mask, "g44_p_function(x,y,z)", m_str_mat_g44_p_function); + mat_g44_p_parser = std::make_unique( + makeParser(m_str_mat_g44_p_function,{"x","y","z"})); + } + + if (pp_mask.query("g12_function(x,y,z)", m_str_mat_g12_function) ) { + mat_g12_s = "parse_g12_function"; + } + if (mat_g12_s == "parse_g12_function") { + Store_parserString(pp_mask, "g12_function(x,y,z)", m_str_mat_g12_function); + mat_g12_parser = std::make_unique( + makeParser(m_str_mat_g12_function,{"x","y","z"})); + } + + if (pp_mask.query("alpha_12_function(x,y,z)", m_str_mat_alpha_12_function) ) { + mat_alpha_12_s = "parse_alpha_12_function"; + } + if (mat_alpha_12_s == "parse_alpha_12_function") { + Store_parserString(pp_mask, "alpha_12_function(x,y,z)", m_str_mat_alpha_12_function); + mat_alpha_12_parser = std::make_unique( + makeParser(m_str_mat_alpha_12_function,{"x","y","z"})); + } + + if (pp_mask.query("alpha_112_function(x,y,z)", m_str_mat_alpha_112_function) ) { + mat_alpha_112_s = "parse_alpha_112_function"; + } + if (mat_alpha_112_s == "parse_alpha_112_function") { + Store_parserString(pp_mask, "alpha_112_function(x,y,z)", m_str_mat_alpha_112_function); + mat_alpha_112_parser = std::make_unique( + makeParser(m_str_mat_alpha_112_function,{"x","y","z"})); + } + + if (pp_mask.query("alpha_123_function(x,y,z)", m_str_mat_alpha_123_function) ) { + mat_alpha_123_s = "parse_alpha_123_function"; + } + if (mat_alpha_123_s == "parse_alpha_123_function") { + Store_parserString(pp_mask, "alpha_123_function(x,y,z)", m_str_mat_alpha_123_function); + mat_alpha_123_parser = std::make_unique( + makeParser(m_str_mat_alpha_123_function,{"x","y","z"})); + } + + const auto& macro_parser_mat_BigGamma = mat_BigGamma_parser->compile<3>(); + const auto& macro_parser_mat_alpha = mat_alpha_parser->compile<3>(); + const auto& macro_parser_mat_beta = mat_beta_parser->compile<3>(); + const auto& macro_parser_mat_gamma = mat_gamma_parser->compile<3>(); + const auto& macro_parser_mat_epsilonX_fe = mat_epsilonX_fe_parser->compile<3>(); + const auto& macro_parser_mat_epsilonZ_fe = mat_epsilonZ_fe_parser->compile<3>(); + const auto& macro_parser_mat_epsilon_de = mat_epsilon_de_parser->compile<3>(); + const auto& macro_parser_mat_epsilon_si = mat_epsilon_si_parser->compile<3>(); + const auto& macro_parser_mat_g11 = mat_g11_parser->compile<3>(); + const auto& macro_parser_mat_g44 = mat_g44_parser->compile<3>(); + const auto& macro_parser_mat_g44_p = mat_g44_p_parser->compile<3>(); + const auto& macro_parser_mat_g12 = mat_g12_parser->compile<3>(); + const auto& macro_parser_mat_alpha_12 = mat_alpha_12_parser->compile<3>(); + const auto& macro_parser_mat_alpha_112 = mat_alpha_112_parser->compile<3>(); + const auto& macro_parser_mat_alpha_123 = mat_alpha_123_parser->compile<3>(); + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_alpha,macro_parser_mat_alpha,mat_alpha_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_BigGamma,macro_parser_mat_BigGamma,mat_BigGamma_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_beta, macro_parser_mat_beta, mat_beta_arr ); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_gamma,macro_parser_mat_gamma,mat_gamma_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_epsilonX_fe,macro_parser_mat_epsilonX_fe,mat_epsilonX_fe_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_epsilonZ_fe,macro_parser_mat_epsilonZ_fe,mat_epsilonZ_fe_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_epsilon_de,macro_parser_mat_epsilon_de,mat_epsilon_de_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_epsilon_si,macro_parser_mat_epsilon_si,mat_epsilon_si_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_g11,macro_parser_mat_g11,mat_g11_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_g44,macro_parser_mat_g44,mat_g44_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_g44_p,macro_parser_mat_g44_p,mat_g44_p_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_g12,macro_parser_mat_g12,mat_g12_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_alpha_12,macro_parser_mat_alpha_12,mat_alpha_12_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_alpha_112,macro_parser_mat_alpha_112,mat_alpha_112_arr); + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv_mat_alpha_123,macro_parser_mat_alpha_123,mat_alpha_123_arr); + }); + + } + BigGamma.FillBoundary(geom.periodicity()); + alpha.FillBoundary(geom.periodicity()); + beta.FillBoundary(geom.periodicity()); + gamma.FillBoundary(geom.periodicity()); + epsilonX_fe.FillBoundary(geom.periodicity()); + epsilonZ_fe.FillBoundary(geom.periodicity()); + epsilon_de.FillBoundary(geom.periodicity()); + epsilon_si.FillBoundary(geom.periodicity()); + g11.FillBoundary(geom.periodicity()); + g44.FillBoundary(geom.periodicity()); + g44_p.FillBoundary(geom.periodicity()); + g12.FillBoundary(geom.periodicity()); + alpha_12.FillBoundary(geom.periodicity()); + alpha_112.FillBoundary(geom.periodicity()); + alpha_123.FillBoundary(geom.periodicity()); +} + +void SetHardToSwitchNucleation(MultiFab& alpha, MultiFab& NucleationMask, const amrex::GpuArray& n_cell, amrex::Real hardswitch_ratio, amrex::Real hardswitch_alpha_ratio) +{ + int seed = random_seed; + + int nprocs = ParallelDescriptor::NProcs(); + + if (prob_type == 1) { + amrex::InitRandom(seed , nprocs, seed ); // give all MPI ranks the same seed + } else { + amrex::InitRandom(seed+ParallelDescriptor::MyProc(), nprocs, seed+ParallelDescriptor::MyProc()); // give all MPI ranks a different seed + } + + int nrand = n_cell[0]*n_cell[2]; + amrex::Gpu::ManagedVector rngs(nrand, 0.0); + + // generate random numbers on the host + for (int i=0; i &mat_alpha_arr = alpha.array(mfi); + const Array4& mask = NucleationMask.array(mfi); + + Real* rng = rngs.data(); + + amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::RandomEngine const& engine) noexcept + { + if (mask(i,j,k) == 0.) { + if (prob_type == 1) { //2D + if (rng[i + k*n_cell[2]] <= hardswitch_ratio){ + mat_alpha_arr(i,j,k) = mat_alpha_arr(i,j,k) * hardswitch_alpha_ratio; + } else { + mat_alpha_arr(i,j,k) = mat_alpha_arr(i,j,k); + } + } else if (prob_type == 2) { //3D + Real rand = Random(engine); + if (rand <= 0.04) { + mat_alpha_arr(i,j,k) = mat_alpha_arr(i,j,k) * 10.0; + } else { + mat_alpha_arr(i,j,k) = mat_alpha_arr(i,j,k); + } + } + } + }); + } + +} diff --git a/Source/Solver/TotalEnergyDensity.H b/Source/Solver/TotalEnergyDensity.H index 15a1326..35e219f 100644 --- a/Source/Solver/TotalEnergyDensity.H +++ b/Source/Solver/TotalEnergyDensity.H @@ -9,10 +9,28 @@ using namespace FerroX; void CalculateTDGL_RHS(Array &GL_rhs, Array &P_old, Array &E, - MultiFab& Gamma, + MultiFab& BigGamma, + MultiFab& alpha, + MultiFab& beta, + MultiFab& gamma, + MultiFab& g11, + MultiFab& g44, + MultiFab& g44_p, + MultiFab& g12, + MultiFab& alpha_12, + MultiFab& alpha_112, + MultiFab& alpha_123, MultiFab& MaterialMask, MultiFab& tphaseMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, const Geometry& geom, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); + +void Compute_P_Sum(const std::array& P, Real& sum); + +void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count); + +void Compute_P_av(const std::array& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real &P_av_z); diff --git a/Source/Solver/TotalEnergyDensity.cpp b/Source/Solver/TotalEnergyDensity.cpp index 53ef0bc..b31bf08 100644 --- a/Source/Solver/TotalEnergyDensity.cpp +++ b/Source/Solver/TotalEnergyDensity.cpp @@ -6,15 +6,38 @@ void CalculateTDGL_RHS(Array &GL_rhs, Array &P_old, Array &E, - MultiFab& Gamma, + MultiFab& BigGamma, + MultiFab& alpha, + MultiFab& beta, + MultiFab& gamma, + MultiFab& g11, + MultiFab& g44, + MultiFab& g44_p, + MultiFab& g12, + MultiFab& alpha_12, + MultiFab& alpha_112, + MultiFab& alpha_123, MultiFab& MaterialMask, MultiFab& tphaseMask, - MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta, + MultiFab& angle_alpha, + MultiFab& angle_beta, + MultiFab& angle_theta, const Geometry& geom, - const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { - // loop over boxes +// Real average_P_r = 0.; +// Real total_P_r = 0.; +// Real FE_index_counter = 0.; +// +// Compute_P_av(P_old, total_P_r, MaterialMask, FE_index_counter, average_P_r); +// +// //Calculate integrated electrode charge (Qe) based on eq 13 of https://pubs.aip.org/aip/jap/article/44/8/3379/6486/Depolarization-fields-in-thin-ferroelectric-films +// Real FE_thickness = FE_hi[2] - FE_lo[2]; +// Real one_m_theta = 1.0 - theta_dep; +// Real E_dep = average_P_r/(epsilon_0*epsilonZ_fe)*one_m_theta; +// + // loop over boxes for ( MFIter mfi(P_old[0]); mfi.isValid(); ++mfi ) { const Box& bx = mfi.validbox(); @@ -31,7 +54,17 @@ void CalculateTDGL_RHS(Array &GL_rhs, const Array4 &Ep = E[0].array(mfi); const Array4 &Eq = E[1].array(mfi); const Array4 &Er = E[2].array(mfi); - const Array4& Gam = Gamma.array(mfi); + const Array4& mat_BigGamma = BigGamma.array(mfi); + const Array4& mat_alpha = alpha.array(mfi); + const Array4& mat_beta = beta.array(mfi); + const Array4& mat_gamma = gamma.array(mfi); + const Array4& mat_g11 = g11.array(mfi); + const Array4& mat_g44 = g44.array(mfi); + const Array4& mat_g44_p = g44_p.array(mfi); + const Array4& mat_g12 = g12.array(mfi); + const Array4& mat_alpha_12 = alpha_12.array(mfi); + const Array4& mat_alpha_112 = alpha_112.array(mfi); + const Array4& mat_alpha_123 = alpha_123.array(mfi); const Array4& mask = MaterialMask.array(mfi); const Array4& tphase = tphaseMask.array(mfi); @@ -73,98 +106,167 @@ void CalculateTDGL_RHS(Array &GL_rhs, R_33 = cos(alpha_rad)*cos(beta_rad); } - Real dFdPp_Landau = alpha*pOld_p(i,j,k) + beta*std::pow(pOld_p(i,j,k),3.) + FerroX::gamma*std::pow(pOld_p(i,j,k),5.) - + 2. * alpha_12 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) - + 2. * alpha_12 * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),2.) - + 4. * alpha_112 * std::pow(pOld_p(i,j,k),3.) * (std::pow(pOld_q(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) - + 2. * alpha_112 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),4.) - + 2. * alpha_112 * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),4.) - + 2. * alpha_123 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); - - Real dFdPq_Landau = alpha*pOld_q(i,j,k) + beta*std::pow(pOld_q(i,j,k),3.) + FerroX::gamma*std::pow(pOld_q(i,j,k),5.) - + 2. * alpha_12 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) - + 2. * alpha_12 * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),2.) - + 4. * alpha_112 * std::pow(pOld_q(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) - + 2. * alpha_112 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),4.) - + 2. * alpha_112 * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),4.) - + 2. * alpha_123 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); + Real dFdPp_Landau = mat_alpha(i,j,k)*pOld_p(i,j,k) + mat_beta(i,j,k)*std::pow(pOld_p(i,j,k),3.) + mat_gamma(i,j,k)*std::pow(pOld_p(i,j,k),5.) + + 2. * mat_alpha_12(i,j,k) * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) + + 2. * mat_alpha_12(i,j,k) * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),2.) + + 4. * mat_alpha_112(i,j,k) * std::pow(pOld_p(i,j,k),3.) * (std::pow(pOld_q(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) + + 2. * mat_alpha_112(i,j,k) * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),4.) + + 2. * mat_alpha_112(i,j,k) * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),4.) + + 2. * mat_alpha_123(i,j,k) * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); + + Real dFdPq_Landau = mat_alpha(i,j,k)*pOld_q(i,j,k) + mat_beta(i,j,k)*std::pow(pOld_q(i,j,k),3.) + mat_gamma(i,j,k)*std::pow(pOld_q(i,j,k),5.) + + 2. * mat_alpha_12(i,j,k) * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) + + 2. * mat_alpha_12(i,j,k) * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),2.) + + 4. * mat_alpha_112(i,j,k) * std::pow(pOld_q(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) + + 2. * mat_alpha_112(i,j,k) * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),4.) + + 2. * mat_alpha_112(i,j,k) * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),4.) + + 2. * mat_alpha_123(i,j,k) * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); - Real dFdPr_Landau = alpha*pOld_r(i,j,k) + beta*std::pow(pOld_r(i,j,k),3.) + FerroX::gamma*std::pow(pOld_r(i,j,k),5.) - + 2. * alpha_12 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) - + 2. * alpha_12 * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),2.) - + 4. * alpha_112 * std::pow(pOld_r(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_q(i,j,k),2.)) - + 2. * alpha_112 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),4.) - + 2. * alpha_112 * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),4.) - + 2. * alpha_123 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_q(i,j,k),2.); - - Real dFdPp_grad = - g11 * DoubleDPDx(pOld_p, mask, i, j, k, dx) - - (g44 + g44_p) * DoubleDPDy(pOld_p, mask, i, j, k, dx) - - (g44 + g44_p) * DoubleDPDz(pOld_p, mask, i, j, k, dx) - - (g12 + g44 - g44_p) * DoubleDPDxDy(pOld_q, mask, i, j, k, dx) // d2P/dxdy - - (g12 + g44 - g44_p) * DoubleDPDxDz(pOld_r, mask, i, j, k, dx); // d2P/dxdz + Real dFdPr_Landau = mat_alpha(i,j,k)*pOld_r(i,j,k) + mat_beta(i,j,k)*std::pow(pOld_r(i,j,k),3.) + mat_gamma(i,j,k)*std::pow(pOld_r(i,j,k),5.) + + 2. * mat_alpha_12(i,j,k) * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) + + 2. * mat_alpha_12(i,j,k) * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),2.) + + 4. * mat_alpha_112(i,j,k) * std::pow(pOld_r(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_q(i,j,k),2.)) + + 2. * mat_alpha_112(i,j,k) * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),4.) + + 2. * mat_alpha_112(i,j,k) * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),4.) + + 2. * mat_alpha_123(i,j,k) * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_q(i,j,k),2.); + + Real dFdPp_grad = - mat_g11(i,j,k) * DoubleDPDx(pOld_p, mask, i, j, k, dx) + - (mat_g44(i,j,k) + mat_g44_p(i,j,k)) * DoubleDPDy(pOld_p, mask, i, j, k, dx) + - (mat_g44(i,j,k) + mat_g44_p(i,j,k)) * DoubleDPDz(pOld_p, mask, i, j, k, dx) + - (mat_g12(i,j,k) + mat_g44(i,j,k) - mat_g44_p(i,j,k)) * DoubleDPDxDy(pOld_q, mask, i, j, k, dx) // d2P/dxdy + - (mat_g12(i,j,k) + mat_g44(i,j,k) - mat_g44_p(i,j,k)) * DoubleDPDxDz(pOld_r, mask, i, j, k, dx); // d2P/dxdz - Real dFdPq_grad = - g11 * DoubleDPDy(pOld_q, mask, i, j, k, dx) - - (g44 - g44_p) * DoubleDPDx(pOld_q, mask, i, j, k, dx) - - (g44 - g44_p) * DoubleDPDz(pOld_q, mask, i, j, k, dx) - - (g12 + g44 + g44_p) * DoubleDPDxDy(pOld_p, mask, i, j, k, dx) // d2P/dxdy - - (g12 + g44 - g44_p) * DoubleDPDyDz(pOld_r, mask, i, j, k, dx);// d2P/dydz + Real dFdPq_grad = - mat_g11(i,j,k) * DoubleDPDy(pOld_q, mask, i, j, k, dx) + - (mat_g44(i,j,k) - mat_g44_p(i,j,k)) * DoubleDPDx(pOld_q, mask, i, j, k, dx) + - (mat_g44(i,j,k) - mat_g44_p(i,j,k)) * DoubleDPDz(pOld_q, mask, i, j, k, dx) + - (mat_g12(i,j,k) + mat_g44(i,j,k) + mat_g44_p(i,j,k)) * DoubleDPDxDy(pOld_p, mask, i, j, k, dx) // d2P/dxdy + - (mat_g12(i,j,k) + mat_g44(i,j,k) - mat_g44_p(i,j,k)) * DoubleDPDyDz(pOld_r, mask, i, j, k, dx);// d2P/dydz - Real dFdPr_grad = - g11 * ( R_31*R_31*DoubleDPDx(pOld_r, mask, i, j, k, dx) + Real dFdPr_grad = - mat_g11(i,j,k) * ( R_31*R_31*DoubleDPDx(pOld_r, mask, i, j, k, dx) +R_32*R_32*DoubleDPDy(pOld_r, mask, i, j, k, dx) +R_33*R_33*DoubleDPDz(pOld_r, mask, i, j, k, dx) +2.*R_31*R_32*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) +2.*R_32*R_33*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) +2.*R_33*R_31*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) - - (g44 - g44_p) * ( R_11*R_11*DoubleDPDx(pOld_r, mask, i, j, k, dx) + - (mat_g44(i,j,k) - mat_g44_p(i,j,k)) * ( R_11*R_11*DoubleDPDx(pOld_r, mask, i, j, k, dx) +R_12*R_12*DoubleDPDy(pOld_r, mask, i, j, k, dx) +R_13*R_13*DoubleDPDz(pOld_r, mask, i, j, k, dx) +2.*R_11*R_12*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) +2.*R_12*R_13*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) +2.*R_13*R_11*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) - - (g44 - g44_p) * ( R_21*R_21*DoubleDPDx(pOld_r, mask, i, j, k, dx) + - (mat_g44(i,j,k) - mat_g44_p(i,j,k)) * ( R_21*R_21*DoubleDPDx(pOld_r, mask, i, j, k, dx) +R_22*R_22*DoubleDPDy(pOld_r, mask, i, j, k, dx) +R_23*R_23*DoubleDPDz(pOld_r, mask, i, j, k, dx) +2.*R_21*R_22*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) +2.*R_22*R_23*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) +2.*R_23*R_21*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) - - (g44 + g44_p + g12) * DoubleDPDyDz(pOld_q, mask, i, j, k, dx) // d2P/dydz - - (g44 + g44_p + g12) * DoubleDPDxDz(pOld_p, mask, i, j, k, dx); // d2P/dxdz + - (mat_g44(i,j,k) + mat_g44_p(i,j,k) + mat_g12(i,j,k)) * DoubleDPDyDz(pOld_q, mask, i, j, k, dx) // d2P/dydz + - (mat_g44(i,j,k) + mat_g44_p(i,j,k) + mat_g12(i,j,k)) * DoubleDPDxDz(pOld_p, mask, i, j, k, dx); // d2P/dxdz - GL_RHS_p(i,j,k) = -1.0 * Gam(i,j,k) * + GL_RHS_p(i,j,k) = -1.0 * mat_BigGamma(i,j,k) * ( dFdPp_Landau + dFdPp_grad - - Ep(i,j,k) + - Ep(i,j,k) ); - GL_RHS_q(i,j,k) = -1.0 * Gam(i,j,k) * + GL_RHS_q(i,j,k) = -1.0 * mat_BigGamma(i,j,k) * ( dFdPq_Landau + dFdPq_grad - - Eq(i,j,k) + - Eq(i,j,k) ); - GL_RHS_r(i,j,k) = -1.0 * Gam(i,j,k) * + GL_RHS_r(i,j,k) = -1.0 * mat_BigGamma(i,j,k) * ( dFdPr_Landau + dFdPr_grad - - Er(i,j,k) + - Er(i,j,k) //+ E_dep ); if (is_polarization_scalar == 1){ - GL_RHS_p(i,j,k) = 0.0; - GL_RHS_q(i,j,k) = 0.0; - } + GL_RHS_p(i,j,k) = 0.0; + GL_RHS_q(i,j,k) = 0.0; + } //set t_phase GL_RHS_r to zero so that it stays zero. It is initialized to zero in t-phase as well //if(x <= t_phase_hi[0] && x >= t_phase_lo[0] && y <= t_phase_hi[1] && y >= t_phase_lo[1] && z <= t_phase_hi[2] && z >= t_phase_lo[2]){ if (tphase(i,j,k) == 1.0){ - GL_RHS_p(i,j,k) = 0.0; - GL_RHS_q(i,j,k) = 0.0; - GL_RHS_r(i,j,k) = 0.0; + GL_RHS_p(i,j,k) = 0.0; + GL_RHS_q(i,j,k) = 0.0; + GL_RHS_r(i,j,k) = 0.0; } }); } } +void Compute_P_Sum(const std::array& P, Real& sum) +{ + + // Initialize to zero + sum = 0.; + + ReduceOps reduce_op; + + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + for (MFIter mfi(P[2],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Box& bx_grid = mfi.validbox(); + auto const& fab = P[2].array(mfi); + + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple + { + return {fab(i,j,k)}; + }); + } + + sum = amrex::get<0>(reduce_data.value()); + ParallelDescriptor::ReduceRealSum(sum); +} + + +void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count) +{ + + // Initialize to zero + count = 0.; + + ReduceOps reduce_op; + + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + for (MFIter mfi(MaterialMask, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Box& bx_grid = mfi.validbox(); + + auto const& fab = MaterialMask.array(mfi); + + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple + { + if(fab(i,j,k) == 0.) { + return {1.}; + } else { + return {0.}; + } + + }); + } + + count = amrex::get<0>(reduce_data.value()); + ParallelDescriptor::ReduceRealSum(count); +} + +void Compute_P_av(const std::array& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real& P_av_z) +{ + Compute_P_Sum(P, sum); + Compute_P_index_Sum(MaterialMask, count); + P_av_z = sum/count; +} diff --git a/Source/main.cpp b/Source/main.cpp index 930caaf..7b3ccfe 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -66,7 +66,21 @@ void main_main (c_FerroX& rFerroX) // Ncomp = number of components for each array int Ncomp = 1; - MultiFab Gamma(ba, dm, Ncomp, Nghost); + MultiFab BigGamma(ba, dm, Ncomp, Nghost); + MultiFab alpha(ba, dm, Ncomp, Nghost); + MultiFab beta(ba, dm, Ncomp, Nghost); + MultiFab gamma(ba, dm, Ncomp, Nghost); + MultiFab epsilonX_fe(ba, dm, Ncomp, Nghost); + MultiFab epsilonZ_fe(ba, dm, Ncomp, Nghost); + MultiFab epsilon_de(ba, dm, Ncomp, Nghost); + MultiFab epsilon_si(ba, dm, Ncomp, Nghost); + MultiFab g11(ba, dm, Ncomp, Nghost); + MultiFab g44(ba, dm, Ncomp, Nghost); + MultiFab g44_p(ba, dm, Ncomp, Nghost); + MultiFab g12(ba, dm, Ncomp, Nghost); + MultiFab alpha_12(ba, dm, Ncomp, Nghost); + MultiFab alpha_112(ba, dm, Ncomp, Nghost); + MultiFab alpha_123(ba, dm, Ncomp, Nghost); Array P_old; for (int dir = 0; dir < AMREX_SPACEDIM; dir++) @@ -145,12 +159,57 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.setVal(0.); PoissonRHS.setVal(0.); tphaseMask.setVal(0.); + MaterialMask.setVal(0.); + BigGamma.setVal(0.); + alpha.setVal(0.); + beta.setVal(0.); + gamma.setVal(0.); + epsilonX_fe.setVal(0.); + epsilonZ_fe.setVal(0.); + epsilon_de.setVal(0.); + epsilon_si.setVal(0.); + g11.setVal(0.); + g44.setVal(0.); + g44_p.setVal(0.); + g12.setVal(0.); + alpha_12.setVal(0.); + alpha_112.setVal(0.); + alpha_123.setVal(0.); angle_alpha.setVal(0.); angle_beta.setVal(0.); angle_theta.setVal(0.); //Initialize material mask InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); + // Initialize material properties + Initialize_MaterialProperties(rFerroX, geom, + BigGamma, + alpha, + beta, + gamma, + epsilonX_fe, + epsilonZ_fe, + epsilon_de, + epsilon_si, + g11, + g44, + g44_p, + g12, + alpha_12, + alpha_112, + alpha_123); + + //check hardswitch_ratio and nucleation_ratio + if (hardswitch_ratio > 1.0 || nucleation_ratio > 1.0 || (hardswitch_ratio + nucleation_ratio) > 1.0) { + amrex::Abort("ERROR: hardswitch_ratio, nucleation_ratio, or their sum exceeds 1.0. Please check input parameters."); + } + + // define hard to switch spots with larger alpha value + if(hardswitch_flag == 1){ + printf("Set HardSwitch with ratio of:%g, with hardswitch_alpha_ratio:%g \n", hardswitch_ratio, hardswitch_alpha_ratio); + SetHardToSwitchNucleation(alpha, MaterialMask, n_cell, hardswitch_ratio, hardswitch_alpha_ratio); + } + //InitializeMaterialMask(rFerroX, geom, MaterialMask); if(Coordinate_Transformation == 1){ Initialize_tphase_Mask(rFerroX, geom, tphaseMask); @@ -178,7 +237,7 @@ void main_main (c_FerroX& rFerroX) beta_face[2].define(convert(ba,IntVect(AMREX_D_DECL(0,0,1))), dm, 1, 0);); // set cell-centered beta coefficient to permittivity based on mask - InitializePermittivity(LinOpBCType_2d, beta_cc, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi); + InitializePermittivity(LinOpBCType_2d, beta_cc, epsilonX_fe, epsilon_de, epsilon_si, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi); eXstatic_MFab_Util::AverageCellCenteredMultiFabToCellFaces(beta_cc, beta_face); // time = starting time in the simulation @@ -190,7 +249,7 @@ void main_main (c_FerroX& rFerroX) int linop_maxorder = 2; int amrlev = 0; //refers to the setcoarsest level of the solve - SetupMLMG(pMLMG, p_mlabec, LinOpBCType_2d, n_cell, beta_face, rFerroX, PoissonPhi, time, info); + SetupMLMG(pMLMG, p_mlabec, LinOpBCType_2d, n_cell, beta_face, P_old, MaterialMask, rFerroX, PoissonPhi, time, info); #ifdef AMREX_USE_EB std::unique_ptr p_mlebabec; @@ -200,8 +259,12 @@ void main_main (c_FerroX& rFerroX) // INITIALIZE P in FE and rho in SC regions //InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, geom, prob_lo, prob_hi);//old - InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi);//mask based - + InitializePandRho(P_old, BigGamma, charge_den, e_den, hole_den, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi);//mask based + if(nucleation_flag == 1){ + printf("Set Nucleation with ratio of:%g, with hardswitch_ratio:%g \n", nucleation_ratio, hardswitch_ratio); + SetNucleation(P_old, MaterialMask, n_cell, hardswitch_ratio, nucleation_ratio); + } + #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_old, charge_den, e_den, hole_den, MaterialMask, @@ -209,7 +272,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_old, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif // Calculate E from Phi @@ -220,7 +283,22 @@ void main_main (c_FerroX& rFerroX) { int plt_step = 0; WritePlotfile(rFerroX, PoissonPhi, PoissonRHS, P_old, E, hole_den, e_den, charge_den, beta_cc, - MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, Phidiff, geom, time, plt_step); + MaterialMask, tphaseMask, + BigGamma, + alpha, + beta, + gamma, + epsilonX_fe, + epsilonZ_fe, + epsilon_de, + epsilon_si, + g11, + g44, + g44_p, + g12, + alpha_12, + alpha_112, + alpha_123, angle_alpha, angle_beta, angle_theta, Phidiff, geom, time, plt_step); } amrex::Print() << "\n ========= Advance Steps ========== \n"<< std::endl; @@ -236,7 +314,7 @@ void main_main (c_FerroX& rFerroX) Real step_strt_time = ParallelDescriptor::second(); // compute f^n = f(P^n,Phi^n) - CalculateTDGL_RHS(GL_rhs, P_old, E, Gamma, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + CalculateTDGL_RHS(GL_rhs, P_old, E, BigGamma, alpha, beta, gamma, g11, g44, g44_p, g12, alpha_12, alpha_112, alpha_123, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); // P^{n+1,*} = P^n + dt * f^n for (int i = 0; i < 3; i++){ @@ -244,6 +322,7 @@ void main_main (c_FerroX& rFerroX) P_new_pre[i].FillBoundary(geom.periodicity()); } + #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_new_pre, charge_den, e_den, hole_den, MaterialMask, @@ -251,7 +330,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_new_pre, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif if (TimeIntegratorOrder == 1) { @@ -263,11 +342,15 @@ void main_main (c_FerroX& rFerroX) P_old[i].FillBoundary(geom.periodicity()); P_new_pre[i].FillBoundary(geom.periodicity()); } + if(nucleation_flag == 1){ + printf("Set Nucleation with ratio of:%g, with hardswitch_ratio:%g \n", nucleation_ratio, hardswitch_ratio); + SetNucleation(P_old, MaterialMask, n_cell, hardswitch_ratio, nucleation_ratio); + } } else { // compute f^{n+1,*} = f(P^{n+1,*},Phi^{n+1,*}) - CalculateTDGL_RHS(GL_rhs_pre, P_new_pre, E, Gamma, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + CalculateTDGL_RHS(GL_rhs_pre, P_new_pre, E, BigGamma, alpha, beta, gamma, g11, g44, g44_p, g12, alpha_12, alpha_112, alpha_123, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); // P^{n+1} = P^n + dt/2 * f^n + dt/2 * f^{n+1,*} for (int i = 0; i < 3; i++){ @@ -282,7 +365,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_new, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif // copy new solution into old solution @@ -315,7 +398,22 @@ void main_main (c_FerroX& rFerroX) { int plt_step = step; WritePlotfile(rFerroX, PoissonPhi, PoissonRHS, P_old, E, hole_den, e_den, charge_den, beta_cc, - MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, Phidiff, geom, time, plt_step); + MaterialMask, tphaseMask, + BigGamma, + alpha, + beta, + gamma, + epsilonX_fe, + epsilonZ_fe, + epsilon_de, + epsilon_si, + g11, + g44, + g44_p, + g12, + alpha_12, + alpha_112, + alpha_123, angle_alpha, angle_beta, angle_theta, Phidiff, geom, time, plt_step); } @@ -332,6 +430,7 @@ void main_main (c_FerroX& rFerroX) amrex::Print() << "step = " << step << ", Phi_Bc_hi = " << Phi_Bc_hi << ", num_Vapp = " << num_Vapp << ", sign = " << sign << std::endl; // Set Dirichlet BC for Phi in z + //SetPhiBC_z(PoissonPhi, P_old, MaterialMask, n_cell, geom); SetPhiBC_z(PoissonPhi, n_cell, geom); // set Dirichlet BC by reading in the ghost cell values @@ -348,7 +447,7 @@ void main_main (c_FerroX& rFerroX) #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, P_old, charge_den, e_den, hole_den, MaterialMask, - angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + angle_alpha, angle_beta, angle_theta, geom, n_cell, prob_lo, prob_hi); #endif }//end inc_step