{"action":"created","issue":{"url":"https://api.github.com/repos/StrikeNP/trac_test/issues/301","repository_url":"https://api.github.com/repos/StrikeNP/trac_test","labels_url":"https://api.github.com/repos/StrikeNP/trac_test/issues/301/labels{/name}","comments_url":"https://api.github.com/repos/StrikeNP/trac_test/issues/301/comments","events_url":"https://api.github.com/repos/StrikeNP/trac_test/issues/301/events","html_url":"https://github.com/StrikeNP/trac_test/issues/301","id":330443994,"node_id":"MDU6SXNzdWUzMzA0NDM5OTQ=","number":301,"title":"Latin hypercube: Sample preferentially within cloud (Trac #301)","user":{"login":"StrikeNP","id":2577059,"node_id":"MDQ6VXNlcjI1NzcwNTk=","avatar_url":"https://avatars0.githubusercontent.com/u/2577059?v=4","gravatar_id":"","url":"https://api.github.com/users/StrikeNP","html_url":"https://github.com/StrikeNP","followers_url":"https://api.github.com/users/StrikeNP/followers","following_url":"https://api.github.com/users/StrikeNP/following{/other_user}","gists_url":"https://api.github.com/users/StrikeNP/gists{/gist_id}","starred_url":"https://api.github.com/users/StrikeNP/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/StrikeNP/subscriptions","organizations_url":"https://api.github.com/users/StrikeNP/orgs","repos_url":"https://api.github.com/users/StrikeNP/repos","events_url":"https://api.github.com/users/StrikeNP/events{/privacy}","received_events_url":"https://api.github.com/users/StrikeNP/received_events","type":"User","site_admin":false},"labels":[{"id":958019852,"node_id":"MDU6TGFiZWw5NTgwMTk4NTI=","url":"https://api.github.com/repos/StrikeNP/trac_test/labels/Migrated%20from%20Trac","name":"Migrated from Trac","color":"FFFFFF","default":false},{"id":958042830,"node_id":"MDU6TGFiZWw5NTgwNDI4MzA=","url":"https://api.github.com/repos/StrikeNP/trac_test/labels/clubb_src","name":"clubb_src","color":"FFFFFF","default":false},{"id":958042829,"node_id":"MDU6TGFiZWw5NTgwNDI4Mjk=","url":"https://api.github.com/repos/StrikeNP/trac_test/labels/dschanen@uwm.edu","name":"dschanen@uwm.edu","color":"FFFFFF","default":false},{"id":958019578,"node_id":"MDU6TGFiZWw5NTgwMTk1Nzg=","url":"https://api.github.com/repos/StrikeNP/trac_test/labels/enhancement","name":"enhancement","color":"a2eeef","default":true}],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":{"url":"https://api.github.com/repos/StrikeNP/trac_test/milestones/8","html_url":"https://github.com/StrikeNP/trac_test/milestone/8","labels_url":"https://api.github.com/repos/StrikeNP/trac_test/milestones/8/labels","id":3408876,"node_id":"MDk6TWlsZXN0b25lMzQwODg3Ng==","number":8,"title":"Improve SILHS","description":null,"creator":{"login":"StrikeNP","id":2577059,"node_id":"MDQ6VXNlcjI1NzcwNTk=","avatar_url":"https://avatars0.githubusercontent.com/u/2577059?v=4","gravatar_id":"","url":"https://api.github.com/users/StrikeNP","html_url":"https://github.com/StrikeNP","followers_url":"https://api.github.com/users/StrikeNP/followers","following_url":"https://api.github.com/users/StrikeNP/following{/other_user}","gists_url":"https://api.github.com/users/StrikeNP/gists{/gist_id}","starred_url":"https://api.github.com/users/StrikeNP/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/StrikeNP/subscriptions","organizations_url":"https://api.github.com/users/StrikeNP/orgs","repos_url":"https://api.github.com/users/StrikeNP/repos","events_url":"https://api.github.com/users/StrikeNP/events{/privacy}","received_events_url":"https://api.github.com/users/StrikeNP/received_events","type":"User","site_admin":false},"open_issues":65,"closed_issues":16,"state":"open","created_at":"2018-06-07T19:59:06Z","updated_at":"2018-06-08T15:37:38Z","due_on":null,"closed_at":null},"comments":5,"created_at":"2018-06-07T21:34:04Z","updated_at":"2018-06-08T15:38:52Z","closed_at":null,"author_association":"OWNER","body":"'''Introduction:'''\n\nThe results of LH sampling seem to converge, but they are too noisy for cases with small cloud fraction, such as RICO.  To avoid noise, we need to sample preferentially from within cloud when cloud fraction is small.\n\n'''Requirements:'''\n\nGenerate two columns: one within cloud at k_s_mellor and one within clear air at k_s_mellor.  From these columns, compute a weighted average.  Eventually, we'll need to generalize for any even number of sample points per grid box and timestep.  However, we'll forego the freedom to use odd numbers of samples with this option.\n\n'''Technical spec:'''\n\n1.  Keep the present code intact.  Enclose all new code in a logical flag, l_lh_cloud_weighted_sampling.  For now, we could assume that we'll only use this with l_lh_vert_overlap = .true., although later we may possibly wish to relax this.  Just in case, let's add some code like the following at the beginning of the latin_hypercube_sampling subroutine:\n\n```text\n    ! Sanity check for l_lh_cloud_weighted_sampling\n    if ( l_lh_cloud_weighted_sampling .and. mod( n_micro_calls, 2 ) /= 0 ) then\n      write(fstderr,*) \"Cloud weighted sampling requires micro calls to be divisible by 2.\"\n      stop \"Fatal error.\"\n    end if\n```\n\n2.  At the k_lh_start, if cloud_frac > 0.5 .or. cloud_frac < 0.01, then choose points randomly, as we do now.  (However, don't leave 0.01 as a magic number; define it as a variable.)\n\n\n3.  On the other hand, if at k_lh_start, 0.01 < cloud_frac < 0.5, then keep choosing points until you find one that lies in the cloudy portion.  E.g., suppose you choose mixture fraction 1 (using the d+1 variate).  Then, if X_u(1,...) > (1-cloud_frac1), then the sample point is within cloud.  If instead you had chosen mixture fraction 2, then use cloud_frac2.  If the point is not within cloud, then choose other values of X_u(d+1) and X_u(1).  When you find a point within cloud, stop.  (This rejection method might be slow; if so, we could choose a faster method.)\n\nInitially we're going to use the rejection method, so that we have something to compare against a more sophisticated weighted method.  We'll also want assertion tests to verify that the cloudy columns are actually cloudy and the clear columns are actually clear.\n  \n4.  Then use the rejection method to keep choosing points until you find one that lies in the clear portion.  \nTo do this we can modify ``` generate_uniform_sample ``` to use a new function called ```choose_permuted_random``` rather than using an inline formula, and this will allow us to generate just one element rather than n_micro_calls * d_variables elements.\nNew code:\n\n```text\n   ! Choose values of sample using permuted vector and random number generator\n    do j = 1,n_micro_calls\n      do k = 1,dp1\n        X_u_one_lev(j,k) = choose_permuted_random( nt_repeat, p_matrix(j,k) )\n      end do\n    end do\n\n!        print*, 'p_matrix(:,1)= ', p_matrix(:,1)\n!        print*, 'p_matrix(:,dp1)= ', p_matrix(:,dp1)\n!        print*, 'X(:,1)= ', X(:,:)\n\n    return\n  end subroutine generate_uniform_sample\n\n!----------------------------------------------------------------------\n  double precision function choose_permuted_random( nt_repeat, p_matrix_element )\n\n    use mt95, only: genrand_real3 ! Procedure(s)\n\n    use mt95, only: genrand_real ! Constants\n\n    implicit none\n\n    integer, intent(in) :: & \n      nt_repeat, &  ! Number of samples before the sequence repeats\n      p_matrix_element       ! Permuted integer\n\n    real(kind=genrand_real) :: rand ! Random float with a range of (0,1)\n\n    ! ---- Begin Code ----\n\n    call genrand_real3( rand ) ! genrand_real3's range is (0,1)\n\n    choose_permuted_random = (1.0d0/nt_repeat)*(p_matrix_element + rand )\n\n    return\n  end function choose_permuted_random\n```\nThen, when we have our X_u_all_levs, we can re-use ```choose_permuted_random``` on k_lh_start for step 5.\n\n5.  Once you have found the two sample points at k_lh_start, proceed as usual, choosing the full column of points.  \n\n6.  In the end, we need to weight the results of the two sample columns by cloud_frac and (1-cloud_frac).\nThis will be tricky.  It looks like we can use the code in ```l_cloud_weighted_averaging```, but in the past that code didn't work as expected.  We'll have to debug this I think.\n\n'''Tests:'''\n\n1.  Ensure that at any timestep, we always have one sample point within cloud and one out of cloud at level k_lh_start.\n\n2. Test whether the new LH sample converges to the analytic solution for KK microphysics for the RF02 and RICO cases.\n\n3. See whether the profiles vary less in time when (l_lh_cloud_weighted_sampling = .true.) than when (l_lh_cloud_weighted_sampling = .false.).  \n\n'''Time estimate:'''\n\nThe tricky part will probably be going from the rejection method to something more sophisticated.\n\n'''Special note'''\n\nWe should keep in mind that this undoes some of the Latin Hypercube sampling aspects of the code (i.e. we're overwriting s with an in-cloud value when t has already been determined using Latin Hypercube sampling).\n\nAttachments:\n\nhttp://carson.math.uwm.edu/trac/clubb/attachment/ticket/301/ql_qr_RICO.png\n\nMigrated from http://carson.math.uwm.edu/trac/clubb/ticket/301\n```json\n{\n    \"status\": \"closed\", \n    \"changetime\": \"2010-04-21T21:38:00\", \n    \"description\": \"'''Introduction:'''\\n\\nThe results of LH sampling seem to converge, but they are too noisy for cases with small cloud fraction, such as RICO.  To avoid noise, we need to sample preferentially from within cloud when cloud fraction is small.\\n\\n'''Requirements:'''\\n\\nGenerate two columns: one within cloud at k_s_mellor and one within clear air at k_s_mellor.  From these columns, compute a weighted average.  Eventually, we'll need to generalize for any even number of sample points per grid box and timestep.  However, we'll forego the freedom to use odd numbers of samples with this option.\\n\\n'''Technical spec:'''\\n\\n1.  Keep the present code intact.  Enclose all new code in a logical flag, l_lh_cloud_weighted_sampling.  For now, we could assume that we'll only use this with l_lh_vert_overlap = .true., although later we may possibly wish to relax this.  Just in case, let's add some code like the following at the beginning of the latin_hypercube_sampling subroutine:\\n{{{\\n    ! Sanity check for l_lh_cloud_weighted_sampling\\n    if ( l_lh_cloud_weighted_sampling .and. mod( n_micro_calls, 2 ) /= 0 ) then\\n      write(fstderr,*) \\\"Cloud weighted sampling requires micro calls to be divisible by 2.\\\"\\n      stop \\\"Fatal error.\\\"\\n    end if\\n}}}\\n\\n2.  At the k_lh_start, if cloud_frac > 0.5 .or. cloud_frac < 0.01, then choose points randomly, as we do now.  (However, don't leave 0.01 as a magic number; define it as a variable.)\\n\\n3.  On the other hand, if at k_lh_start, 0.01 < cloud_frac < 0.5, then keep choosing points until you find one that lies in the cloudy portion.  E.g., suppose you choose mixture fraction 1 (using the d+1 variate).  Then, if X_u(1,...) > (1-cloud_frac1), then the sample point is within cloud.  If instead you had chosen mixture fraction 2, then use cloud_frac2.  If the point is not within cloud, then choose other values of X_u(d+1) and X_u(1).  When you find a point within cloud, stop.  (This rejection method might be slow; if so, we could choose a faster method.)\\nInitially we're going to use the rejection method, so that we have something to compare against a more sophisticated weighted method.  We'll also want assertion tests to verify that the cloudy columns are actually cloudy and the clear columns are actually clear.\\n  \\n4.  Then use the rejection method to keep choosing points until you find one that lies in the clear portion.  \\nTo do this we can modify {{{ generate_uniform_sample }}} to use a new function called {{{choose_permuted_random}}} rather than using an inline formula, and this will allow us to generate just one element rather than n_micro_calls * d_variables elements.\\nNew code:\\n{{{\\n   ! Choose values of sample using permuted vector and random number generator\\n    do j = 1,n_micro_calls\\n      do k = 1,dp1\\n        X_u_one_lev(j,k) = choose_permuted_random( nt_repeat, p_matrix(j,k) )\\n      end do\\n    end do\\n\\n!        print*, 'p_matrix(:,1)= ', p_matrix(:,1)\\n!        print*, 'p_matrix(:,dp1)= ', p_matrix(:,dp1)\\n!        print*, 'X(:,1)= ', X(:,:)\\n\\n    return\\n  end subroutine generate_uniform_sample\\n\\n!----------------------------------------------------------------------\\n  double precision function choose_permuted_random( nt_repeat, p_matrix_element )\\n\\n    use mt95, only: genrand_real3 ! Procedure(s)\\n\\n    use mt95, only: genrand_real ! Constants\\n\\n    implicit none\\n\\n    integer, intent(in) :: & \\n      nt_repeat, &  ! Number of samples before the sequence repeats\\n      p_matrix_element       ! Permuted integer\\n\\n    real(kind=genrand_real) :: rand ! Random float with a range of (0,1)\\n\\n    ! ---- Begin Code ----\\n\\n    call genrand_real3( rand ) ! genrand_real3's range is (0,1)\\n\\n    choose_permuted_random = (1.0d0/nt_repeat)*(p_matrix_element + rand )\\n\\n    return\\n  end function choose_permuted_random\\n}}}\\nThen, when we have our X_u_all_levs, we can re-use {{{choose_permuted_random}}} on k_lh_start for step 5.\\n\\n5.  Once you have found the two sample points at k_lh_start, proceed as usual, choosing the full column of points.  \\n\\n6.  In the end, we need to weight the results of the two sample columns by cloud_frac and (1-cloud_frac).\\nThis will be tricky.  It looks like we can use the code in {{{l_cloud_weighted_averaging}}}, but in the past that code didn't work as expected.  We'll have to debug this I think.\\n\\n'''Tests:'''\\n\\n1.  Ensure that at any timestep, we always have one sample point within cloud and one out of cloud at level k_lh_start.\\n\\n2. Test whether the new LH sample converges to the analytic solution for KK microphysics for the RF02 and RICO cases.\\n\\n3. See whether the profiles vary less in time when (l_lh_cloud_weighted_sampling = .true.) than when (l_lh_cloud_weighted_sampling = .false.).  \\n\\n'''Time estimate:'''\\n\\nThe tricky part will probably be going from the rejection method to something more sophisticated.\\n\\n'''Special note'''\\n\\nWe should keep in mind that this undoes some of the Latin Hypercube sampling aspects of the code (i.e. we're overwriting s with an in-cloud value when t has already been determined using Latin Hypercube sampling).\", \n    \"reporter\": \"vlarson@uwm.edu\", \n    \"cc\": \"vlarson@uwm.edu\", \n    \"resolution\": \"Verified by V. Larson\", \n    \"_ts\": \"1271885880000000\", \n    \"component\": \"clubb_src\", \n    \"summary\": \"Latin hypercube: Sample preferentially within cloud\", \n    \"priority\": \"major\", \n    \"keywords\": \"\", \n    \"time\": \"2010-03-15T17:10:16\", \n    \"milestone\": \"Improve SILHS\", \n    \"owner\": \"dschanen@uwm.edu\", \n    \"type\": \"enhancement\"\n}\n```\n"},"comment":{"url":"https://api.github.com/repos/StrikeNP/trac_test/issues/comments/395800535","html_url":"https://github.com/StrikeNP/trac_test/issues/301#issuecomment-395800535","issue_url":"https://api.github.com/repos/StrikeNP/trac_test/issues/301","id":395800535,"node_id":"MDEyOklzc3VlQ29tbWVudDM5NTgwMDUzNQ==","user":{"login":"StrikeNP","id":2577059,"node_id":"MDQ6VXNlcjI1NzcwNTk=","avatar_url":"https://avatars0.githubusercontent.com/u/2577059?v=4","gravatar_id":"","url":"https://api.github.com/users/StrikeNP","html_url":"https://github.com/StrikeNP","followers_url":"https://api.github.com/users/StrikeNP/followers","following_url":"https://api.github.com/users/StrikeNP/following{/other_user}","gists_url":"https://api.github.com/users/StrikeNP/gists{/gist_id}","starred_url":"https://api.github.com/users/StrikeNP/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/StrikeNP/subscriptions","organizations_url":"https://api.github.com/users/StrikeNP/orgs","repos_url":"https://api.github.com/users/StrikeNP/repos","events_url":"https://api.github.com/users/StrikeNP/events{/privacy}","received_events_url":"https://api.github.com/users/StrikeNP/received_events","type":"User","site_admin":false},"created_at":"2018-06-08T15:38:52Z","updated_at":"2018-06-08T15:38:52Z","author_association":"OWNER","body":"Trac update at 2010/03/19 at 20:33\n\ndschanen@uwm.edu commented:\n\nReplying to [comment:5 dschanen@…]:\n> Looking back at this comment it doesn't make any of sense.  DYCOMS II RF02 does appear to converge.  What I meant to say is that the time averaged lh_rcm doesn't look like rcm, and time averaged lh_rrainm doesn't look like rrainm in the RICO case (and so the tendencies look different as well, but we would expect that for different sample points).  The RF02 simulation looks reasonable, but since we're in the > 50% range almost everywhere it's not a very good test case.\n\nWe determined that this had to do with the weights that were being applied to the individual terms in the sum.  If RICO is  simulated with the maximum overlap assumption for all variates and average over time, then we get results that look ok, but there is a lot of oscillation in lh_rcm at individual timestep when n_micro_calls = 2.  The next steps are:\n1. Check that variances like lh_rcp2_zt look correct with a large number of sample points.  This should cut down on the variability of lh_rcm from timestep to timestep too.\n2. Verify the results when only the s_mellor variable is assumed to be maximally overlapped.  In the past lh_rrainm looked noisy.\n3. Re-simulate RF02 and verify it still works as expected.\n4. Run RICO interactively with 2 sample points, and max correlation for all variates.\n5. Use plotgen to compare the simulation in 4 to an analytic solution.  Past simulations using random overlap had extremely noisy rain rate plots.\n6. Commit the working version.\n7. Re-factor to make the code more readable (e.g. add functions where we have if ... then .. else if now).\n8. Re-produce the results from before the code was re-factored.\n9. Re-commit the cleaned up and working code.\n\n"}}
{"actor":{"display_login":"StrikeNP"}}