/DFT_pnp module | toolbox |=================================== DFT_pnp ================================== | | front end procedures: | | - build_pnp | - solve_pnp | 200 dict dup begin | DFT_pnp |=============================== build PNP system ============================ | DFT_window pnppdict | -- /build_pnp { /pnppdict name |-------------- permanent locals pnppdict begin |-- minimal window is 3 dup 3 lt { pop 3 } if /Nwindow name |-- count number of equations channel /Poisson_counter Nmobile Nconfined add /l array def _channel 0 Poisson_counter copy pop /NPoisson 1 def 0 1 Nconfined 1 sub { Nmobile add /kspecies name selfFac kspecies get 1e-10 gt { NPoisson Poisson_counter kspecies put /NPoisson NPoisson 1 add def } if } for /Neqns NPoisson Nmobile add Nconfined add def |-- Jacobian dimensions /Nbands Nwindow 3 add Neqns mul 1 sub def /JacobianRank Ngrid 2 sub Neqns mul def end |-- globals channel [ /CL /CR /FLUX ] { Nmobile /d array def } forall [ /Es /Us ] { [ Nmobile Nconfined add { Ngrid /d array } repeat ] def } forall /Vs [ NPoisson { Ngrid /d array } repeat ] def _channel |-- permanent locals pnppdict begin [ /W_F /W_K ] { Ngrid /d array def } forall [ /bt_1 /bt1 /wB_1 /wB1 /w_C_1 /w_C1 ] { Ngrid 2 sub /d array def } forall /Jacobian [ Nbands { JacobianRank /d array } repeat ] def /Lower [ Nbands 2 div { JacobianRank /d array } repeat ] def [ /RHS /Delta ] { JacobianRank /d array def } forall /LU_indices JacobianRank /l array def /tEX_hs Ngrid /d array def Nconfined 1 ge { /Zmat [ JacobianRank { Nconfined /d array } repeat ] def /ZT [ Nconfined { JacobianRank /d array } repeat ] def /Hmat [ Nconfined { Nconfined /d array } repeat ] def /tempmat [ Nconfined { Nconfined /d array } repeat ] def /VT [ Nconfined { JacobianRank /d array } repeat ] def /col_vec [ JacobianRank { 1 /d array } repeat ] def /small_col_vec1 [ Nconfined { 1 /d array } repeat ] def /small_col_vec2 [ Nconfined { 1 /d array } repeat ] def /aux_indices Nconfined /l array def } if end |-- initialize fixed difference coefficients h1 bt1 copy h add bt_1 copy h mul -1 pwr pop bt1 h1 mul -1 pwr pop } bind def |-- choppers /h { H interior } bind def /h1 { H 2 H length 2 sub getinterval } bind def /wa { WA interior } bind def /wa_1 { WA 0 Ngrid 2 sub getinterval } bind def /wa1 { WA 2 Ngrid 2 sub getinterval } bind def /wb { WB interior } bind def /wc { WC interior } bind def /wd { WD interior } bind def /w_F_1 { W_F 0 Ngrid 2 sub getinterval } bind def /w_F { W_F interior } bind def /w_F1 { W_F 2 Ngrid 2 sub getinterval } bind def /w_K { W_K interior } bind def /w_K_1 { W_K 0 Ngrid 2 sub getinterval } bind def /w_K1 { W_K 2 Ngrid 2 sub getinterval } bind def /x { X interior } bind def /x1 { X 2 Ngrid 2 sub getinterval } bind def /x_1 { X 0 Ngrid 2 sub getinterval } bind def /area { Area interior } bind def /eps { EPS interior } bind def /v { V interior } bind def /v_1 { V 0 v length getinterval } bind def /v1 { V 2 v length getinterval } bind def |================================== solve_pnp_dft ============================ | | - co-solves the PNP and DFT to determine voltage, concentration and excess | chemical potential profiles under flux. Particle fluxes and net current | are computed | | - expects an estimate of V and Us | | Execute in a 'stopped' context. Normal termination signals convergence | (that is, the maximal change in any element of the electric and (log) | concentration arrays is smaller than the 'tolerance'). If this does | not occur in 'maxiter' iterations, a 'stop' is executed. |------ defaults of control variables /maxiter 100 def /verbose true def /tolerance 1e-7 def /dampinglimit 1.0 def /solve_pnp_dft { /PoissonF epsilon0 kT_e mul e0 div Avogadro div SystemLength dup mul div def prime_Jac_invar comp_Cs pchem comp_EXs _pchem comp_Es /Newtoncount 0 def |------------------- Newton iteration false maxiter { /Newtoncount Newtoncount 1 add def pchem comp_dEXs _pchem build_Jacobian_RHS |-- solve the Newton system Jacobian Lower Nbands 2 div LU_indices bandLU not { (\n ** bandLU **) toconsole false exit } if RHS bandBS Delta copy pop Nconfined 1 ge { solve_aux_probs } if |-- find maximal absolute change in solution (V or Us) RHS maximum /RHSmax name verbose { RHSmax _ pop } if |-- backstep if change is greater than 1 log unit RHSmax dampinglimit gt { Delta RHSmax div dampinglimit mul pop } if |-- chicken out if change is excessive | RHSmax 1000.0 ge { stop } if |---------------- update V and Us using Newton result in Delta 0 1 NPoisson 1 sub { /kPoisson name Vs kPoisson get interior Delta kPoisson Neqns wa extract add pop } for 0 1 Nmobile Nconfined add 1 sub { /kspecies name Us kspecies get interior Delta kspecies NPoisson add Neqns wa extract add pop } for |-- re-compute Cs, EXs, and Es comp_Cs pchem comp_EXs _pchem comp_Es |-- exit test for Newton iteration: change of V, Us < tolerance RHSmax tolerance lt { pop true exit } if } repeat | Newton iteration not { stop } if |-- compute fluxes comp_fluxes verbose { Newtoncount _ pop } if } bind def |---------------------------- compute Cs /comp_Cs { 0 1 Nmobile 1 sub { /kspecies name Us kspecies get Cs kspecies get copy exp pop } for 0 1 Nconfined 1 sub { Nmobile add /kspecies name Us kspecies get kspecies confinedrange getinterval Cs kspecies get kspecies confinedrange getinterval copy exp pop } for } bind def |------------------------------ compute Es /comp_Es { 0 1 Nmobile 1 sub { /kspecies name Vs 0 get Es kspecies get copy Z kspecies get mul E0s kspecies get add EXs kspecies get add pop } for 0 1 Nconfined 1 sub { Nmobile add /kspecies name Vs Poisson_counter kspecies get get /V name /interv { kspecies confinedrange getinterval } def 0.0 Es kspecies get copy pop V interv Es kspecies get interv copy Z kspecies get mul E0s kspecies get interv add EXs kspecies get interv add pop } for } bind def |----------------------------- compute fluxes and net current /comp_fluxes { 0.0 0 1 Nmobile 1 sub { /kspecies name Us kspecies get 0 get Es kspecies get 0 get add exp Us kspecies get last get Es kspecies get last get add exp sub Es kspecies get WA copy exp Ds kspecies get div Area div DX mul integrateRS last get div Avogadro mul dup FLUX kspecies put Z kspecies get mul add } for e0 mul channel /NetCurrent name _channel } bind def |============================ Supporting procedures ========================= |---------------------------- solve auxiliary problems for confined particles | | When there are P confined particles there are P Ax=b problems to solve | as outlined in Press pgs. 75-77. Before calling this routine, | Press Eq. (2.7.20) has already been solved and the answer put into Delta. | In this routine, Eqs. (2.7.17) are solved using the already LU-factored | "Jacobian" and RHS as the right-hand side for all aux problems. | At the end of the routine, the solutions to the aux problems are combined | to give the actual solution to Jacobain * x = residues. The result | is copied into BOTH Delta and RHS, since both are needed to finish | the Newton iteration. /solve_aux_probs { |-- solve Eqs. (2.7.17) for Z matrix 0 1 Nconfined 1 sub { Nmobile add /kspecies name /left_idx kspecies confinedrange pop def /interv { kspecies confinedrange getinterval } def |-- determine u's 0.0 RHS copy pop 0.0 WA copy pop beta_max WA left_idx put wa RHS kspecies NPoisson add Neqns dilute pop |-- solve for Z's Jacobian Lower Nbands 2 div LU_indices RHS bandBS ZT kspecies Nmobile sub get copy pop |-- determine v's 0.0 WA copy pop Area interv WA interv copy DX interv mul Cs kspecies get interv mul pop WA Avogadro mul pop WA left_idx get 1.0 sub WA left_idx put 0.0 VT kspecies Nmobile sub get copy wa exch kspecies NPoisson add Neqns dilute pop } for Zmat ZT mattranspose pop tempmat VT Zmat matmul pop 0 1 tempmat length 1 sub { /k name tempmat k get k get 1.0 add tempmat k get k put } for tempmat aux_indices Hmat invertLU not { (\n ** invertLU **) toconsole exit } if pop col_vec Zmat small_col_vec1 Hmat small_col_vec2 VT Delta make_col_vec matmul matmul matmul make_arr_RHS neg Delta add Delta copy pop } bind def |-- convert array with length of RHS into matrix with 1 column (col_vec) /make_col_vec { /vec name 0 1 vec length 1 sub { /k name vec k get col_vec k get 0 put } for col_vec } bind def |-- put column vector with length of RHS into RHS array /make_arr_RHS { /vec name 0 1 vec length 1 sub { /k name vec k get 0 get RHS k put } for RHS } bind def |---------------------------- prime Jacobian invariables | At this time, these invariables include coefficients for the Poisson | equation that need not be recomputed inside the Newton iteration. /prime_Jac_invar { Area WA copy EPS mul pop | A wa_1 wB_1 copy wa add area div bt_1 mul PoissonF mul pop | wB_1 wa1 wB1 copy wa add area div bt1 mul PoissonF mul pop | wB1 bt1 maximum /beta_max name | for later scaling of confined particle eqns } bind def |---------------------------- build Jacobian and right-hand side of | Newton system /build_Jacobian_RHS { /kdiag Nbands 2 div def |-- null Jacobian and RHS arrays 0.0 RHS copy pop 0.0 Jacobian { copy } forall pop |-- Poisson(s) for mean field and fields of confined particles 0 1 NPoisson 1 sub { /kPoisson name Vs kPoisson get /V name 0.0 wa copy pop | sum[z*exp(u)] 0 1 Nmobile Nconfined add 1 sub { /kspecies name /omega 0.0 kPoisson 1 ge Poisson_counter kspecies get kPoisson eq and { pop selfFac kspecies get } if def Cs kspecies get interior wc copy 1.0 omega sub mul Z kspecies get mul 0 kPoisson kspecies NPoisson add to_Jac wa wc add pop } for wB_1 -1 kPoisson kPoisson to_Jac wB_1 wc copy wB1 add neg 0 kPoisson kPoisson to_Jac wB1 1 kPoisson kPoisson to_Jac wB_1 wb copy v_1 mul wc v mul add wB1 wc copy v1 mul add wa add Nmobile Nconfined add 1 Nspecies 1 sub { /kspecies name Cs kspecies get interior wc copy Z kspecies get mul add } for RHS kPoisson Neqns dilute pop } for | end Poisson loop |-- Nernst-Planck's for mobile particles 0 1 Nmobile 1 sub { /kspecies name |-- make w_F, w_C, w_K |-- scale the NP equations to order 1 by dividing each by the | minimum of diffusion coefficient * area /Fscale 1.0 0.0 Ds kspecies get W_F copy Area mul extrema pop def W_F Cs kspecies get mul Fscale div pop w_F_1 w_C_1 copy w_F add bt_1 mul pop w_F1 w_C1 copy w_F add bt1 mul pop Us kspecies get W_K copy Es kspecies get add pop |-- Jacobian elements w_C_1 wc copy w_C1 add neg pop w_C_1 wa copy Z kspecies get mul -1 kspecies NPoisson add 0 to_Jac wc wa copy Z kspecies get mul 0 kspecies NPoisson add 0 to_Jac w_C1 wa copy Z kspecies get mul 1 kspecies NPoisson add 0 to_Jac w_K_1 wb copy w_K sub bt_1 mul w_F_1 mul w_C_1 add -1 kspecies NPoisson add dup to_Jac w_C_1 -1 fill_dEXs w_K_1 wb copy w_K sub bt_1 mul w_K1 wd copy w_K sub bt1 mul add w_F mul wc add 0 kspecies NPoisson add dup to_Jac wc 0 fill_dEXs w_K1 wb copy w_K sub bt1 mul w_F1 mul w_C1 add 1 kspecies NPoisson add dup to_Jac w_C1 1 fill_dEXs |-- residue elements wc w_K mul w_C_1 w_K_1 mul add w_C1 w_K1 mul add RHS kspecies NPoisson add Neqns dilute pop } for | end mobile species loop |-- Boltzmann's for confined particles 0 1 Nconfined 1 sub { Nmobile add /kspecies name /left_idx kspecies confinedrange pop def |-- RHS elements /interv { kspecies confinedrange exch 1 sub exch getinterval } def |-- Us and Es are 0 outside the confining interval Us kspecies get W_K copy Es kspecies get add WA copy pop wa interv w_K_1 interv sub pop /interv { kspecies confinedrange getinterval } def Area interv WB interv copy DX interv mul Cs kspecies get interv mul integrateRS last get Avogadro mul total_confined kspecies Nmobile sub get sub WA left_idx put wa beta_max mul RHS kspecies NPoisson add Neqns dilute pop |-- Jacobian elements /interv { kspecies confinedrange exch 1 add exch 1 sub getinterval } def 0.0 WA copy pop Z kspecies get WA interv copy pop wa beta_max mul 0 kspecies NPoisson add Poisson_counter kspecies get to_Jac wa neg -1 kspecies NPoisson add Poisson_counter kspecies get to_Jac 0.0 WA copy pop beta_max neg WA interv copy pop wa -1 kspecies NPoisson add dup to_Jac beta_max WA copy pop wa 0 kspecies NPoisson add dup to_Jac beta_max 0 fill_dEXs beta_max neg -1 fill_dEXs } for | end confined species loop RHS neg pop } bind def |------- fill in dEX terms | factor offset | -- /fill_dEXs { /ofs name /factor name 0 1 Nmobile Nconfined add 1 sub { /kspecies2 name dEXs kspecies get kspecies2 get /lanes name lanes length /nlanes name /kcenter nlanes 2 div def kcenter neg 1 kcenter { /klane name 0.0 WA copy pop klane getlane kspecies isconfined { carve_confined WA interv copy } { carve_mobile wa copy } ifelse factor mul pop wa klane ofs add kspecies NPoisson add kspecies2 NPoisson add to_Jac } for } for } bind def |------- get a lane from dEX: lane_ofs | Lane /getlane { kcenter add lanes exch get } bind def |------- carve out shifted interior subarrays | Array shift | sub_array /carve_mobile { ofs 1 add Ngrid 2 sub getinterval } bind def /carve_confined { kspecies confinedrange exch 1 add ofs add exch 1 sub getinterval } bind def |----- Jacobian insertion | array l i j | -- | | l - negative or positive (node offset) | i - 0,1,.. (equation index) | j - 0,1,.. (index of independent variable) | | Tiles the 'array' into a bandarray of the Jacobian. /Jaccount 0 def /to_Jac { /jJ name /iJ name /lJ name /arr name /kband kdiag Neqns lJ mul add jJ add iJ sub def lJ 0 lt { /fband Neqns lJ abs mul iJ add def /arr arr lJ abs arr length lJ abs sub getinterval def } if lJ 0 eq { /fband iJ def } if lJ 0 gt { /fband iJ def /arr arr 0 arr length lJ sub getinterval def } if arr Jacobian kband get fband Neqns dilute_add pop /Jaccount Jaccount 1 add def } bind def |========================================================================= end _module