Skip to content

SCRATCHPAD.md

Latest commit

 

History

History
211 lines (166 loc) · 18 KB

SCRATCHPAD.md

File metadata and controls

211 lines (166 loc) · 18 KB
Raw

Intended Use

  • NOTE: This is a temporary storage space for todo items, notes, code, etc
  • It is regularly cleared, do not keep important stuff here!

[0.60177399 0.99 0.01 ]

Issues to be filed

  • connect up "get_formula" method for solution phases
  • Future issue: figure out if we can deactivate ordering in olivine when Mn, Ni, Co are missing
  • enki-community as a subrepo of ThermoEngineLite to increase visibility (ask Cian)

Simplified Data for equilibrium calculations

  • Create CompactEquilState
  • Examine how History works in magmaforge (wrt EquilState objects)
  • Write conversion method to extract CompactEquilState from EquilState
  • Retrieve and save CompactEquilState list for history
  • Save history to json
  • Convert history into pandas dataframe
  • Save HistoryTable as csv

MagmaForge structure

  • split primary quantities and dependent quantities (e.g., density)

  • split primary quantities and user output

  • choice of independent variables-- S vs T, P vs V

    • P, T, extensive moles of endmembers for every phase
    • internally switch where needed
    • can have multiple iterations of a phase existing at once
    • fO2 is a dependent user quantity-- can calculate from PTX but process may have variables associated (buffered/not, which buffer)

  • other process variables-- frac/equil xtal, how much to keep

  • we're thinking that structured arrays are our best bet

  • gather results stage at end to put everything in 2D structured array

    • and then pandas for i/o interface layer

  • async operation in magmaforge-- progress bars and partial results

  • metadata for phases

  • transaction receipts / ledger

  • list of entitities as phase-samples in the simulation

    • each entity has a different size based on number of components
    • BUT, we could store them in pools of same size (e.g. all pure together, all ternary together etc.)
    • Does this interact at all with missing components? (for the future)

  • test endmembers of Opx on server and ThermoEngineLite

    • endmember of Opx are identical to Cpx in ThermoEngineLite... adjustments to H, S, (and presumably V) are never made

  • look at taylor_soln in notebook to see why adjustments are not made...

Store in Issues later?

Pigeonite?!?

  • do we need to incorporate pigeonite param values into a 3rd pyroxene phases???

Endmember comparison script

  • when comparing with objc, we need to test pure endmembers as well

Managing missing components in coder code

  • rewind coder to previous state w/o (non-generalized) modifications made for Olivine

  • Rebuild phase library with coder changes

  • assess failing tests situation

    • only missing component tests failed

  • identify a specific failing coder test with missing components

  • add code via strings to solution notebook for Olivine (and have coder write them to files)

  • test effect for olivine (does it pass now)

  • add code via strings to solution notebook for Cpx ...

  • Deactivate nudge comp and show that equlibrate begins to fail

  • Adjust docstrings for ordering functions in add_expressions_to_model and OrderingFunctions class

  • change ordering functions tuple to dictionary input and rebuild every phase

Notes on enforcing non-zero site mole fractions using coder

  • define cation site list here
  • Add test for site_fractions exists (prevent compile error)
  • define nudge_scale
  • for each cation site:
  • if < TINY:
  •  [x] eval cation-site-dependence on n's

  •  [x] nudge all n's according to deriv 

  •     e.g. n[0] += nudge_scale*deriv[0] ... n[6]
  • reset n1-n6 values to their nudged values from updated n[] array

Managing missing components in equilibrate code

  • modify C-code in coder_templates for missing components in complex solutions
  • in coder.py: skip calculations for missing components in complex solutiosn
  • test new phases in equilibrate
  • repeat for simple solutions
    • NOTE we are currently skipping Feldspar missing component tests b/c they require rebuilding phase, and the pipeline won't do that for a week
  • Make it easy to trigger a rebuild of all phases
    • add commit-prefix [rebuild-cache] to manually trigger rebuilding all phases
  • Profile equilibrate as is
    • Cannot distinguish where in the algorithm it is thrashing
  • Lightly refactor equilibrate by extracting methods, so we can usefully profile algorithm
  • [/] Remaining coder issues with missing components
    • [/] deactivate nudge_composition
    • allow missing component as a feasible comp in coder
    • write olivine test w/ missing components (pure Fa)
    • show that order parameters must be solved w/ missing components accounted for in order param minimization
    • figure out how to manually set some order params to default values
    • Just modify Olivine notebook (for initial order guess)
    • [-] test in c code
    • [-] modify coder.py and coder_templates to accept order param changes

NOTES on code/data efficiency and cache usage

  • it appears that the code is engineered in a really inefficent way
  • inlining every function (done by coder/sympy) is resulting in really large functions that cannot easily be held in the cache and are causing non-neglible level of instruction cache misses
  • also, the data is not laid out optimally, so there are non-negligible data cache misses
  • Critically, the frontend cycles are idle 16% of the time, reflecting front-end starvation, where the instructions cannot be delivered fast enough (because they are so verbose and cache misses are eating time).
    • shrinking code size will help with this
  • perhaps even more critical are the lack of memory for order parameter values. They are recalculated from scratch every time. This is wasted computation, since they should be able start from their previous values to find the optimal ordering to minimize the energy
  • This would likely remove most of the order parameter solving time, which I think is significant for ordered phases
  • Also, only writing the gibbs energy function to jax, and then caclulating high order derivs using jax should further optimize the code
    • This would remove most of the compile time, reducing compile times to minutes instead of hours.
  • Direct implementation of stateless energy function would allow separate storage of order param values.
  • Then they can be merely refined at subsequent iterations
  • Refining data storage, and including order param values, will allow this recall of previous ordered state, and reduce ordering calculation burden.
  • Also, we could use exact symbolic expressions for following ordering manifold to fully avoid ordering calculation of derivative values (just like now), and thus we would only take the dgds and d2gds2 matrices and perform matrix operations to caculate high order derivs following manifold
    • this is same as it is done now, but the important distinction is that we would NOT rely on jax to numerically evalulate derivs and follow the manifold.
  • a well laid-out data structure for each phase would allow better storage of results and reanalysis later (e.g. in magmaforge)
    • we might (?) also store implied cation site fractions.
    • Or just make it really easy to recalculate them from the moles of endmembers and order params
  • uniformatize phase modeling systems:
    • by maintaining loops in the c-code, rather than unrolling every loop, the code for phases could be shared, even for phases with different numbers of components
  • We could migrate to an entity component system where phases are described as a combination of factors that contribute to the Gibbs energy
    • it is not clear whether this will be amazing or a bridge too far
    • How would this interact with sympy-driven approach?
      • can sympy even write out code including loops over arrays?
      • this may just not be necessary... wait for last
    • sharing of code across phases would shrink compile times
    • gibbs energy expressions would be a combination of reusable systems, e.g:
      • ideal mixing of endmembers
      • ideal mixing entropy of distinct units (components, species, or cations on site, does not matter)
      • regular excess mixing contribution to enthalpy
      • excess energy of compositional dependent species or of site mole species
      • endmember energy terms (not all need be present!, only small fraction have ordering or landau transition, but we currently check every phase to see if values are zero)
      • Possibly we can mix and match standard terms that apply across many contexts with terms specific to a particular phase, where generalization is unnecessary
    • optimizing data layout for these reusable systems, by grouping data types together into numpy structured array (akin to structs in C)
      • the equilibrate code could then consider structures of arrays (of these base structs), to allow cache optimization
      • We can split out arrays for present_phases and absent_phases, avoiding unnecessary loads
  • Parallelization:
    • this becomes possible across phases, for refining any order parameters
    • and then parallel computation across both energy terms and across phases

Steps

  • Use C-profiler on C-code to determine time spent in order_s function
  • try running sympy2jax printer
  • try taking deriv of non-ordered phase models
  • try derivs of jax code for ordered phase (ignoring internal order params)
  • validate numerical results for high order derivs from jax vs sympy code

Cache and CPU performance stats

Performance counter stats for 'python run_baseline_magmaforge.py':

  223974814772      cycles:u                                                      (30.02%)
  322738624247      instructions:u            #    1.44  insn per cycle         
                                              #    0.11  stalled cycles per insn  (30.04%)
     711874376      L1-icache-load-misses:u   #    2.38% of all L1-icache accesses  (30.06%)
   29911950697      L1-icache-loads:u                                             (30.07%)
    4149391508      L1-dcache-load-misses:u   #    3.67% of all L1-dcache accesses  (30.06%)
  112936612355      L1-dcache-loads:u                                             (30.08%)

LLC-loads:u
LLC-load-misses:u
787967121 cache-misses:u (30.10%) 294850626 branch-misses:u # 0.48% of all branches (30.10%) 176890627601 all_dc_accesses:u (30.10%) 0 all_tlbs_flushed:u (30.11%) 2675185991 l1_dtlb_misses:u (30.13%) 4147104926 l2_cache_accesses_from_dc_misses:u (30.11%) 11155398073 l2_cache_accesses_from_ic_misses:u (30.10%) 65438682 l2_pf_miss_l2_l3:u (30.10%) 18811502452 l2_latency.l2_cycles_waiting_on_fills:u (30.07%) 61515857378 branch-instructions:u (30.05%) 295308833 branch-misses:u # 0.48% of all branches (30.04%) 36245095855 idle-cycles-frontend:u # 16.18% frontend cycles idle (30.04%) 712187456 L1-dcache-prefetches:u (30.01%) 36365929417 stalled-cycles-frontend:u # 16.24% frontend cycles idle (30.00%)

  29.650411402 seconds time elapsed

  50.944646000 seconds user
   5.873210000 seconds sys

result = (1.0/4.0)R(-n1 - n2 - n3 - n4 - n5 - n6 - n7)(-4n6protected_log(n6/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 2(n3 + n4)(protected_log((n3 + n4)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2) - 2(protected_log((2n1 + 2n2 + n6 + 2n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(2n1 + 2n2 + n6 + 2*n7) +

2*(protected_log((n3 + n4 + 2n5 - n6)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(-n3 - n4 - 2*n5 + n6)

  • 2*(protected_log((n3 + n4 + 2n5 + n6)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(n3 + n4 + 2n5 + n6) - 2(protected_log((2n3 + 2n4 + 2n5 - n6)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 2M_LN2)(-2n3 - 2n4 - 2n5 + n6) + 2*(protected_log((2n1 + 2n2 + n3 + n4 + n6 + 2n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(2n1 + 2n2 + n3 + n4 + n6 + 2n7) + 2(protected_log((2n1 + 2n2 + n3 + n4 + 2n5 + 2n6 + 2n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(2n1 + 2n2 + n3 + n4 + 2n5 + 2n6 + 2n7) + 2(protected_log((-n1s2 - n2s2 - n3s2 - n4s2 - n5s2 - n6s2 - n7s2 + n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(n1s2 + n2s2 + n3s2 + n4s2 + n5s2 + n6s2 + n7s2 - n7) - 2(protected_log((n1s2 + n2s2 + n3s2 + n4s2 + n5s2 + n6s2 + n7s2 + n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(n1s2 + n2s2 + n3s2 + n4s2 + n5s2 + n6s2 + n7s2 + n7) + (protected_log((-2n1s1 - 2n2s1 - 2n3s1 - 2n4s1 - 2n5s1 + 2n5 - 2n6s1 + n6 - 2n7s1)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 2M_LN2)(2n1s1 + 2n2s1 + 2n3s1 + 2n4s1 + 2n5s1 - 2n5 + 2n6s1 - n6 + 2n7s1) - (protected_log((2n1s1 + 2n2s1 + 2n3s1 + 2n4s1 + 2n5s1 + 2n5 + 2n6s1 + n6 + 2n7s1)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 2M_LN2)(2n1s1 + 2n2s1 + 2n3s1 + 2n4s1 + 2n5s1 + 2n5 + 2n6s1 + n6 + 2n7s1) + 2(protected_log((-n1s2 - n2s2 + 2n2 - n3s2 - n4s2 - n5s2 - n6s2 - n7s2 - n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(n1s2 + n2s2 - 2n2 + n3s2 + n4s2 + n5s2 + n6s2 + n7s2 + n7) + (protected_log((-2n1s1 - 2n2s1 - 2n3s1 - 2n4s1 + 4n4 - 2n5s1 + 2n5 - 2n6s1 - n6 - 2n7s1)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 3M_LN2)(2n1s1 + 2n2s1 + 2n3s1 + 2n4s1 - 4n4 + 2n5s1 - 2n5 + 2n6s1 + n6 + 2n7s1) - (protected_log((2n1s1 + 2n2s1 + 2n3s1 + 4n3 + 2n4s1 + 2n5s1 + 2n5 + 2n6s1 - n6 + 2n7s1)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 3M_LN2)(2n1s1 + 2n2s1 + 2n3s1 + 4n3 + 2n4s1 + 2n5s1 + 2n5 + 2n6s1 - n6 + 2n7s1) - 2(protected_log((n1s2 + 2n1 + n2s2 + n3s2 + n3 + n4s2 + n4 + n5s2 + n6s2 + n6 + n7s2 + 3n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - M_LN2)(n1s2 + 2n1 + n2s2 + n3s2 + n3 + n4s2 + n4 + n5s2 + n6s2 + n6 + n7s2 + 3n7) - 4(n1 + n2 + n3 + n4 + n5)protected_log((n1 + n2 + n3 + n4 + n5)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + 4(n1 + n2 + n3 + n4 + n5 + n7)*protected_log((n1 + n2 + n3 + n4 + n5 + n7)/(n1 + n2 + n3 + n4 + n5 + n6 + n7)))/(n1 + n2 + n3 + n4 + n5 + n6 + n7)

  • (n1 + n2 + n3 + n4 + n5 + n6 + n7)(-1.0/8.0S27n3s2/(n1 + n2 + n3 + n4 + n5 + n6 + n7) - 1.0/8.0S27n4s2/(n1 + n2 + n3 + n4 + n5 + n6 + n7) - 1.0/4.0S27n5s2/(n1 + n2 + n3 + n4 + n5 + n6 + n7) - 1.0/8.0S27n6s2/(n1 + n2 + n3 + n4 + n5 + n6 + n7) - dS1 - 1.0/2.0dS27n2((n7)(n7))/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/8.0dS27n3((n7)(n7))/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/8.0dS27n3n7s2/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/8.0dS27n4*((n7)(n7))/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/8.0dS27n4n7s2/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/4.0dS27n5*((n7)(n7))/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/4.0dS27n5n7s2/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/8.0dS27n6*((n7)(n7))/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) - 1.0/8.0dS27n6n7s2/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + (1.0/4.0)dS27((n7)(n7)(n7))/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + (1.0/4.0)dS27((n7)(n7))s2/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + n2n7*(-1.0/2.0S27 + pS1 - pS2)/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + n2*(dS1 - dS2 - (endmember[0].dmu0dT)(T, P) + (endmember[1].dmu0dT)(T, P))/(n1 + n2 + n3 + n4 + n5 + n6 + n7) + n3n7(-1.0/8.0S27 + pS1 - pS3)/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + n3*(dS1 - dS3 - (endmember[0].dmu0dT)(T, P) + (endmember[2].dmu0dT)(T, P))/(n1 + n2 + n3 + n4 + n5 + n6 + n7) + n4n7(-1.0/8.0S27 + pS1 - pS4)/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + n4*(dS1 - dS4 - (endmember[0].dmu0dT)(T, P) + (endmember[3].dmu0dT)(T, P))/(n1 + n2 + n3 + n4 + n5 + n6 + n7) + n5n7(-1.0/4.0S27 + pS1 - pS5)/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + n5*(dS1 - dS5 - (endmember[0].dmu0dT)(T, P) + (endmember[4].dmu0dT)(T, P))/(n1 + n2 + n3 + n4 + n5 + n6 + n7) + n6n7(-1.0/8.0S27 + (1.0/2.0)pS1 + (1.0/2.0)pS3 - 1.0/2.0pS4 + (1.0/2.0)pS5 - pSjd)/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + n6((1.0/2.0)dS1 + (1.0/2.0)dS3 - 1.0/2.0dS4 + (1.0/2.0)dS5 - dSjd - 1.0/2.0(endmember[0].dmu0dT)(T, P) - 1.0/2.0(endmember[2].dmu0dT)(T, P) + (1.0/2.0)(endmember[3].dmu0dT)(T, P) - 1.0/2.0(endmember[4].dmu0dT)(T, P) + (endmember[5].dmu0dT)(T, P))/(n1 + n2 + n3 + n4 + n5 + n6 + n7) + ((n7)(n7))((1.0/4.0)S27 - 1.0/4.0dS27 + pS1 - pS7)/((n1 + n2 + n3 + n4 + n5 + n6 + n7)(n1 + n2 + n3 + n4 + n5 + n6 + n7)) + n7s2*((1.0/4.0)S27 - 1.0/4.0dS27)/(n1 + n2 + n3 + n4 + n5 + n6 + n7) + n7*(-1.0/4.0*S27 + dS1 - dS7 - pS1 - (*endmember[0].dmu0dT)(T, P) + (*endmember[6].dmu0dT)(T, P))/(n1 + n2 + n3 + n4 + n5 + n6 + n7) + (*endmember[0].dmu0dT)(T, P));

Try replacing this line in build solution: Sconfig_Na = R*Sconf_2cats_coupled(xM1_Al+xM1_Fe3+xM1_Ti, xM2_Na).subs(xM1_Al + xM1_Fe3 + xM1_Ti, 1-xM1_Fe2-xM1_Mg)