Tmfp by atheulings · Pull Request #12 · eScatter/cstool

atheulings · 2017-06-30T12:58:33Z

Merge branch tmfp to master in order to implement the calculation of the electron range used in Geant4.

lvkessel · 2017-06-30T13:52:54Z

We are really calculating the inelastic imfp twice, which seems kind of pointless... but I suppose it fits in with the philosophy of cstool-is-a-library-not-an-application.

Anyway, I do have a few points:

It is probably nicer to use a function that simply integrates, without providing the icdf. Right now, you are using compute_tcs_icdf from cstool.compile, and throwing away the dummy icdf. That's OK for testing purposes, but unnecessary for production code.
Comments on lines 244, 245 are completely unrelated to the code that surrounds them. (OK, the relationship is a copy-paste thing, but the units don't even check out anymore)
Loop in lines 247 -- 257: this can be done in a more elegant / modern way. At least make sure that j, omega and omega_old only exist within the loop scope. For example, looping over for j in range(len(e_ran)), and breaking when e_ran[j] becomes too large, is generally considered "better". You could make things a lot nicer by using the syntax for j, omega in enumerate(e_ran). Then you still need to do something about omega_old, which is just e_ran[j-1] (unless j == 0). That relationship is not very clear from the code as it is right now.
Line 254: this can be replaced by +=.
You drop the units in lines 223 and 246, then reinsert them in line 263. Seems quite pointless (and very ugly) to me. Lines 227/228 should of course be multiplied by a length.

lvkessel · 2017-06-30T15:06:56Z

Also, just a general thing: why do we always use 129 values along the energy axis? If it's just an arbitrary choice, why not pick a round number?

"dv1" is not really an acronym. The "1" is because it is a first attempt. The D is probably for Delft, and the V could perhaps mean "Version".

atheulings · 2017-07-06T10:37:42Z

I have resolved all the comments on the code, please check this.

I do not know why we always use 129 values along the energy axis. Maybe @kerimarat or @slokhorst know?

lvkessel

OK, I did a more line-by-line style review here. These are mostly style-related suggestions, but in some cases I suggest a significantly different program flow. Do as you wish though, these are just opinions.

lvkessel · 2017-07-06T11:24:58Z

examples/cs.py

+            tmprange[i] = 1/tcs
+            for j, omega in enumerate(e_ran):
+                if (omega < K/2 and omega < K-s.band_structure.fermi):
+                    if j==0:


This may be more of a personal opinion, but I would prefer omega_old = 0 if j==0 else e_ran[j-1] here. These ternary operators tend to lead to excessively long lines, so one should be careful with them, but I think this is a case where it helps the code clarity.

lvkessel · 2017-07-06T11:25:05Z

examples/cs.py

+                s.elf_file.get_min_energy_interval())
+            tmprange[i] = 1/tcs
+            for j, omega in enumerate(e_ran):
+                if (omega < K/2 and omega < K-s.band_structure.fermi):


I would suggest the following:

if omega >= K/2 or omega >= K - s.band_structure.fermi: break

and no else clause for the actual code. This pattern is pretty common.

Again, this comment is mostly related to the code that comes after it.

lvkessel · 2017-07-06T11:25:07Z

examples/cs.py

+    tmprange = np.zeros(e_ran.shape) * units.m
+    print("# Computing inelastic electron range.")
+    for i, K in enumerate(e_ran):
+        ran_energies[i] = K.to('eV')


Why this? Does supplying data=e_ran.to('eV') in line 217 not already do this?

lvkessel · 2017-07-06T11:25:09Z

examples/cs.py

+        if ( K >= s.band_structure.barrier):
+            # make sure the tmfp is monotonously increasing for energies above
+            # the vacuum potential barrier
+            if (tmfp < tmfpmax):


Suggestion: tmfpmax = np.maximum(tmfp_max, tmfp) instead of these four lines and change line 239 to ran_tmfp[i] = tmfp_max.

So do not change the tmfp variable, conceptually it will always hold the tmfp for the current kinetic energy. Meanwhile, tmfp_max is the largest one you have found (bar anything below the vacuum barrier), and you make it clear that that's what you want to write to ran_tmfp.

"These four lines" relates to the line I commented on and the three lines BELOW it. github doesn't allow for multiline comments, sorry for the confusion.

lvkessel · 2017-07-06T11:25:12Z

examples/cs.py

    return compute_tcs_icdf(integrant, 0*units('rad'), np.pi*units('rad'), P)

+def compute_elastic_tcs(dcs, P):
+    def integrant(theta):


I agree with you that we should be consistent in spelling "integrand" incorrectly.

lvkessel · 2017-07-06T11:25:13Z

examples/cs.py

 from cstool.phonon import phonon_cs_fn
 from cstool.inelastic import inelastic_cs_fn
 from cstool.compile import compute_tcs_icdf
+from cstool.compile import compute_tcs


These two lines can be merged in one. Would make it more consistent with the other imports.

lvkessel · 2017-07-06T11:25:49Z

cstool/compile.py

+def compute_tcs(f, a, b, P, sampling=100000):
+    x = np.linspace(a, b.to(a.units), sampling) * a.units
+    y = f(x)
+    cf = np.r_[0, np.cumsum((x[1:] - x[:-1]) * (y[:-1] + y[1:]) / 2.0)]


You are still storing all intermediate integration results here, while you only need the last. Actually, since this is simply an n-point trapezoid rule, you could just return np.trapz(y, x) here -- shorter, intent is clear, and probably faster.

lvkessel · 2017-07-06T11:36:53Z

As for using 129 values, looking at the commit history it seems to have been introduced by @jhidding.

slokhorst · 2017-07-06T13:41:13Z

The 129 energy values are just a random choice I think. In e-scatter, we use 1024 values, so I guess it would make sense to match that here.

lvkessel · 2017-07-06T14:21:10Z

At least 1024 is a round number. That lets me sleep better than 129.

As the lower energy bounds are different between cstool and e-scatter, there is no point in using 1024 here to "match" e-scatter. But I think we will all agree that cstool should provide at least the resolution that our simulator intends to use. Anyway, we're getting a bit off-topic here.

slokhorst · 2017-07-13T14:10:24Z

Don't forget to fix #13 before/after merge

lvkessel · 2017-07-13T14:18:10Z

Absolutely, thanks for creating that issue. If nobody else fixes it, I'll do it as soon as I get back. It's independent from the tmfp thing though.

slokhorst · 2017-07-13T15:17:00Z

@lvkessel You're right. This shouldn't cause a merge conflict, so I've fixed it in the master branch: d111a3d.

…anch model.

jhidding · 2019-09-30T13:59:23Z

At least 1024 is a round number. That lets me sleep better than 129.

As the lower energy bounds are different between cstool and e-scatter, there is no point in using 1024 here to "match" e-scatter. But I think we will all agree that cstool should provide at least the resolution that our simulator intends to use. Anyway, we're getting a bit off-topic here.

Historical reference: The 129 was probably chosen as (128 + 1), creating more or less round numbers, since the end-point of the interval is also included in the range.

atheulings added 2 commits June 29, 2017 10:48

Write transport_mfp to hdf5.

edf02f4

Write electron_range to hdf5.

27dce18

atheulings requested a review from lvkessel June 30, 2017 12:58

Don't remove units

36ab8cc

Luc van Kessel and others added 3 commits July 4, 2017 17:09

Add "dv1" model.

46be971

"dv1" is not really an acronym. The "1" is because it is a first attempt. The D is probably for Delft, and the V could perhaps mean "Version".

Add function that returns only tcs and no icdf.

1ebdd6b

Clean up code.

ee994b7

lvkessel reviewed Jul 6, 2017

View reviewed changes

atheulings added 7 commits July 21, 2017 13:02

Update material parameters and add alumina.

e9623a7

Add reference for acoustic deformation potential of silicon.

e7b7556

Multiply single phonon model with 3 as first order approx. of dual br…

bcaf77b

…anch model.

Phonon parameters from Landolt-Bornstein.

c08a5f3

Clean-up silicon.yaml material file.

a0fd5d0

Correct DOS mass for silicon.

1342c39

Update README.rst

7ff2314

Comments

Conversation

atheulings commented Jun 30, 2017

Uh oh!

lvkessel commented Jun 30, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lvkessel commented Jun 30, 2017

Uh oh!

atheulings commented Jul 6, 2017

Uh oh!

lvkessel left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lvkessel Jul 6, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lvkessel commented Jul 6, 2017

Uh oh!

slokhorst commented Jul 6, 2017

Uh oh!

lvkessel commented Jul 6, 2017

Uh oh!

slokhorst commented Jul 13, 2017

Uh oh!

lvkessel commented Jul 13, 2017

Uh oh!

slokhorst commented Jul 13, 2017

Uh oh!

jhidding commented Sep 30, 2019

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

lvkessel commented Jun 30, 2017 •

edited

Loading

lvkessel Jul 6, 2017 •

edited

Loading