Set aerosol distribution params with command line arguments

[claresinger]

This PR addresses issue #116. It currently includes the addition of the blk2m microphysics, so that PR should be merged first and I will update this one with any changes.

One question @pdziekan, what is the difference between dry_distros.emplace and dry_distros.push_back. I found one example in the libcloudph++ kinematic setup that uses push_back and sets two modes with different kappa.

[pdziekan]

@claresinger the difference between push_back and emplace is subtle, both can be used to add elements to dry_distros: push_back adds a copy of the argument to the container emplace constructs a new element in the container, arguments are passed to the constructor

You can define distributions with different kappas by emplacing multiple elements to dry_distros.

[claresinger]

Problem solved. It was an issue with units!

Ignore everything below here

Thanks @pdziekan. I am getting this compile error currently.

/home/travis/build/igfuw/UWLCM/src/opts/opts_blk_2m.hpp:64:5: error: no matching function for call to ‘std::vector<libcloudphxx::blk_2m::opts_t<float>::lognormal_mode_t, std::allocator<libcloudphxx::blk_2m::opts_t<float>::lognormal_mode_t> >::push_back(<brace-enclosed initializer list>)’
2382     rt_params.cloudph_opts.dry_distros.push_back({

This is where I try and initialize the aerosol distribution for bulk_2m. https://github.com/claresinger/UWLCM/blob/a5b0abe2227d1a2359627f113d337ed0254239ae/src/opts/opts_blk_2m.hpp#L50-L76

I was using the example from libcloudph as a template https://github.com/igfuw/libcloudphxx/blob/b88d610774a5113925463171efb7046aa8a07180/models/kinematic_2D/src/opts_blk_2m.hpp#L49-L60

Is the problem that I am using push_back or is there a different cause for this error? Thanks for your help!

[claresinger]

@pdziekan In order for the tests to pass, I can't output these aerosol distribution parameters from user_params, but I would like to include them in the const.h5 file. I think we just need to update the ref files in the tests, but I'm not sure how to do that or if I should leave that to you. Thanks! https://github.com/claresinger/UWLCM/blob/3582442c8ba8f3d8375227ce45d6a002f163267e/src/solvers/slvr_common.hpp#L137-L146

      /* TODO: need to update ref files in tests to include this in output
      this->record_aux_const("mean_rd1", "aerosol_dist_params", params.user_params.mean_rd1 / si::metres);  
      this->record_aux_const("sdev_rd1", "aerosol_dist_params", params.user_params.sdev_rd1);
      this->record_aux_const("n1_stp", "aerosol_dist_params", params.user_params.n1_stp * si::cubic_metres);
      this->record_aux_const("kappa1", "aerosol_dist_params", params.user_params.kappa1);
      this->record_aux_const("mean_rd2", "aerosol_dist_params", params.user_params.mean_rd2 / si::metres);  
      this->record_aux_const("sdev_rd2", "aerosol_dist_params", params.user_params.sdev_rd2);
      this->record_aux_const("n2_stp", "aerosol_dist_params", params.user_params.n2_stp * si::cubic_metres);
      this->record_aux_const("kappa2", "aerosol_dist_params", params.user_params.kappa2);
      */

[pdziekan]

@claresinger you can update redata yourself:

0. Clear the "output" directory in api tests build directory.
1. Run only api iles tests with ctest -R api_test_iles
2. The api_test_iles_diff test will fail, because new consts are not found in refdata const.h5. It is good to look at the output of the api_test_iles_diff test to make sure that the only differences found are the ones expected (i.e. lack of new constants in refdata). Unfortunately, output of the h5diff command, which is used to find differences, is cluttered by messages about variables that have no differences, so it is difficult to spot the ones that actually are different.
3. Copy contents of the "output" directory to refdata_iles
4. Repeat steps 0-3 for the Smagorinsky api tests (replace _iles with _smg)

[claresinger]

@pdziekan - circling back to how to emplace multiple dry_distros. You said "You can define distributions with different kappas by emplacing multiple elements to dry_distros." What is the syntax for this? I tried this, but it doesn't seem to work correctly. Thanks!

      rt_params.cloudph_opts_init.dry_distros.emplace(
        user_params.kappa1,
        std::make_shared<setup::log_dry_radii<thrust_real_t>> (
          user_params.mean_rd1,
          thrust_real_t(1.0e-6) * si::metres,
          user_params.sdev_rd1,
          thrust_real_t(1.2),
          user_params.n1_stp,
          thrust_real_t(0) / si::cubic_metres
        )
      );
      rt_params.cloudph_opts_init.dry_distros.emplace(
        user_params.kappa2,
        std::make_shared<setup::log_dry_radii<thrust_real_t>> (
          thrust_real_t(1.0e-6) * si::metres,
          user_params.mean_rd2,
          thrust_real_t(1.2),
          user_params.sdev_rd2,
          thrust_real_t(0) / si::cubic_metres,
          user_params.n2_stp
        )
      );

[pdziekan]

@claresinger you've coded it correctly. Why do you think it does not work as expected?

To diagnose this issue, you can use the diag_kappa_rng(min,max) function to calculate moments only over super-droplets with kappa between min and max. For example, you can calculate concentration of aerosols (in 1/kg) with kappa between 0.5 and 1 with: diag_kappa_rng(0.5,1); diag_wet_mom(0);

[claresinger]

@pdziekan, I say it doesn't work because when I run a simulation with the aerosol initialized as above vs. like this, I get different results (when kappa1 = kappa2 -- of course otherwise I'd expect a different result).

      rt_params.cloudph_opts_init.dry_distros.emplace(
        user_params.kappa1,
        std::make_shared<setup::log_dry_radii<thrust_real_t>> (
          user_params.mean_rd1,
          user_params.mean_rd2,
          user_params.sdev_rd1,
          user_params.sdev_rd2,
          user_params.n1_stp,
          user_params.n2_stp
        )
      );

If I try and emplace two distributions with the same kappa will that cause an issue?

It seems this is the case. When you try and emplace a second distribution with the same kappa it doesn't get added.

emplacing two dist with different kappa
kappa1, n1 (cm-3) =  [0.61] [125.]
kappa2, n2 (cm-3) =  [0.611] [1.e-06]
mean nd (cm-3) =  102.04923
mean nd for kappa < 0.5 (cm-3) =  0.0
mean nd for kappa > 0.5 (cm-3) =  102.04923

emplacing two dist with same kappa to mimic dycoms default
kappa1, n1 (cm-3) =  [0.61] [125.]
kappa2, n2 (cm-3) =  [0.61] [65.]
mean nd (cm-3) =  102.004845
mean nd for kappa < 0.5 (cm-3) =  0.0
mean nd for kappa > 0.5 (cm-3) =  102.004845

Here are two short simulations I did. In the first I use kappa1=0.61 and kappa2=0.611 so they are slightly different, but not in a physically relevant way. In the second they are both kappa=0.61. For the first simulation I set the second distribution to only have 1e-6 cm-3 vs. a typical concentration in the second simulation of 65 cm-3. The resulting number concentration of dry particles are nearly identical.

[pdziekan]

@claresinger opts_init.dry_distros is a map, so it is not possible to put two elements with the same key (kappa).

Sorry for not pointing this out earlier.

[claresinger]

@pdziekan The tests pass for me locally, but are failing on travis. The only difference is coming from the hgt_fctr in const.h5 where locally I get -0 (which I then put in the refdata) and on travis I get -3.8e-40 which is essentially zero. But I guess since it is computing the relative difference, that is still 1 (since (x-0)/x=1). Do you know how to avoid this issue with h5diff?

dataset: </hgt_fctr> and </hgt_fctr>
size:           [4]           [4]
position        hgt_fctr        hgt_fctr        difference      relative       
------------------------------------------------------------------------
[ 3 ]          -3.863183684e-40 -0              3.863183684e-40 1             
1 differences found

[claresinger]

@pdziekan I just want to confirm, that I won't have the same problem even if kappa1 = kappa2 when using bulk 2m because the chem_b parameter is not used as a key for the map. Is that correct? Thanks! https://github.com/claresinger/UWLCM/blob/add_cmdline_aerosols/src/opts/opts_blk_2m.hpp

  if (user_params.n1_stp * si::cubic_metres >= 0 && user_params.n2_stp * si::cubic_metres >= 0) {
    rt_params.cloudph_opts.dry_distros.push_back({
      .mean_rd = user_params.mean_rd1 / si::metres,
      .sdev_rd = user_params.sdev_rd1,
      .N_stp   = user_params.n1_stp * si::cubic_metres,
      .chem_b  = user_params.kappa1
    });
    rt_params.cloudph_opts.dry_distros.push_back({
      .mean_rd = user_params.mean_rd2 / si::metres,
      .sdev_rd = user_params.sdev_rd2,
      .N_stp   = user_params.n2_stp * si::cubic_metres,
      .chem_b  = user_params.kappa2
    });

[claresinger]

Other than this problem with hgt_fctr and the tests, everything looks good. I've confirmed that the results from a short simulation when I let the aerosol parameters default to the values specified by the case vs. when I set them to be equal to those values give identical results.

[pdziekan]

@pdziekan I just want to confirm, that I won't have the same problem even if kappa1 = kappa2 when using bulk 2m because the chem_b parameter is not used as a key for the map. Is that correct? Thanks!

Yes, that's correct.

@pdziekan The tests pass for me locally, but are failing on travis. The only difference is coming from the hgt_fctr in const.h5 where locally I get -0 (which I then put in the refdata) and on travis I get -3.8e-40 which is essentially zero. But I guess since it is computing the relative difference, that is still 1 (since (x-0)/x=1). Do you know how to avoid this issue with h5diff?

Which build configuration you used when producing new refdata? On travis, CMAKE_BUILD_TYPE=RelWithDebInfo and you should use the same config to make refdata. This may affect the precision of math functions and may be the reason why you get 0 vs -3.8e-40 on travis. In particular, Debug configuration could give more precise values, i.e. 0.

If this doesn't help, then I don't have any simple solution to the h5diff testing. In h5diff you can set tolerance for absolute error with the --delta option. I don't know what happens if you use both --relative and --delta in one h5diff call. I expect that it is either not permitted or that both conditions need to be satisfied, which is not what we want here. You can try combining two h5diff calls with the logical OR operator || in bash. One call would test for relative difference in all variables except hgt_fctr. You can exclude variables with the --compare option. The other call would test for absolute difference in hgt_fctr only. You can specify variable to be tested with: h5diff [OPTIONS] file1 file2 hgt_fctr

Other than this problem with hgt_fctr and the tests, everything looks good. I've confirmed that the results from a short simulation when I let the aerosol parameters default to the values specified by the case vs. when I set them to be equal to those values give identical results.

Nice!

[claresinger]

@pdziekan Good catch about the build type. I think it's ready to merge now if it looks good to you!

[claresinger]

@pdziekan -- The changes in refdata are from the blk_2m tests and also from changing the const.h5 file to include the information about the aerosol size distribution.

One other note, there is some weird randomness with when the moist_thermal tests pass. The diff test sometimes fails because of differences in the cloud fraction. But then just by re-triggering the build it will sometimes pass instead. I'm not sure what is causing this, maybe an issue with the tolerance just being too small. But it would be nice if the tests would reproducibly either pass or fail.

[pdziekan]

@claresinger yes, I know that randomness in the moist_thermal test is annoying. It would be much nicer to have exactly reproducible results. Unfortunately, even for the same rng seed, different compilers/CUDA versions/build types/number of threads/etc. give different results. That's why there is some tolerance for randomness in this test. Feel free to add a PR increasing tolerance whenever you see the test fail.