From dda9a924e6522e0f21dd6f1052c8346dadf0dba1 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Wed, 7 Jun 2023 16:49:09 -0400 Subject: [PATCH] apply mixing volume changes to allocatable_arrays branch --- .gitignore | 2 ++ examples/global-lunar-bombardment/qmctest2.in | 8 +++++ examples/global-lunar-bombardment/qmctest3.in | 3 ++ src/regolith/regolith_mix.f90 | 29 ++++++++----------- src/util/util_traverse_pop_array.f90 | 22 ++++++++++---- 5 files changed, 41 insertions(+), 23 deletions(-) create mode 100644 examples/global-lunar-bombardment/qmctest2.in create mode 100644 examples/global-lunar-bombardment/qmctest3.in diff --git a/.gitignore b/.gitignore index bc510bb1..3b3bb35a 100644 --- a/.gitignore +++ b/.gitignore @@ -76,3 +76,5 @@ python/ctem/ctem.egg-info/ *.i90 INSTALL + +examples/global-lunar-bombardment/temp.in diff --git a/examples/global-lunar-bombardment/qmctest2.in b/examples/global-lunar-bombardment/qmctest2.in new file mode 100644 index 00000000..0c90086d --- /dev/null +++ b/examples/global-lunar-bombardment/qmctest2.in @@ -0,0 +1,8 @@ +#Dcrat(m) vel(m/s) ang(deg) xoffset(m) yoffset(m) t(Ga) +2050000 18300.0 45.0 -2908888.89 -1607137.57 0.6 +658628 18300.0 45.0 570078.99 770213.11 0.5 +1012000 18300.0 45.0 -297849.189 1121963.97 0.2 +623000 18300.0 45.0 1064349.59 -473044.26 0.4 +809000 18300.0 45.0 1670883.00 509433.00 0.3 +639000 18300.0 45.0 -1874653.62 -609499.34 0.1 +233931 18300.0 45.0 -533541.38 1523910.2 0.05 \ No newline at end of file diff --git a/examples/global-lunar-bombardment/qmctest3.in b/examples/global-lunar-bombardment/qmctest3.in new file mode 100644 index 00000000..a61f2281 --- /dev/null +++ b/examples/global-lunar-bombardment/qmctest3.in @@ -0,0 +1,3 @@ +#Dcrat(m) vel(m/s) ang(deg) xoffset(m) yoffset(m) t(Ga) +886977 18300.0 45.0 1482572.48 -62471.4873 0.02 +10000 18300.0 45.0 70000.0 10000.0 0.006 \ No newline at end of file diff --git a/src/regolith/regolith_mix.f90 b/src/regolith/regolith_mix.f90 index 775391a9..7dec2759 100644 --- a/src/regolith/regolith_mix.f90 +++ b/src/regolith/regolith_mix.f90 @@ -49,8 +49,8 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain) newlayer%age(:) = 0.0_DP allocate(newlayer%meltdist(1+domain%rcnum)) allocate(newlayer%distvol(1+domain%rcnum)) - newlayer%meltdist(:) = 0.0_SP - newlayer%distvol(:) = 0.0_SP + newlayer%meltdist(:) = 0.0_DP + newlayer%distvol(:) = 0.0_DP newlayer%ejm = 0.0_DP newlayer%ejmf = 0.0_DP newlayer%meltvolume = 0.0_DP @@ -62,27 +62,15 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain) do i = N,1,-1 newlayer%thickness = newlayer%thickness + poppedarray(i)%thickness newlayer%comp = newlayer%comp + poppedarray(i)%thickness * poppedarray(i)%comp - !newlayer%meltfrac = newlayer%meltfrac + poppedarray(i)%thickness * poppedarray(i)%meltfrac newlayer%age(:) = newlayer%age(:) + poppedarray(i)%age(:) - !newlayer%meltdist(:) = newlayer%meltdist(:) + poppedarray(i)%thickness * poppedarray(i)%meltdist(:) - newlayer%distvol(:) = newlayer%distvol(:) + poppedarray(i)%thickness * poppedarray(i)%distvol(:) - newlayer%ejm = newlayer%ejm + poppedarray(i)%thickness * poppedarray(i)%ejm - !newlayer%ejmf = newlayer%ejmf + poppedarray(i)%thickness * poppedarray(i)%ejmf - newlayer%meltvolume = newlayer%meltvolume + poppedarray(i)%thickness * poppedarray(i)%meltvolume - !newlayer%totvolume = newlayer%totvolume + poppedarray(i)%thickness * poppedarray(i)%totvolume + newlayer%distvol(:) = newlayer%distvol(:) + poppedarray(i)%distvol(:) + newlayer%ejm = newlayer%ejm + poppedarray(i)%ejm + newlayer%meltvolume = newlayer%meltvolume + poppedarray(i)%meltvolume end do ! Get average values of composition and melt fraction newlayer%comp = newlayer%comp / newlayer%thickness - !newlayer%meltfrac = newlayer%meltfrac / newlayer%thickness - !newlayer%meltdist(:) = newlayer%meltdist(:) / newlayer%thickness - newlayer%distvol(:) = newlayer%distvol(:) / newlayer%thickness - newlayer%ejm = newlayer%ejm / newlayer%thickness - !newlayer%ejmf = newlayer%ejmf / newlayer%thickness - newlayer%meltvolume = newlayer%meltvolume / newlayer%thickness - !newlayer%totvolume = newlayer%totvolume / newlayer%thickness - !test forcing melt fraction to equal melt volume / total volume newlayer%totvolume = newlayer%thickness * user%pix * user%pix newlayer%meltfrac = newlayer%meltvolume / newlayer%totvolume newlayer%meltdist(:) = newlayer%distvol(:) / newlayer%totvolume @@ -91,5 +79,12 @@ subroutine regolith_mix(user,surfi,mixing_depth,domain) call util_push_array(surfi%regolayer, newlayer) !call util_destroy_list(poppedlist_top) + + ! do i = N,1,-1 + ! if (abs(surfi%regolayer(i)%meltvolume - sum(surfi%regolayer(i)%distvol) > 1e-5)) then + ! write(*,*) "melt array =/= melt value!", domain%nqmc, domain%currentqmc, abs(surfi%regolayer(i)%meltvolume - sum(surfi%regolayer(i)%distvol)) + ! end if + ! end do + return end subroutine regolith_mix diff --git a/src/util/util_traverse_pop_array.f90 b/src/util/util_traverse_pop_array.f90 index 3ba03cbc..fc670800 100644 --- a/src/util/util_traverse_pop_array.f90 +++ b/src/util/util_traverse_pop_array.f90 @@ -65,19 +65,29 @@ subroutine util_traverse_pop_array(user,regolayer,traverse_depth,poppedarray) allocate(poppedarray,source=regolayer(maxi:N)) + allocate(oldregodata,source=regolayer(maxi:maxi)) !for #1 element of poppedarray, shrink thickness by whatever was lefr over. In corresponding maxi of regolayer, also need to change that. poppedarray(1)%thickness = poppedarray(1)%thickness - depth_diff regolayer(maxi)%thickness = regolayer(maxi)%thickness - poppedarray(1)%thickness + poppedarray(1)%meltvolume = (poppedarray(1)%thickness/oldregodata(1)%thickness) * oldregodata(1)%meltvolume + poppedarray(1)%distvol(:) = (poppedarray(1)%thickness/oldregodata(1)%thickness) * oldregodata(1)%distvol(:) + poppedarray(1)%ejm = (poppedarray(1)%thickness/oldregodata(1)%thickness) * oldregodata(1)%ejm + regolayer(maxi)%meltvolume = (regolayer(maxi)%thickness/oldregodata(1)%thickness) * oldregodata(1)%meltvolume + regolayer(maxi)%distvol(:) = (regolayer(maxi)%thickness/oldregodata(1)%thickness) * oldregodata(1)%distvol(:) + regolayer(maxi)%ejm = (regolayer(maxi)%thickness/oldregodata(1)%thickness) * oldregodata(1)%ejm + poppedarray(1)%totvolume = poppedarray(1)%thickness * user%pix * user%pix - poppedarray(1)%meltvolume = poppedarray(1)%meltfrac * poppedarray(1)%totvolume - poppedarray(1)%distvol(:) = poppedarray(1)%meltdist(:) * poppedarray(1)%totvolume - poppedarray(1)%ejm = poppedarray(1)%ejmf * poppedarray(1)%totvolume + if(poppedarray(1)%totvolume > 0) then + poppedarray(1)%meltfrac = poppedarray(1)%meltvolume / poppedarray(1)%totvolume + poppedarray(1)%ejmf = poppedarray(1)%ejm / poppedarray(1)%totvolume + end if regolayer(maxi)%totvolume = regolayer(maxi)%thickness * user%pix * user%pix - regolayer(maxi)%distvol(:) = regolayer(maxi)%meltdist(:) * regolayer(maxi)%totvolume - regolayer(maxi)%meltvolume = regolayer(maxi)%meltfrac * regolayer(maxi)%totvolume - regolayer(maxi)%ejm = regolayer(maxi)%ejmf * regolayer(maxi)%totvolume + regolayer(maxi)%meltfrac = regolayer(maxi)%meltvolume / regolayer(maxi)%totvolume + regolayer(maxi)%ejmf = regolayer(maxi)%ejm / regolayer(maxi)%totvolume + + deallocate(oldregodata) ! copy regolayer from 1 to maxi to temp variable, then deallocate regolayer, then movealloc templayer onto regolayer <--may need temp array allocate(oldregodata,source=regolayer(1:maxi))