CRTM Monthly Meeting Protocol
Core Topic of the Meeting: Monte Carlo Radiative Transfer by Igor Polonsky (AER)
Date: 2021-09-30 Time: 15:00 EST
Location: Google hangout
Invited Speakers: Igor Polonsky (AER)
Meeting Chair: Benjamin Johnson (JCSDA)
Keeper of the Minutes: Patrick Stegmann (JCSDA)
Attendees: Patrick Stegmann, Cheng Dang, Igor Polonsky, Benjamin Johnson, Andrew Tangborn, Bryan Karpowicz, Cheng Dang, Haidao Lin, Haixia Liu, Hongli Wang, Jim Jung, Kevin Garrett, Ming Chen, Mingjing Tong, Nick Nalli, Quanhua Liu, Andrew Collard, Quanhua Liu, Sarah Lu, Scott Sieron, Shih-Wei Wei, yangqiu, Yingtao Mao
Agenda Item 1: | V2.4.1 Release Planning: bug fixes, features and testing |
Discussion: | Ben: Just a brief CRTM update. We’re delayed for the 2.4.1 release, but we want to make sure that we have a stable scalar release. Currently we have a compilation error with Intel Fortran that randomly pops up. Igor: Have you tried compiling it with gfortran in debug mode? Ben: We had a software engineer working with us for a while who was on the developer call with the Intel team and they seemed to agree that it’s a compiler issue. Patrick, which version did you use? Patrick: Version 18.0.3 on S4 and version 19.1.3 on Discover. Bryan: Dan mentioned this to me a while ago - that ADA module - it happens to some of the MATMULs in the ADA module. The compiler has issues with that. Ben: I make a note to myself about that. Quanhua: How about gfortran? Ben: That works fine. Igor: gfortran handles memory differently and you can insist on strict language. Ben: I have some time tomorrow and I can sort out tomorrow where it lies. Quanhua: I can also check out the Hera machine. Ben: We’ll make an issue for this bug. The other thing is that we want to switch away from the autotools build and go to a pure cmake build. Quanhua: Is the binary file support still available? Ben: Yes, but you need netCDF as a build requirement. Igor: Actually, I just looked into ADA and my advice is not to put addition in the multiplication operation. The compiler might run out of memory. Quanhua: I remember in the transmittance module somewhere in ODPS there is a strange subroutine that causes problems. Ben: The other subroutine with problems is the ODPS predictor routine. Yingtao: That’s a separate issue. Quanhua: Perhaps not. It could tell you the error comes from routine A, but it’s actually caused by routine B. Igor: Intel has this strange idea of getting all arrays into the stack. So they have the flag to get the arrays into the heap. Ben: Yes, there are things like heap32.
Ben: There is one more issue I that I came across with the autotools build in the configuration step. It can’t find the system netCDF .so file. It doesn’t happen on all machines and it doesn’t happen consistently. Those are the things that I am going to deal with in the next few days. So, when do you think you want to contribute code? Quanhua: Maybe next week. First, I will use Github to get my work there and then merge v3.0 into v.2.4.1. Ben: There will be a lot of manual merging conflicts you’ll have. We’ll send you some instructions Mark on what commands you need to use.
|
Result: | N/A |
Tasks: | N/A |
Responsible People: | N/A |
Deadline: | N/A |
Agenda Item 2: | V2.4.0 vs. v2.3.0 Comparison Update |
Discussion: | Ben: About a month ago, Scott has done some O-B comparisons to see if the impact of the version change from CRTM v2.3 to v2.4. Is that all you have right now, or is there something more you want to talk about? Scott: We were looking at the cloud fraction. Ben: I will post the link for the corresponding issue in the chat. I think the differences look very good. Quanhua: So, the update between the CRTM versions is just the surface? Ben: There’s also the reflection correction and cloud fraction. Haixia: Can I interrupt? Scott, so the control you are using, using CRTM 2.3.0 means we are using cloud fraction = 1. When we are using v2.4.0, we have to provide the GFDL cloud fraction. It’s not overcast cloud fraction anymore. What Scott just mentioned is that we will do another experiment to force v2.4.0 to have a cloud fraction of 1, so we can have an apple to apple comparison. Ben: Yes, that is the cloud fraction fix we made. Quanhua: So, that means the impact is likely coming from the cloud fraction, not the model. Ben: Scott, can you do map-like plots so we can identify geographical regions of interest? Scott: I’m not quite sure. Quanhua: What Ben means is that right now you are showing the global distribution. Andrew: Some tools exist for radiosondes, but they may be unequally distributed. Ben: I guess I am a little bit confused. Could you make a temperature map at 800 mbar? Andrew: That’s straightforward. Mingjing has a comment. Mingjing: I have an initial evaluation with our shield model. I also forced the cloud fraction to one and for the conventional observation the differences are not large, but with v2.4.0 more observations are assimilated. However, it seems to be a positive signal for the new version. But for the cloud channels the bias is larger. There is more positive bias. Maybe you can help me understand. Ben: If you have some plots, that might be helpful to understand. Mingjing: Yes, after I fix the surface problem in our model. Ben: Igor was scheduled to give a talk but he had to leave early. |
Result: | N/A |
Tasks: | N/A |
Responsible People: | N/A |
Deadline: | N/A |
Agenda Item 3: | Coefficient generation update by Patrick Stegmann |
Discussion: | Ben: Patrick, do you have updates on the coefficient generation? Patrick: Yes, we did have the first coefficient generation user tutorial on Wednesday. The tutorial was focused on MW instruments, with the goal of enabling CRTM users to create their own CRTM MW instrument coefficients. As a follow-up, all participants will also receive a brief overview with concise instructions on how follow the coefficient generation workflow. The next tutorial in the series will focus on using the CRTM Rosenkranz utility as a standalone and for comparisons with e.g. monoRTM. After that we will have a tutorial on IR instrument coefficient generation, but this will be more involved because of the LBLRTM issues. Lastly, Yingtao has agreed to work on the NLTE task and perhaps he wants to share some updates? Yingtao: Yes, I am working on restoring the NLTE code in the repository and will likely provide a delivery next week or the week after. I was thinking that the NLTE task only includes new coefficients and not restoring the code, thus the delay.
|
Result: | N/A |
Tasks: | N/A |
Responsible People: | N/A |
Deadline: | N/A |
Agenda Item 4: | Updates by Cheng |
Discussion: | Ben: Cheng, is there anything you want to talk about? Cheng: I had an issue regarding the versioning of the aerosol coefficients. I thought this was something we could implement in v.2.4.1. If there are no further comments I will just stick with a specific solution. Ben: An issue is also the attribution because of the RTTOV problem where people computed their own coefficients and blamed RTTOV for their deficiancies. The metadata will contain the version of the CRTM for which it was developed. Quanhua: v2.4.1 means release number for the CRTM? Ben: Yes. Yingtao: I saw that Dang Cheng suggested a solution where you have the starting version of the CRTM and the range of CRTM models for which it works. Ben: Yes, so it gets to be a bit tricky. There’s a couple of issues here. The source files have a specific version inside them. The other part is the version of the tables itself. We may need to go through a few internal revisions. Yingtao: There are a lot of numbers to go through for the user. We could have a configuration file with a CRTM release. Cheng: For the current control we have a kind of release management. I think this is still functional. But for example, if we have a new table, should we just name it after the new release? I think we should use the CRTM release number as a threshold for the coefficient data. Or maybe we can include something to help us determine that this is the correct table. |
Result: | N/A |
Tasks: | N/A |
Responsible People: | N/A |
Deadline: | N/A |
Agenda Item 5: | Open discussions |
Discussion: | Ben: Any other topics that people would like to discuss? Quanhua: One question, do we need to have some kind of user utilities with the CRTM. E.g. in the past, in order to compute the surface emissivity, the user needs to compute some internal variables from the CRTM. In case of the reflectance this is more complicated. Thus, we could think about a user utility, and set optional parameters to compute that. Yingtao: I agree. For instance, for the TOA reflectance, this is not part of the core radiative transfer utility. So, we can add some tools for that. Quanhua: Yes, or for people who want weighting functions. Ben: The Navy uses pyCRTM for that quite extensively. Bryan: Yeah, NRL uses it pretty regularly. Yingtao: Bryan, does your pyCRTM include the standard profile for that? Bryan: What I do is using the check_crtm.py as a template. Yingtao: Some people just want to use the standard profile as an input, but it’s not easily possible. We implemented a user utility to do that to have an easy test. Bryan: Sure. Yingtao: Some people just want to play with it. Ben: That’s sort of the spirit of pyCRTM. I will get some patches from Ed. Quanhua: Does the JCSDA have their own website to allow people to run the CRTM online? Ben: The answer is yes, but that doesn’t exist right now. I would love to do that but our team is so small that we need to focus on core applications. There will be a JEDI application and there are two groups working on near-realtime applications on the JCSDA website. Quanhua: Ten years ago, the CRTM had this functionality, but due to the NOAA security policy we were not allowed to have a PHP machine online. Ben: I did something similar during my PhD but currently we don’t have the bandwidth. Yingtao: That would be a good functionality, to quickly look at the weighting functions for a specific instrument. Ben: Any other topics? One question about the GOES-T ABI coefficients? Yingtao: GOES-18? Last time there was a channel that was a little bit different between our result and the one from Wisconsin. Quanhua: That’s GOES-17. Yingtao: GOES-18 probably also has the same problem. Haixia: Just a follow-up question about GOES-18, what channel has differences? Yingtao: Channel 12, ozone channel? Haixia: No, I am talking about GOES-18. Yingtao: That’s similar. Haixia: Do we have an updated launch date for GOES-17? Ben: January 8 according to GOES-R.gov. But it’s probably going to be delayed because of the shutdown. Igor is back, do you want to do your presentation now?
|
Result: | N/A |
Tasks: | N/A |
Responsible People: | N/A |
Deadline: | N/A |
Agenda Item 6: | Science Update from Igor Polonsky |
Discussion: | Igor: Yes, it’s going to be about Monte Carlo in atmospheric optics. I’ll give you some history about Monte Carlo. You can interrupt if you have questions. We talk about why Monte Carlo is so popular and useful. Then I talk about how Monte Carlo can estimate Jacobians. I’ll close with some references on free codes. I didn’t make any conclusions. So, Monte Carlo historically is just derived from the fact that many physical phenomena can be described through statistical estimation. The name of Monte Carlo was suggested by Nicolas Metropolis who worked at Los Alamos. These methods were essential in the 50s to develop the hydrogen bomb. What is a simple application of Monte Carlo? The computation of integrals, e.g. when you want to compute a surface. Bryan: Igor, your slides aren’t advancing. Igor: Ah, ok, I will do it outside of application mode. Anyway, we can compute the integral under a curve by putting a lot of points inside the curve. According to a simple analysis, the error is just the inverse square root of the number of points. To increase accuracy by 10 times, you need to increase the number of photons by 100. This was a real problem 10 years ago. Sometimes it took a few days to run these calculations on supercomputers. How can we apply Monte Carlo to real physics? On a quantum level, all processes are random. That’s why the most trivial and straightforward application of Monte Carlo to atmospheric optics is just to follow nature. We trace photon trajectories from sources like the sun or emission. Then we through dice to see if it gets absorbed or scattered and we include refraction and Brouillion and Raman scattering. It’s very simple to add new physics. The biggest problem is actually specifying new physics. For more complicated problems with dependence of properties on many coordinates, you need to evaluate extinction or optical thickness in the direction of propagation. The more complicated your geometry is, the more complicated is the photon tracing. In simple cases the implementation is straightforward. We just emulate some randomness. We just check how far the photon goes before it gets absorbed. We convert optical to geometric distance and that is the biggest slowdown for Monte Carlo computations. For any event we ask for a random number. Technically your Monte Carlo code should include all complexities. In Monte Carlo, the implementation of polarization for instance is trivial. You just add polarization to the photon state. Adding physics to Monte Carlo is not complicated. Adding Zeeman for example is simple, it means your photon changes frequency at the throw of a dice. Why is the Monte Carlo method still not that accurate compared to a real measurement? When we have a measurement, we have a very high number of photons from the sun, but in a Monte Carlo code, we can only handle a much smaller number of photons, which means that the error is proportional to that. Quanhua: Igor, the photon number you are talking about is the number received by the instrument, not from the sun? Igor: Yes, but the photons in the detector also come from the sun. The Monte Carlo simulation just has a much smaller number of photons. There are ways to mitigate this problem. The direct Monte Carlo method is not efficient because a lot of photons are wasted by not getting measured by the instrument. The method here is called local estimate. The method does not behave well in a statistical sense, but it converges to the true result. Here we would like to focus on trajectories that make a contribution to the detector. At every point of our trajectory we assign a weight to the photon and check if the photon makes a contribution to the detector. And in this method as I said we assign a photon weight. Now if the dice show that the photon is going to be absorbed, this means that the absorption just eats part of the photon, not the whole photon. And at any point we just trace the contribution to the detector. Finally, Jacobians. The problem with them is that with Monte Carlo you can’t use Finite Differences. Your estimation of differences will include parts, one difference from the photon path and the other from the actual changes in the media. So, there are two methods to estimate Jacobians with Monte Carlo. One method is called correlated sampling. In the first method you use the same photon trace for both methods. When you use weights, you use the same weights for the perturbed case and you can use finite differences. Lastly, there are some 1D and 3D codes available online. I would like to recommend some books, first the one by Marchuck et al. 1980, second by Dupree and Fraley, 2002. This one is more focused on particle physics and for nuclear reactors. I also added a paper by Deutschmann on Monte Carlo with Jacobian estimation. This is my short introduction to Monte Carlo. Ben: Thank you. There is also another category of models, called reverse Monte Carlo. But you can’t keep track of polarization changes. Igor: Allow me to disagree. First you have the adjoint method based on the reciprocity theorem. The angular width of the sun is also very limited, so you have about the same problem. The adjoint method is very powerful when you have a large source, like an instrument looking at the earth. When you put it this way, Monte Carlo really cries for any symmetries you can apply. The biggest challenge of Monte Carlo is to trace only useful photons. Ming Cheng: I have a few questions. There are different random number generators. In your experience, is there a difference between machines? Igor: Thank you for the question. One time at Los Alamos when I was working with Anthony Davis on Lidar, his code had a problem. That time I found out that not all random numbers are made equal. First, all random number generators you use need to have a long period. Ming: That question is, when you generate this kind of random number, do you have some idea of the probability model that best fits your current model. Igor: For Monte Carlo, the more you know, the better. If you have a prior estimate of your probability distribution you should use it. You can also create a Monte Carlo method that also updates your probability distribution. Ming: What’s your experience with moving from 1D to 3D? Igor: It depends on what kind of physical property you would like to measure. If it’s fluxes, you just need a small number of photons. In that sense it’s similar to the discrete ordinate method. For an instrument radiance you need more photons. Quanhua: For thermal emission of the earth, Monte Carlo is fine with the backtrace method. I can quickly show the screen what we did 10 years ago. Do you see the paper? This is a 3D model we used. You do not need many photons there. For the polarization you can use the backtrace method to find where it was absorbed/emitted. Then you can do a forward trace to account for the polarization and it’s very accurate. For example, here, we compared Monte Carlo and a 1D model and the difference is either 0.1 or 0.2. But in the past the computers were too slow for satellite processing. Ming: Do you think scattering becomes very complicated? Quanhua: No, this model is very easy for the user, similar to the CRTM. If you want to know the polarization, you just use the thermal trace. The limitation is that it’s only for thermal emission. Igor: Let me disagree. If you have to deal with an SRF you will spend a lot of time and you are dealing with a really large source, so you are implicitly using the adjoint method. Quanhua: We call this the backward-forward Monte Carlo method. Backward for the absorption and forward for the polarization state. Ben: I posted a paper by Kummerow in the link that also describes this method. Quanhua: I have contact with Kummerow because he worked with our university at the time. We were probably a bit earlier. Ben: We are pretty far past the time, just real quick, Igor do you have a plan for Monte Carlo? Igor: Just for fun. Quanhua: Igor, if you are interested and have time I think it’s a good topic. For a 2D satellite image, the Monte Carlo model could be faster. Igor: There is one complication: Monte Carlo is usually monochromatic and we can’t apply the CRTM transmittance approach to Monte Carlo. It’s non-trivial. I would be more comfortable with something like OSS. Quanhua: If you do the forward model in any case for each small box of a cubical grid. Many photons don’t go to the CRTM. But many photons might go to another pixel of a 2D image. That’s where Monte Carlo could be fast. Yingtao: Do you want to use Monte Carlo as a benchmark? Quanhua: The Monte Carlo people already looked at the geometric impact of 3D effects. Ben: That’s relevant for aircraft. Igor: Then you can look at turbulence. Yingtao: And bulk cloud effects. Igor: There were papers in the 90s that related the cloud structure and turbulence to the radiation field. Ben: We could also combine machine learning and Monte Carlo. Quanhua: Maybe for today’s technology Monte Carlo for a 2D image might not be slower than a 1D model. Igor: I disagree. The geometrical inhomogeneity will require a huge look up table. Ben: Alright, I think we’ll wrap it up. |
Result: | N/A |
Tasks: | N/A |
Responsible People: | N/A |
Deadline: | N/A |
16:51 EST End of meeting.
Google Hangouts Chat Comments:
Benjamin Johnson - NOAA Affiliate
3:06 PM
https://github.com/JCSDA-internal/crtm/issues/160
You
3:22 PM
I think Igor needs to leave the meeting early today at 3:30.
Benjamin Johnson - NOAA Affiliate
3:24 PM
https://github.com/NOAA-EMC/GSI/issues/188
Igor Polonsky
3:31 PM
have to leave. I have another telecon to join
Benjamin Johnson - NOAA Affiliate
3:32 PM
Thanks Igor!
Cheng Dang
3:40 PM
https://github.com/JCSDA-internal/crtm/issues/231
Benjamin Johnson - NOAA Affiliate
4:04 PM
brian you might want to just interrupt, he might not be able to see hands.
bryan*
Nick Nalli - NOAA Affiliate
4:04 PM
What slide are we on?
Benjamin Johnson - NOAA Affiliate
4:39 PM
Kummerow and Roberti did something similar
Sarah Lu - NOAA Affiliate
4:40 PM
Thank Igor for the presentation. Sorry I need to leave the meeting now.
Benjamin Johnson - NOAA Affiliate
4:40 PM
http://rain.atmos.colostate.edu/research/pubs/roberti1999.pdf
thanks Sarah
yanqiu
4:46 PM
has to leave. thanks for the presentation!