diff --git a/Project.toml b/Project.toml index 55b6a74..d41276a 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Jan Swierczek-Jereczek"] version = "0.1.0" [deps] +AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -16,7 +17,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] CUDA = "^4.0" @@ -30,7 +30,6 @@ JLD2 = "^0.4" OrdinaryDiffEq = "^6" Reexport = "^1.2" Statistics = "^1.9" -StatsBase = "~0.33" julia = "^1.9" [extras] diff --git a/docs/build/APIref.html b/docs/build/APIref.html index 4c82551..9dd8b90 100644 --- a/docs/build/APIref.html +++ b/docs/build/APIref.html @@ -1,2 +1,2 @@ -API reference · FastIsostasy.jl
Edit on GitHub

API reference

Basic structs

FastIsostasy.MultilayerEarth — Type
MultilayerEarth

Return a struct containing all information related to the radially layered structure of the solid Earth and its parameters.

source
Missing docstring.

Missing docstring for RefGeoState. Check Documenter's build log for details.

Missing docstring.

Missing docstring for FastisoResults. Check Documenter's build log for details.

Mechanics

Missing docstring.

Missing docstring for forward_isostasy. Check Documenter's build log for details.

Missing docstring.

Missing docstring for init_results. Check Documenter's build log for details.

Missing docstring.

Missing docstring for forwardstep_isostasy!. Check Documenter's build log for details.

Missing docstring.

Missing docstring for simple_euler!. Check Documenter's build log for details.

Missing docstring.

Missing docstring for apply_bc. Check Documenter's build log for details.

Missing docstring.

Missing docstring for compute_elastic_response. Check Documenter's build log for details.

Sea-level

FastIsostasy.update_geostate! — Function
update_geostate!(sstruct::SuperStruct, u::Matrix, H_ice::Matrix)

Update the ::GeoState computing the current geoid perturbation, the sea-level changes and the load columns for the next time step of the isostasy integration.

source
FastIsostasy.update_geoid! — Function
update_geoid!(sstruct::SuperStruct)

Update the geoid of a ::GeoState by convoluting the Green's function with the load change.

source
FastIsostasy.get_geoidgreen — Function
get_geoidgreen(sstruct::SuperStruct)

Return the Green's function used to compute the changes in geoid.

Reference

Coulon et al. 2021.

source
FastIsostasy.update_sealevel! — Function
update_sealevel!(sstruct::SuperStruct)

Update the sea-level ::GeoState by adding the various contributions.

Reference

Coulon et al. (2021), Figure 1.

source
FastIsostasy.update_slc! — Function
update_slc!(sstruct::SuperStruct)

Update the sea-level contribution of melting above floatation, density correction and potential ocean volume.

Reference

Goelzer et al. (2020), eq. (12).

source
FastIsostasy.update_V_af! — Function
update_V_af!(sstruct::SuperStruct)

Update the ice volume above floatation.

Reference

Goelzer et al. (2020), eq. (13). Note: we do not use eq. (1) as it is only a special case of eq. (13) that does not allow a correct representation of external sea-level forcings.

source
FastIsostasy.update_slc_af! — Function
update_slc_af!(sstruct::SuperStruct)

Update the sea-level contribution of ice above floatation.

Reference

Goelzer et al. (2020), eq. (2).

source
FastIsostasy.update_V_pov! — Function
update_V_pov!(sstruct::SuperStruct)

Update the potential ocean volume.

Reference

Goelzer et al. (2020), eq. (14). Note: we do not use eq. (8) as it is only a special case of eq. (14) that does not allow a correct representation of external sea-level forcinsstruct.geostate.

source
FastIsostasy.update_slc_pov! — Function
update_slc_pov!(sstruct::SuperStruct)

Update the sea-level contribution associated with the potential ocean volume.

Reference

Goelzer et al. (2020), eq. (9).

source
FastIsostasy.update_V_den! — Function
update_V_den!(sstruct::SuperStruct)

Update the ocean volume associated with the density correction.

Reference

Goelzer et al. (2020), eq. (10).

source
FastIsostasy.update_slc_den! — Function
update_slc_den!(sstruct::SuperStruct)

Update the sea-level contribution associated with the density correction.

Reference

Goelzer et al. (2020), eq. (11).

source

Parameter inversion

Utils

Missing docstring.

Missing docstring for matrify_constant. Check Documenter's build log for details.

Missing docstring.

Missing docstring for sphericaldistance. Check Documenter's build log for details.

Missing docstring.

Missing docstring for scalefactor. Check Documenter's build log for details.

FastIsostasy.latlon2stereo — Function
latlon2stereo(lat, lon, lat0, lon0)

Compute stereographic projection (x,y) for a given latitude lat longitude lon, reference latitude lat0 and reference longitude lon0. Optionally one can provide lat::AbstractMatrix and lon::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be provided in degrees! Reference: John P. Snyder (1987), p. 157, eq. (21-2), (21-3), (21-4).

source
FastIsostasy.stereo2latlon — Function
stereo2latlon(x, y, lat0, lon0)

Compute the inverse stereographic projection (lat, lon) based on Cartesian coordinates (x,y) and for a given reference latitude lat0 and reference longitude lon0. Optionally one can provide x::AbstractMatrix and y::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be para elloprovided in degrees!

Convert stereographic (x,y)-coordinates to latitude-longitude. Reference: John P. Snyder (1987), p. 159, eq. (20-14), (20-15), (20-18), (21-15).

source
FastIsostasy.get_rigidity — Function
get_rigidity(t::T, E::T, nu::T) where {T<:AbstractFloat}

Compute rigidity D based on thickness t, Young modulus E and Poisson ration nu.

source
Missing docstring.

Missing docstring for get_effective_viscosity. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_differential_fourier. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_viscosity_ratio. Check Documenter's build log for details.

Missing docstring.

Missing docstring for three_layer_scaling. Check Documenter's build log for details.

FastIsostasy.loginterp_viscosity — Function
loginterp_viscosity(tvec, layer_viscosities, layers_thickness, pseudodiff)

Compute a log-interpolator of the equivalent viscosity from provided viscosity fields layer_viscosities at time stamps tvec.

source
Missing docstring.

Missing docstring for hyperbolic_channel_coeffs. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_greenintegrand_coeffs. Check Documenter's build log for details.

Missing docstring.

Missing docstring for build_greenintegrand. Check Documenter's build log for details.

FastIsostasy.get_elasticgreen — Function
get_elasticgreen(Omega, quad_support, quad_coeffs)

Integrate load response over field by using 2D quadrature with specified support points and associated coefficients.

source
FastIsostasy.quadrature1D — Function
quadrature1D(f, n, x1, x2)

Compute 1D Gauss-Legendre quadrature of f between x1 and x2 based on n support points.

source
Missing docstring.

Missing docstring for quadrature2D. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_normalized_lin_transform. Check Documenter's build log for details.

Missing docstring.

Missing docstring for normalized_lin_transform. Check Documenter's build log for details.

+API reference · FastIsostasy.jl
Edit on GitHub

API reference

Basic structs

FastIsostasy.MultilayerEarth — Type
MultilayerEarth

Return a struct containing all information related to the radially layered structure of the solid Earth and its parameters.

source
Missing docstring.

Missing docstring for RefGeoState. Check Documenter's build log for details.

Missing docstring.

Missing docstring for FastisoResults. Check Documenter's build log for details.

Mechanics

Missing docstring.

Missing docstring for forward_isostasy. Check Documenter's build log for details.

Missing docstring.

Missing docstring for init_results. Check Documenter's build log for details.

Missing docstring.

Missing docstring for forwardstep_isostasy!. Check Documenter's build log for details.

Missing docstring.

Missing docstring for simple_euler!. Check Documenter's build log for details.

Missing docstring.

Missing docstring for apply_bc. Check Documenter's build log for details.

Missing docstring.

Missing docstring for compute_elastic_response. Check Documenter's build log for details.

Sea-level

FastIsostasy.update_geostate! — Function
update_geostate!(sstruct::SuperStruct, u::Matrix, H_ice::Matrix)

Update the ::GeoState computing the current geoid perturbation, the sea-level changes and the load columns for the next time step of the isostasy integration.

source
FastIsostasy.update_geoid! — Function
update_geoid!(sstruct::SuperStruct)

Update the geoid of a ::GeoState by convoluting the Green's function with the load change.

source
FastIsostasy.get_geoidgreen — Function
get_geoidgreen(sstruct::SuperStruct)

Return the Green's function used to compute the changes in geoid.

Reference

Coulon et al. 2021.

source
FastIsostasy.update_sealevel! — Function
update_sealevel!(sstruct::SuperStruct)

Update the sea-level ::GeoState by adding the various contributions.

Reference

Coulon et al. (2021), Figure 1.

source
FastIsostasy.update_slc! — Function
update_slc!(sstruct::SuperStruct)

Update the sea-level contribution of melting above floatation, density correction and potential ocean volume.

Reference

Goelzer et al. (2020), eq. (12).

source
FastIsostasy.update_V_af! — Function
update_V_af!(sstruct::SuperStruct)

Update the ice volume above floatation.

Reference

Goelzer et al. (2020), eq. (13). Note: we do not use eq. (1) as it is only a special case of eq. (13) that does not allow a correct representation of external sea-level forcings.

source
FastIsostasy.update_slc_af! — Function
update_slc_af!(sstruct::SuperStruct)

Update the sea-level contribution of ice above floatation.

Reference

Goelzer et al. (2020), eq. (2).

source
FastIsostasy.update_V_pov! — Function
update_V_pov!(sstruct::SuperStruct)

Update the potential ocean volume.

Reference

Goelzer et al. (2020), eq. (14). Note: we do not use eq. (8) as it is only a special case of eq. (14) that does not allow a correct representation of external sea-level forcinsstruct.geostate.

source
FastIsostasy.update_slc_pov! — Function
update_slc_pov!(sstruct::SuperStruct)

Update the sea-level contribution associated with the potential ocean volume.

Reference

Goelzer et al. (2020), eq. (9).

source
FastIsostasy.update_V_den! — Function
update_V_den!(sstruct::SuperStruct)

Update the ocean volume associated with the density correction.

Reference

Goelzer et al. (2020), eq. (10).

source
FastIsostasy.update_slc_den! — Function
update_slc_den!(sstruct::SuperStruct)

Update the sea-level contribution associated with the density correction.

Reference

Goelzer et al. (2020), eq. (11).

source

Parameter inversion

Utils

Missing docstring.

Missing docstring for matrify. Check Documenter's build log for details.

Missing docstring.

Missing docstring for sphericaldistance. Check Documenter's build log for details.

Missing docstring.

Missing docstring for scalefactor. Check Documenter's build log for details.

FastIsostasy.latlon2stereo — Function
latlon2stereo(lat, lon, lat0, lon0)

Compute stereographic projection (x,y) for a given latitude lat longitude lon, reference latitude lat0 and reference longitude lon0. Optionally one can provide lat::AbstractMatrix and lon::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be provided in degrees! Reference: John P. Snyder (1987), p. 157, eq. (21-2), (21-3), (21-4).

source
FastIsostasy.stereo2latlon — Function
stereo2latlon(x, y, lat0, lon0)

Compute the inverse stereographic projection (lat, lon) based on Cartesian coordinates (x,y) and for a given reference latitude lat0 and reference longitude lon0. Optionally one can provide x::AbstractMatrix and y::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be para elloprovided in degrees!

Convert stereographic (x,y)-coordinates to latitude-longitude. Reference: John P. Snyder (1987), p. 159, eq. (20-14), (20-15), (20-18), (21-15).

source
FastIsostasy.get_rigidity — Function
get_rigidity(t::T, E::T, nu::T) where {T<:AbstractFloat}

Compute rigidity D based on thickness t, Young modulus E and Poisson ration nu.

source
Missing docstring.

Missing docstring for get_effective_viscosity. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_differential_fourier. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_viscosity_ratio. Check Documenter's build log for details.

Missing docstring.

Missing docstring for three_layer_scaling. Check Documenter's build log for details.

FastIsostasy.loginterp_viscosity — Function
loginterp_viscosity(tvec, layer_viscosities, layers_thickness, pseudodiff)

Compute a log-interpolator of the equivalent viscosity from provided viscosity fields layer_viscosities at time stamps tvec.

source
Missing docstring.

Missing docstring for hyperbolic_channel_coeffs. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_greenintegrand_coeffs. Check Documenter's build log for details.

Missing docstring.

Missing docstring for build_greenintegrand. Check Documenter's build log for details.

FastIsostasy.get_elasticgreen — Function
get_elasticgreen(Omega, quad_support, quad_coeffs)

Integrate load response over field by using 2D quadrature with specified support points and associated coefficients.

source
FastIsostasy.quadrature1D — Function
quadrature1D(f, n, x1, x2)

Compute 1D Gauss-Legendre quadrature of f between x1 and x2 based on n support points.

source
Missing docstring.

Missing docstring for quadrature2D. Check Documenter's build log for details.

Missing docstring.

Missing docstring for get_normalized_lin_transform. Check Documenter's build log for details.

Missing docstring.

Missing docstring for normalized_lin_transform. Check Documenter's build log for details.

diff --git a/docs/build/search_index.js b/docs/build/search_index.js index e46a3fd..2aeb073 100644 --- a/docs/build/search_index.js +++ b/docs/build/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"examples.html#Examples","page":"Examples","title":"Examples","text":"","category":"section"},{"location":"examples.html#Multi-layer-Earth","page":"Examples","title":"Multi-layer Earth","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"FastIsostasy relies on a (polar) stereographic projection that allows to treat the radially-layered, onion-like structure of the solid Earth as a superposition of horizontal layers. Furthermore, FastIsostasy reduces this 3D problem into a 2D problem by collapsing the depth dimension by computing an effective viscosity field. The user is required to provide the 3D information, which will then be used under the hood to compute the effective viscosity. This tutorial shows such an example.","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"We want to render a situation similar to the one depicted below:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"(Image: Schematic representation of the three-layer set-up.)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Initializing a MultilayerEarth with parameters corresponding to this situation automatically computes the conversion from a 3D to a 2D problem. This can be simply executed by running:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = 3000e3 # (m) half-width of the domain\nn = 8 # implies an NxN grid with N = 2^n = 256.\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()\n\nlv = [1e19, 1e21] # (Pa*s)\nlb = [88e3, 400e3] # (m)\np = MultilayerEarth(Omega, c, layer_viscosities = lv, layer_boundaries = lb)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"The next section shows how to use the now obtained p::MultilayerEarth for actual GIA computation.","category":"page"},{"location":"examples.html#Simple-load-and-geometry","page":"Examples","title":"Simple load and geometry","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"We now apply a constant load, here a cylinder of ice with radius R = 1000 km and thickness H = 1 km, over the domain introduced in Multi-layer Earth. To obtain the bedrock displacement over time and store it at time steps specified by a vector t_out, we can use the convenience function fastisostasy:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"R = 1000e3 # ice disc radius (m)\nH = 1000.0 # ice disc thickness (m)\nHice = uniform_ice_cylinder(Omega, R, H)\nt_out = years2seconds.([0.0, 100.0, 500.0, 1500.0, 5000.0, 10_000.0, 50_000.0])\n\nresults = fastisostasy(t_out, Omega, c, p, Hice)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Yes, that was it! You can now easily access the elastic and viscous displacement by calling results.elastic or results.viscous. For the present case, the latter can be compared to an analytic solution that is known for this particular case. Let's look at the accuracy of our numerical scheme over time by running following plotting commands:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"M = Omega.N ÷ 2\nx = diag(Omega.X)[1:M]\ny = diag(Omega.Y)[1:M]\nr = sqrt.( x .^ 2 + y .^ 2 )\n\nfig = Figure()\nax3 = Axis(fig[1, 1])\ncolors = [:gray80, :gray65, :gray50, :gray35, :gray20, :gray5]\n\nfor i in eachindex(t_out)\n t = t_out[i]\n analytic_solution_r(r) = analytic_solution(r, t, c, p, H, R, analytic_support)\n u_analytic = analytic_solution_r.( r )\n u_numeric = diag(results.viscous[i])\n lines!(ax3, x, u_numeric[1:M], color = colors[i], linewidth = 5)\n lines!(ax3, x, u_analytic, color = colors[i], linewidth = 5, linestyle = :dash)\nend\nfig","category":"page"},{"location":"examples.html#Time-changing-load","page":"Examples","title":"Time-changing load","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"That looks pretty good! One might however object that the convenience function fastisostasy ends up being not so convenient as soon as the ice load changes over time. This case can however be easily handled by providing snapshots of the ice thickness and their assoicated time. By passing this to fastisostasy, an interpolator is created and called within the time integration. Let's create a tool example where the thickness of the ice cylinder asymptotically grows from 0 to 1 km, this can be implemented by:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"normalized_asymptote(t) = 1 - exp(-t)\nt_Hice_asymptotic = 0:10.0:t_out[end]\nHice_asymptotic = [normalized_asymptote(t) .* Hice for t in t_Hice_asymptote]\nresults = fastisostasy(t_out, Omega, c, p, t_Hice_asymptotic, Hice_asymptotic)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"This concept will also apply to the upper-mantle viscosity in future versions, as it can change over time.","category":"page"},{"location":"examples.html#GPU-support","page":"Examples","title":"GPU support","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"For about n 6, the previous example can be computed even faster by using GPU parallelism. It could not represent less work from the user's perspective, as it boils down to calling the ComputationDomain with an extra keyword argument:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Omega = ComputationDomain(W, n, use_cuda=true)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"That's it, nothing more! One could suggest you lay back but your computation might be completed too soon for that.","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"info: Only CUDA supported!\nFor now only Nvidia GPUs are supported (sorry Mac users, destiny is taking its revenge upon you for being so cool) and there is no plan of extending this compatibility at this point.","category":"page"},{"location":"examples.html#Simple-load-and-geometry-DIY","page":"Examples","title":"Simple load and geometry - DIY","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Nonetheless, as any high-level convenience function, fastisostasy has limitations. An ice-sheet modeller typically wants to embed FastIsostasy within a time-stepping loop. This can be easily done by getting familiar with some intermediate-level functions:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = 3000e3 # (m) half-width of the domain\nn = 8 # implies an NxN grid with N = 2^n = 256.\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()\np = MultilayerEarth(Omega, c)\n\nR = 1000e3 # ice disc radius (m)\nH = 1000.0 # ice disc thickness (m)\nHice = uniform_ice_cylinder(Omega, R, H)\nt_out = years2seconds.([0.0, 100.0, 500.0, 1500.0, 5000.0, 10_000.0, 50_000.0])\n\nresults = fastisostasy(t_out, Omega, c, p, Hice)","category":"page"},{"location":"examples.html#GIA-following-Antarctic-deglaciation","page":"Examples","title":"GIA following Antarctic deglaciation","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"We now want to provide a tough example that presents:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"a heterogeneous lithosphere thickness\na heterogeneous upper-mantle viscosity\nvarious viscous channels\na more elaborate load that evolves over time\nchanges in the sea-level","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"For this we run a deglaciation of Antarctica, based on the ice thickness estimated in GLAC1D.","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = 3000e3 # (m) half-width of the domain\nn = 8 # implies an NxN grid with N = 2^n = 256.\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()","category":"page"},{"location":"examples.html#Inversion-of-solid-Earth-parameters","page":"Examples","title":"Inversion of solid-Earth parameters","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"FastIsostasy.jl relies on simplification of the full problem and might therefore need a calibration step to match the output of a 3D GIA model. By means of an unscented Kalman inversion, one can e.g. infer the appropriate effective upper-mantle viscosity based on the response of a 3D GIA model to a given load. Whereas this is know to be a tedious step, FastIsostasy is developped to ease the procedure by providing a convenience struct Paraminversion that can be run by:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = T(3000e3) # half-length of the square domain (m)\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()\n\nlb = [88e3, 180e3, 280e3, 400e3]\nlv = get_wiens_layervisc(Omega)\np = MultilayerEarth(\n Omega,\n c,\n layer_boundaries = lb,\n layer_viscosities = lv,\n)\nground_truth = copy(p.effective_viscosity)\n\nR = T(2000e3) # ice disc radius (m)\nH = T(1000) # ice disc thickness (m)\nHcylinder = uniform_ice_cylinder(Omega, R, H)\nt_out = years2seconds.(0.0:1_000.0:2_000.0)\n\nt1 = time()\nresults = fastisostasy(t_out, Omega, c, p, Hcylinder, ODEsolver=BS3(), interactive_geostate=false)\nt_fastiso = time() - t1\nprintln(\"Took $t_fastiso seconds!\")\nprintln(\"-------------------------------------\")\n\ntinv = t_out[2:end]\nHice = [Hcylinder for t in tinv]\nY = results.viscous[2:end]\nparaminv = ParamInversion(Omega, c, p, tinv, Y, Hice)\npriors, ukiobj = perform(paraminv)\nlogeta, Gx, e_mean, e_sort = extract_inversion(priors, ukiobj, paraminv)","category":"page"},{"location":"collected_garbage.html","page":"-","title":"-","text":"","category":"page"},{"location":"collected_garbage.html","page":"-","title":"-","text":"","category":"page"},{"location":"collected_garbage.html#Arrays","page":"-","title":"Arrays","text":"","category":"section"},{"location":"collected_garbage.html","page":"-","title":"-","text":"Array dimensions correspond to the spatial dimension of the variable they describe. If they evolve over time, they are stored as vector of arrays. For instance, the vertical displacement of the bedrock is a 2D variable that evolves over time. Therefore, it is stored in a Vector{Matrix}.","category":"page"},{"location":"APIref.html#API-reference","page":"API reference","title":"API reference","text":"","category":"section"},{"location":"APIref.html#Basic-structs","page":"API reference","title":"Basic structs","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"ComputationDomain\nPhysicalConstants\nMultilayerEarth\nRefSealevelState\nSealevelState\nPrecomputedFastiso\nSuperStruct\nFastisoResults","category":"page"},{"location":"APIref.html#FastIsostasy.ComputationDomain","page":"API reference","title":"FastIsostasy.ComputationDomain","text":"ComputationDomain\n\nReturn a struct containing all information related to geometry of the domain and potentially used parallelism.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.PhysicalConstants","page":"API reference","title":"FastIsostasy.PhysicalConstants","text":"PhysicalConstants\n\nReturn a struct containing important physical constants.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.MultilayerEarth","page":"API reference","title":"FastIsostasy.MultilayerEarth","text":"MultilayerEarth\n\nReturn a struct containing all information related to the radially layered structure of the solid Earth and its parameters.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.GeoState","page":"API reference","title":"FastIsostasy.GeoState","text":"GeoState\n\nReturn a mutable struct containing the geostate which will be updated over the simulation.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.PrecomputedFastiso","page":"API reference","title":"FastIsostasy.PrecomputedFastiso","text":"PrecomputedFastiso\n\nReturn a struct containing all the pre-computed terms needed for the forward integration of the model.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.SuperStruct","page":"API reference","title":"FastIsostasy.SuperStruct","text":"SuperStruct\n\nReturn a struct containing all the other structs needed for the forward integration of the model.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#Mechanics","page":"API reference","title":"Mechanics","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"ice_load\nfastisostasy\ninit_superstruct\nforward_isostasy\ninit_results\nforwardstep_isostasy!\ndudt_isostasy!\nsimple_euler!\napply_bc\ncompute_elastic_response","category":"page"},{"location":"APIref.html#FastIsostasy.ice_load","page":"API reference","title":"FastIsostasy.ice_load","text":"ice_load(c::PhysicalConstants{T}, H::T)\n\nCompute ice load based on ice thickness.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.fastisostasy","page":"API reference","title":"FastIsostasy.fastisostasy","text":"fastisostasy()\n\nMain function. List of all available solvers here.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.init_superstruct","page":"API reference","title":"FastIsostasy.init_superstruct","text":"init_superstruct()\n\nInit a SuperStruct containing all necessary fields for forward integration in time.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.dudt_isostasy!","page":"API reference","title":"FastIsostasy.dudt_isostasy!","text":"dudt_isostasy!()\n\nUpdate the displacement rate dudt.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#Sea-level","page":"API reference","title":"Sea-level","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"update_geostate!\nupdate_geoid!\nget_loadchange\nget_geoidgreen\nupdate_loadcolumns!\nupdate_sealevel!\nupdate_slc!\nupdate_V_af!\nupdate_slc_af!\nupdate_V_pov!\nupdate_slc_pov!\nupdate_V_den!\nupdate_slc_den!","category":"page"},{"location":"APIref.html#FastIsostasy.update_geostate!","page":"API reference","title":"FastIsostasy.update_geostate!","text":"update_geostate!(sstruct::SuperStruct, u::Matrix, H_ice::Matrix)\n\nUpdate the ::GeoState computing the current geoid perturbation, the sea-level changes and the load columns for the next time step of the isostasy integration.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_geoid!","page":"API reference","title":"FastIsostasy.update_geoid!","text":"update_geoid!(sstruct::SuperStruct)\n\nUpdate the geoid of a ::GeoState by convoluting the Green's function with the load change.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_loadchange","page":"API reference","title":"FastIsostasy.get_loadchange","text":"get_loadchange(sstruct::SuperStruct)\n\nCompute the load change compared to the reference configuration after Coulon et al. 2021. Difference to Coulon: HERE WE SCALE WITH CORRECTION FACTOR FOR PROJECTION!\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_geoidgreen","page":"API reference","title":"FastIsostasy.get_geoidgreen","text":"get_geoidgreen(sstruct::SuperStruct)\n\nReturn the Green's function used to compute the changes in geoid.\n\nReference\n\nCoulon et al. 2021.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_loadcolumns!","page":"API reference","title":"FastIsostasy.update_loadcolumns!","text":"update_loadcolumns!(sstruct::SuperStruct, u::XMatrix, H_ice::XMatrix)\n\nUpdate the load columns of a ::GeoState.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_sealevel!","page":"API reference","title":"FastIsostasy.update_sealevel!","text":"update_sealevel!(sstruct::SuperStruct)\n\nUpdate the sea-level ::GeoState by adding the various contributions.\n\nReference\n\nCoulon et al. (2021), Figure 1.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc!","page":"API reference","title":"FastIsostasy.update_slc!","text":"update_slc!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution of melting above floatation, density correction and potential ocean volume.\n\nReference\n\nGoelzer et al. (2020), eq. (12).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_V_af!","page":"API reference","title":"FastIsostasy.update_V_af!","text":"update_V_af!(sstruct::SuperStruct)\n\nUpdate the ice volume above floatation.\n\nReference\n\nGoelzer et al. (2020), eq. (13). Note: we do not use eq. (1) as it is only a special case of eq. (13) that does not allow a correct representation of external sea-level forcings.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc_af!","page":"API reference","title":"FastIsostasy.update_slc_af!","text":"update_slc_af!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution of ice above floatation.\n\nReference\n\nGoelzer et al. (2020), eq. (2).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_V_pov!","page":"API reference","title":"FastIsostasy.update_V_pov!","text":"update_V_pov!(sstruct::SuperStruct)\n\nUpdate the potential ocean volume.\n\nReference\n\nGoelzer et al. (2020), eq. (14). Note: we do not use eq. (8) as it is only a special case of eq. (14) that does not allow a correct representation of external sea-level forcinsstruct.geostate.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc_pov!","page":"API reference","title":"FastIsostasy.update_slc_pov!","text":"update_slc_pov!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution associated with the potential ocean volume.\n\nReference\n\nGoelzer et al. (2020), eq. (9).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_V_den!","page":"API reference","title":"FastIsostasy.update_V_den!","text":"update_V_den!(sstruct::SuperStruct)\n\nUpdate the ocean volume associated with the density correction.\n\nReference\n\nGoelzer et al. (2020), eq. (10).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc_den!","page":"API reference","title":"FastIsostasy.update_slc_den!","text":"update_slc_den!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution associated with the density correction.\n\nReference\n\nGoelzer et al. (2020), eq. (11).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#Parameter-inversion","page":"API reference","title":"Parameter inversion","text":"","category":"section"},{"location":"APIref.html#Utils","page":"API reference","title":"Utils","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"years2seconds\nseconds2years\nm_per_sec2mm_per_yr\nmatrify_vectorconstant\nmatrify_constant\nget_r\nmeshgrid\ndist2angulardist\nsphericaldistance\nscalefactor\nlatlon2stereo\nstereo2latlon\nget_rigidity\nget_effective_viscosity\nget_differential_fourier\nget_viscosity_ratio\nthree_layer_scaling\nloginterp_viscosity\nhyperbolic_channel_coeffs\nget_greenintegrand_coeffs\nbuild_greenintegrand\nget_elasticgreen\nget_quad_coeffs\nquadrature1D\nquadrature2D\nget_normalized_lin_transform\nnormalized_lin_transform\nkernelpromote","category":"page"},{"location":"APIref.html#FastIsostasy.years2seconds","page":"API reference","title":"FastIsostasy.years2seconds","text":"years2seconds(t::Real)\n\nConvert input time t from years to seconds.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.seconds2years","page":"API reference","title":"FastIsostasy.seconds2years","text":"seconds2years(t::Real)\n\nConvert input time t from seconds to years.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.m_per_sec2mm_per_yr","page":"API reference","title":"FastIsostasy.m_per_sec2mm_per_yr","text":"m_per_sec2mm_per_yr(dudt::Real)\n\nConvert displacement rate dudt from $ m \\, s^{-1} $ to $ mm \\, \\mathrm{yr}^{-1} $.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.matrify_vectorconstant","page":"API reference","title":"FastIsostasy.matrify_vectorconstant","text":"matrify_vectorconstant(x, N)\n\nGenerate a vector of constant matrices from a vector of constants.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_r","page":"API reference","title":"FastIsostasy.get_r","text":"get_r(x::T, y::T) where {T<:Real}\n\nGet euclidean distance of point (x, y) to origin.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.meshgrid","page":"API reference","title":"FastIsostasy.meshgrid","text":"meshgrid(x, y)\n\nReturn a 2D meshgrid spanned by x, y.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.dist2angulardist","page":"API reference","title":"FastIsostasy.dist2angulardist","text":"dist2angulardist(r::Real)\n\nConvert Euclidean to angular distance along great circle.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.latlon2stereo","page":"API reference","title":"FastIsostasy.latlon2stereo","text":"latlon2stereo(lat, lon, lat0, lon0)\n\nCompute stereographic projection (x,y) for a given latitude lat longitude lon, reference latitude lat0 and reference longitude lon0. Optionally one can provide lat::AbstractMatrix and lon::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be provided in degrees! Reference: John P. Snyder (1987), p. 157, eq. (21-2), (21-3), (21-4).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.stereo2latlon","page":"API reference","title":"FastIsostasy.stereo2latlon","text":"stereo2latlon(x, y, lat0, lon0)\n\nCompute the inverse stereographic projection (lat, lon) based on Cartesian coordinates (x,y) and for a given reference latitude lat0 and reference longitude lon0. Optionally one can provide x::AbstractMatrix and y::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be para elloprovided in degrees!\n\nConvert stereographic (x,y)-coordinates to latitude-longitude. Reference: John P. Snyder (1987), p. 159, eq. (20-14), (20-15), (20-18), (21-15).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_rigidity","page":"API reference","title":"FastIsostasy.get_rigidity","text":"get_rigidity(t::T, E::T, nu::T) where {T<:AbstractFloat}\n\nCompute rigidity D based on thickness t, Young modulus E and Poisson ration nu.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.loginterp_viscosity","page":"API reference","title":"FastIsostasy.loginterp_viscosity","text":"loginterp_viscosity(tvec, layer_viscosities, layers_thickness, pseudodiff)\n\nCompute a log-interpolator of the equivalent viscosity from provided viscosity fields layer_viscosities at time stamps tvec.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_elasticgreen","page":"API reference","title":"FastIsostasy.get_elasticgreen","text":"get_elasticgreen(Omega, quad_support, quad_coeffs)\n\nIntegrate load response over field by using 2D quadrature with specified support points and associated coefficients.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_quad_coeffs","page":"API reference","title":"FastIsostasy.get_quad_coeffs","text":"get_quad_coeffs(T, n)\n\nReturn support points and associated coefficients with specified Type for Gauss-Legendre quadrature.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.quadrature1D","page":"API reference","title":"FastIsostasy.quadrature1D","text":"quadrature1D(f, n, x1, x2)\n\nCompute 1D Gauss-Legendre quadrature of f between x1 and x2 based on n support points.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.kernelpromote","page":"API reference","title":"FastIsostasy.kernelpromote","text":"kernelpromote(X, arraykernel)\n\nPromote X to the kernel (Array or CuArray) specified by arraykernel.\n\n\n\n\n\n","category":"function"},{"location":"index.html#Introduction","page":"Introduction","title":"Introduction","text":"","category":"section"},{"location":"index.html#Getting-started","page":"Introduction","title":"Getting started","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FastIsostasy.jl is work under devlopment and is not a registered julia package so far. To install it, please run:","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"using Pkg\nPkg.add(\"https://github.com/JanJereczek/FastIsostasy.jl\")","category":"page"},{"location":"index.html#FastIsostasy.jl-–-For-whom?","page":"Introduction","title":"FastIsostasy.jl – For whom?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"This package is mainly addressed to ice sheet modellers looking for a regional model of glacial isostatic adjustment (GIA) that (1) captures the 3D structure of solid-Earth parameters, (2) computes an approximation of the sea-level equation, (3) runs kiloyear simulations on high resolution within minutes (without the need of HPC hardware) and (4) comes with ready-to-use calibration tools. For GIA \"purists\", this package is likely to miss interesting processes but we belive that the ridiculous run-time of FastIsostasy.jl can help them to perform some fast prototypting of a problem they might then transfer to a more comprehensive model.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"tip: Star us on GitHub!\nIf you have found this library useful, please consider starring it on GitHub. This gives us a lower bound of the satisfied user count.","category":"page"},{"location":"index.html#How-to-read-the-docs?","page":"Introduction","title":"How to read the docs?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"If you already know about GIA, skip to Overview of GIA for ice-sheet simulation. If you are already familiar with the complexity range of GIA models, skip to Why FastIsostasy?. If you want to have a more thorough but still very accessbile introduction to GIA, we highly recommend reading Whitehouse et al. 2018. If you want to get started right away, feel free to directly go to the Examples. If you face any problem using the code or want to know more about the functionalities of the package, visit the API reference. If you face a problem you cannot solve, please open a GitHub issue with a minimal and reproduceable example.","category":"page"},{"location":"index.html#What-is-glacial-isostatic-adjustment?","page":"Introduction","title":"What is glacial isostatic adjustment?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"The evolution of cryosphere components leads to changes in the ice and liquid water column and therefore in the vertical load applied upon the solid Earth. Glacial isostatic adjustment (GIA) denotes the mechanical response of the solid Earth, which is characterized by its vertical and horizontal displacement. GIA models usually encompass related processes, such as the resulting changes in sea-surface height and sea level.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"The magnitude and time scale of GIA depends on the applied load and on solid-Earth parameters, here assumed to be the density, the viscosity and the lithospheric thickness. These parameters display a radial and sometimes also a lateral variability, further jointly denoted by parameter \"heterogeneity\". For further details, please refer to Wiens et al. 2021 and Ivins et al. 2023.","category":"page"},{"location":"index.html#Why-should-we-care?","page":"Introduction","title":"Why should we care?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"GIA is known to present many feedbacks on ice-sheet evolution. Their net effect is negative, meaning that GIA inhibits ice-sheet growth and retreat. In other words, it tends to stabilize a given state and is therefore particularly important in the context of paleo-climate and climate change.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"The speed and magnitude of anthropogenic warming is a potential threat to the Greenland and the West-Antarctic ice sheets. They both represent an ice volume that could lead to multi-meter sea-level rise. The effect of GIA in this context appears to be particularly relevant - not only from a theoretical but also from a practical perspective, as a large portion of human livelihoods are concentrated along coasts.","category":"page"},{"location":"index.html#Motivation","page":"Introduction","title":"Motivation","text":"","category":"section"},{"location":"index.html#Overview-of-GIA-for-ice-sheet-simulation","page":"Introduction","title":"Overview of GIA for ice-sheet simulation","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"GIA models present a wide range of complexity, which can only be briefly mentioned here. On the lower end, models such as the Elastic-Lithopshere/Viscous-Asthenopshere are (1) cheap to run and (2) easy to implement, which has made them popular within the ice-sheet modelling community. They present some acceptable limitations such as (3) regionally approximating a global problem and (4) lacking the radially layered structure of the solid Earth. However, some limitations have shown to be too important to be overlooked – mainly the fact that (5) the heterogeneity of the lithospheric thickness and upper-mantle viscosity cannot be represented.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"On the higher end of the complexity spectrum, we find the 3D GIA models which address all the limitations of low-complexity models but are (1) expensive to run, (2) more tedious to couple to an ice-sheet model and (3) generally lack a well-documented and open-source code base. Due to these drawbacks, they do not represent a standard tool within the ice-sheet modelling community. Nonetheless, they are becoming increasingly used, as for instance in Gomez et al. 2018 and Van Calcar et al. 2023.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"We here willingly omit to speak about 1D GIA models, as they lack the representation of heterogeneous solid-Earth parameters.","category":"page"},{"location":"index.html#Where-is-FastIsosatsy.jl-on-the-complexity-range?","page":"Introduction","title":"Where is FastIsosatsy.jl on the complexity range?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"Although they are increasingly being coupled to ice-sheet models, we believe that the expense of 3D GIA models can be avoided while still addressing the aforementioned limitations of simplistic models. Models specifically designed for ice-sheet modelling, such as Bueler et al. 2007 and Coulon et al. 2021, have shown first improvements in closing the gap between simplistic and expensive models. FastIsostasy continues this work by generalizing both of these contributions into one, while benchmarking results against 1D and 3D GIA models.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FastIsostasy heavily relies on the Fast-Fourier Transform (FFT), as (1) its central PDE is solved by applying a Fourier collocation scheme and (2) important diagnostic fields are computed by matrix convolutions which can famously be accelerated by the use of FFT. FFT therefore inspired the name \"FastIsostasy\", along with a GitHub repository that eased the first steps of this package. The use of a performant language such as julia, as well as supporting performance-relevant computations on GPU allows FastIsostasy to live up to the expectations of low computation time.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"We believe that FastIsostasy drastically reduces the burdens associated with using a 3D GIA model while offering all the complexity needed for ice-sheet modelling. As targeted and efficient climate-change mitigation relies on a good representation of important mechanisms in numerical models, we believe that this can be a significant contribution for future research.","category":"page"},{"location":"index.html#Technical-details","page":"Introduction","title":"Technical details","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"In case you wonder, FastIsostasy.jl:","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"Takes all parameters in SI units. This might be made more flexible in future by the use of Unitful.jl.\nRelies on a regular, square grid as those typically used for finite-difference schemes.\nHas a hybrid approach to solving its underlying PDE: while some terms are evaluated by finite differences, the usual expense of such method is avoided by applying a Fourier collocation scheme.\nFor now only supports square domains with the number of points being a power of 2. This can accelerate computations of FFTs but will be made more flexible in future work.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FastIsostasy.jl largely relies on following packages:","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FFTW.jl\nInterpolations.jl\nCUDA.jl\nDSP.jl\nKalmanEnsembleProcesses.jl","category":"page"},{"location":"index.html#References","page":"Introduction","title":"References","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"Whitehouse et al. 2018\nWiens et al. 2021\nIvins et al. 2023.\nGomez et al. 2018\nVan Calcar et al. 2023\nBueler et al. 2007\nCoulon et al. 2021","category":"page"}] +[{"location":"examples.html#Examples","page":"Examples","title":"Examples","text":"","category":"section"},{"location":"examples.html#Multi-layer-Earth","page":"Examples","title":"Multi-layer Earth","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"FastIsostasy relies on a (polar) stereographic projection that allows to treat the radially-layered, onion-like structure of the solid Earth as a superposition of horizontal layers. Furthermore, FastIsostasy reduces this 3D problem into a 2D problem by collapsing the depth dimension by computing an effective viscosity field. The user is required to provide the 3D information, which will then be used under the hood to compute the effective viscosity. This tutorial shows such an example.","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"We want to render a situation similar to the one depicted below:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"(Image: Schematic representation of the three-layer set-up.)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Initializing a MultilayerEarth with parameters corresponding to this situation automatically computes the conversion from a 3D to a 2D problem. This can be simply executed by running:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = 3000e3 # (m) half-width of the domain\nn = 8 # implies an NxN grid with N = 2^n = 256.\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()\n\nlv = [1e19, 1e21] # (Pa*s)\nlb = [88e3, 400e3] # (m)\np = MultilayerEarth(Omega, c, layer_viscosities = lv, layer_boundaries = lb)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"The next section shows how to use the now obtained p::MultilayerEarth for actual GIA computation.","category":"page"},{"location":"examples.html#Simple-load-and-geometry","page":"Examples","title":"Simple load and geometry","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"We now apply a constant load, here a cylinder of ice with radius R = 1000 km and thickness H = 1 km, over the domain introduced in Multi-layer Earth. To obtain the bedrock displacement over time and store it at time steps specified by a vector t_out, we can use the convenience function fastisostasy:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"R = 1000e3 # ice disc radius (m)\nH = 1000.0 # ice disc thickness (m)\nHice = uniform_ice_cylinder(Omega, R, H)\nt_out = years2seconds.([0.0, 100.0, 500.0, 1500.0, 5000.0, 10_000.0, 50_000.0])\n\nresults = fastisostasy(t_out, Omega, c, p, Hice)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Yes, that was it! You can now easily access the elastic and viscous displacement by calling results.elastic or results.viscous. For the present case, the latter can be compared to an analytic solution that is known for this particular case. Let's look at the accuracy of our numerical scheme over time by running following plotting commands:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"M = Omega.N ÷ 2\nx = diag(Omega.X)[1:M]\ny = diag(Omega.Y)[1:M]\nr = sqrt.( x .^ 2 + y .^ 2 )\n\nfig = Figure()\nax3 = Axis(fig[1, 1])\ncolors = [:gray80, :gray65, :gray50, :gray35, :gray20, :gray5]\n\nfor i in eachindex(t_out)\n t = t_out[i]\n analytic_solution_r(r) = analytic_solution(r, t, c, p, H, R, analytic_support)\n u_analytic = analytic_solution_r.( r )\n u_numeric = diag(results.viscous[i])\n lines!(ax3, x, u_numeric[1:M], color = colors[i], linewidth = 5)\n lines!(ax3, x, u_analytic, color = colors[i], linewidth = 5, linestyle = :dash)\nend\nfig","category":"page"},{"location":"examples.html#Time-changing-load","page":"Examples","title":"Time-changing load","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"That looks pretty good! One might however object that the convenience function fastisostasy ends up being not so convenient as soon as the ice load changes over time. This case can however be easily handled by providing snapshots of the ice thickness and their assoicated time. By passing this to fastisostasy, an interpolator is created and called within the time integration. Let's create a tool example where the thickness of the ice cylinder asymptotically grows from 0 to 1 km, this can be implemented by:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"normalized_asymptote(t) = 1 - exp(-t)\nt_Hice_asymptotic = 0:10.0:t_out[end]\nHice_asymptotic = [normalized_asymptote(t) .* Hice for t in t_Hice_asymptote]\nresults = fastisostasy(t_out, Omega, c, p, t_Hice_asymptotic, Hice_asymptotic)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"This concept will also apply to the upper-mantle viscosity in future versions, as it can change over time.","category":"page"},{"location":"examples.html#GPU-support","page":"Examples","title":"GPU support","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"For about n 6, the previous example can be computed even faster by using GPU parallelism. It could not represent less work from the user's perspective, as it boils down to calling the ComputationDomain with an extra keyword argument:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Omega = ComputationDomain(W, n, use_cuda=true)","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"That's it, nothing more! One could suggest you lay back but your computation might be completed too soon for that.","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"info: Only CUDA supported!\nFor now only Nvidia GPUs are supported (sorry Mac users, destiny is taking its revenge upon you for being so cool) and there is no plan of extending this compatibility at this point.","category":"page"},{"location":"examples.html#Simple-load-and-geometry-DIY","page":"Examples","title":"Simple load and geometry - DIY","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"Nonetheless, as any high-level convenience function, fastisostasy has limitations. An ice-sheet modeller typically wants to embed FastIsostasy within a time-stepping loop. This can be easily done by getting familiar with some intermediate-level functions:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = 3000e3 # (m) half-width of the domain\nn = 8 # implies an NxN grid with N = 2^n = 256.\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()\np = MultilayerEarth(Omega, c)\n\nR = 1000e3 # ice disc radius (m)\nH = 1000.0 # ice disc thickness (m)\nHice = uniform_ice_cylinder(Omega, R, H)\nt_out = years2seconds.([0.0, 100.0, 500.0, 1500.0, 5000.0, 10_000.0, 50_000.0])\n\nresults = fastisostasy(t_out, Omega, c, p, Hice)","category":"page"},{"location":"examples.html#GIA-following-Antarctic-deglaciation","page":"Examples","title":"GIA following Antarctic deglaciation","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"We now want to provide a tough example that presents:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"a heterogeneous lithosphere thickness\na heterogeneous upper-mantle viscosity\nvarious viscous channels\na more elaborate load that evolves over time\nchanges in the sea-level","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"For this we run a deglaciation of Antarctica, based on the ice thickness estimated in GLAC1D.","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = 3000e3 # (m) half-width of the domain\nn = 8 # implies an NxN grid with N = 2^n = 256.\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()","category":"page"},{"location":"examples.html#Inversion-of-solid-Earth-parameters","page":"Examples","title":"Inversion of solid-Earth parameters","text":"","category":"section"},{"location":"examples.html","page":"Examples","title":"Examples","text":"FastIsostasy.jl relies on simplification of the full problem and might therefore need a calibration step to match the output of a 3D GIA model. By means of an unscented Kalman inversion, one can e.g. infer the appropriate effective upper-mantle viscosity based on the response of a 3D GIA model to a given load. Whereas this is know to be a tedious step, FastIsostasy is developped to ease the procedure by providing a convenience struct Paraminversion that can be run by:","category":"page"},{"location":"examples.html","page":"Examples","title":"Examples","text":"W = T(3000e3) # half-length of the square domain (m)\nOmega = ComputationDomain(W, n)\nc = PhysicalConstants()\n\nlb = [88e3, 180e3, 280e3, 400e3]\nlv = get_wiens_layervisc(Omega)\np = MultilayerEarth(\n Omega,\n c,\n layer_boundaries = lb,\n layer_viscosities = lv,\n)\nground_truth = copy(p.effective_viscosity)\n\nR = T(2000e3) # ice disc radius (m)\nH = T(1000) # ice disc thickness (m)\nHcylinder = uniform_ice_cylinder(Omega, R, H)\nt_out = years2seconds.(0.0:1_000.0:2_000.0)\n\nt1 = time()\nresults = fastisostasy(t_out, Omega, c, p, Hcylinder, ODEsolver=BS3(), interactive_geostate=false)\nt_fastiso = time() - t1\nprintln(\"Took $t_fastiso seconds!\")\nprintln(\"-------------------------------------\")\n\ntinv = t_out[2:end]\nHice = [Hcylinder for t in tinv]\nY = results.viscous[2:end]\nparaminv = ParamInversion(Omega, c, p, tinv, Y, Hice)\npriors, ukiobj = perform(paraminv)\nlogeta, Gx, e_mean, e_sort = extract_inversion(priors, ukiobj, paraminv)","category":"page"},{"location":"collected_garbage.html","page":"-","title":"-","text":"","category":"page"},{"location":"collected_garbage.html","page":"-","title":"-","text":"","category":"page"},{"location":"collected_garbage.html#Arrays","page":"-","title":"Arrays","text":"","category":"section"},{"location":"collected_garbage.html","page":"-","title":"-","text":"Array dimensions correspond to the spatial dimension of the variable they describe. If they evolve over time, they are stored as vector of arrays. For instance, the vertical displacement of the bedrock is a 2D variable that evolves over time. Therefore, it is stored in a Vector{Matrix}.","category":"page"},{"location":"APIref.html#API-reference","page":"API reference","title":"API reference","text":"","category":"section"},{"location":"APIref.html#Basic-structs","page":"API reference","title":"Basic structs","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"ComputationDomain\nPhysicalConstants\nMultilayerEarth\nRefSealevelState\nSealevelState\nPrecomputedFastiso\nSuperStruct\nFastisoResults","category":"page"},{"location":"APIref.html#FastIsostasy.ComputationDomain","page":"API reference","title":"FastIsostasy.ComputationDomain","text":"ComputationDomain\n\nReturn a struct containing all information related to geometry of the domain and potentially used parallelism.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.PhysicalConstants","page":"API reference","title":"FastIsostasy.PhysicalConstants","text":"PhysicalConstants\n\nReturn a struct containing important physical constants.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.MultilayerEarth","page":"API reference","title":"FastIsostasy.MultilayerEarth","text":"MultilayerEarth\n\nReturn a struct containing all information related to the radially layered structure of the solid Earth and its parameters.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.GeoState","page":"API reference","title":"FastIsostasy.GeoState","text":"GeoState\n\nReturn a mutable struct containing the geostate which will be updated over the simulation.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.PrecomputedFastiso","page":"API reference","title":"FastIsostasy.PrecomputedFastiso","text":"PrecomputedFastiso\n\nReturn a struct containing all the pre-computed terms needed for the forward integration of the model.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#FastIsostasy.SuperStruct","page":"API reference","title":"FastIsostasy.SuperStruct","text":"SuperStruct\n\nReturn a struct containing all the other structs needed for the forward integration of the model.\n\n\n\n\n\n","category":"type"},{"location":"APIref.html#Mechanics","page":"API reference","title":"Mechanics","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"ice_load\nfastisostasy\ninit_superstruct\nforward_isostasy\ninit_results\nforwardstep_isostasy!\ndudt_isostasy!\nsimple_euler!\napply_bc\ncompute_elastic_response","category":"page"},{"location":"APIref.html#FastIsostasy.ice_load","page":"API reference","title":"FastIsostasy.ice_load","text":"ice_load(c::PhysicalConstants{T}, H::T)\n\nCompute ice load based on ice thickness.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.fastisostasy","page":"API reference","title":"FastIsostasy.fastisostasy","text":"fastisostasy()\n\nMain function. List of all available solvers here.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.init_superstruct","page":"API reference","title":"FastIsostasy.init_superstruct","text":"init_superstruct()\n\nInit a SuperStruct containing all necessary fields for forward integration in time.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.dudt_isostasy!","page":"API reference","title":"FastIsostasy.dudt_isostasy!","text":"dudt_isostasy!()\n\nUpdate the displacement rate dudt.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#Sea-level","page":"API reference","title":"Sea-level","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"update_geostate!\nupdate_geoid!\nget_loadchange\nget_geoidgreen\nupdate_loadcolumns!\nupdate_sealevel!\nupdate_slc!\nupdate_V_af!\nupdate_slc_af!\nupdate_V_pov!\nupdate_slc_pov!\nupdate_V_den!\nupdate_slc_den!","category":"page"},{"location":"APIref.html#FastIsostasy.update_geostate!","page":"API reference","title":"FastIsostasy.update_geostate!","text":"update_geostate!(sstruct::SuperStruct, u::Matrix, H_ice::Matrix)\n\nUpdate the ::GeoState computing the current geoid perturbation, the sea-level changes and the load columns for the next time step of the isostasy integration.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_geoid!","page":"API reference","title":"FastIsostasy.update_geoid!","text":"update_geoid!(sstruct::SuperStruct)\n\nUpdate the geoid of a ::GeoState by convoluting the Green's function with the load change.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_loadchange","page":"API reference","title":"FastIsostasy.get_loadchange","text":"get_loadchange(sstruct::SuperStruct)\n\nCompute the load change compared to the reference configuration after Coulon et al. 2021. Difference to Coulon: HERE WE SCALE WITH CORRECTION FACTOR FOR PROJECTION!\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_geoidgreen","page":"API reference","title":"FastIsostasy.get_geoidgreen","text":"get_geoidgreen(sstruct::SuperStruct)\n\nReturn the Green's function used to compute the changes in geoid.\n\nReference\n\nCoulon et al. 2021.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_loadcolumns!","page":"API reference","title":"FastIsostasy.update_loadcolumns!","text":"update_loadcolumns!(sstruct::SuperStruct, u::XMatrix, H_ice::XMatrix)\n\nUpdate the load columns of a ::GeoState.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_sealevel!","page":"API reference","title":"FastIsostasy.update_sealevel!","text":"update_sealevel!(sstruct::SuperStruct)\n\nUpdate the sea-level ::GeoState by adding the various contributions.\n\nReference\n\nCoulon et al. (2021), Figure 1.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc!","page":"API reference","title":"FastIsostasy.update_slc!","text":"update_slc!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution of melting above floatation, density correction and potential ocean volume.\n\nReference\n\nGoelzer et al. (2020), eq. (12).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_V_af!","page":"API reference","title":"FastIsostasy.update_V_af!","text":"update_V_af!(sstruct::SuperStruct)\n\nUpdate the ice volume above floatation.\n\nReference\n\nGoelzer et al. (2020), eq. (13). Note: we do not use eq. (1) as it is only a special case of eq. (13) that does not allow a correct representation of external sea-level forcings.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc_af!","page":"API reference","title":"FastIsostasy.update_slc_af!","text":"update_slc_af!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution of ice above floatation.\n\nReference\n\nGoelzer et al. (2020), eq. (2).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_V_pov!","page":"API reference","title":"FastIsostasy.update_V_pov!","text":"update_V_pov!(sstruct::SuperStruct)\n\nUpdate the potential ocean volume.\n\nReference\n\nGoelzer et al. (2020), eq. (14). Note: we do not use eq. (8) as it is only a special case of eq. (14) that does not allow a correct representation of external sea-level forcinsstruct.geostate.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc_pov!","page":"API reference","title":"FastIsostasy.update_slc_pov!","text":"update_slc_pov!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution associated with the potential ocean volume.\n\nReference\n\nGoelzer et al. (2020), eq. (9).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_V_den!","page":"API reference","title":"FastIsostasy.update_V_den!","text":"update_V_den!(sstruct::SuperStruct)\n\nUpdate the ocean volume associated with the density correction.\n\nReference\n\nGoelzer et al. (2020), eq. (10).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.update_slc_den!","page":"API reference","title":"FastIsostasy.update_slc_den!","text":"update_slc_den!(sstruct::SuperStruct)\n\nUpdate the sea-level contribution associated with the density correction.\n\nReference\n\nGoelzer et al. (2020), eq. (11).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#Parameter-inversion","page":"API reference","title":"Parameter inversion","text":"","category":"section"},{"location":"APIref.html#Utils","page":"API reference","title":"Utils","text":"","category":"section"},{"location":"APIref.html","page":"API reference","title":"API reference","text":"years2seconds\nseconds2years\nm_per_sec2mm_per_yr\nmatrify_vectorconstant\nmatrify_constant\nget_r\nmeshgrid\ndist2angulardist\nsphericaldistance\nscalefactor\nlatlon2stereo\nstereo2latlon\nget_rigidity\nget_effective_viscosity\nget_differential_fourier\nget_viscosity_ratio\nthree_layer_scaling\nloginterp_viscosity\nhyperbolic_channel_coeffs\nget_greenintegrand_coeffs\nbuild_greenintegrand\nget_elasticgreen\nget_quad_coeffs\nquadrature1D\nquadrature2D\nget_normalized_lin_transform\nnormalized_lin_transform\nkernelpromote","category":"page"},{"location":"APIref.html#FastIsostasy.years2seconds","page":"API reference","title":"FastIsostasy.years2seconds","text":"years2seconds(t::Real)\n\nConvert input time t from years to seconds.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.seconds2years","page":"API reference","title":"FastIsostasy.seconds2years","text":"seconds2years(t::Real)\n\nConvert input time t from seconds to years.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.m_per_sec2mm_per_yr","page":"API reference","title":"FastIsostasy.m_per_sec2mm_per_yr","text":"m_per_sec2mm_per_yr(dudt::Real)\n\nConvert displacement rate dudt from $ m \\, s^{-1} $ to $ mm \\, \\mathrm{yr}^{-1} $.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.matrify","page":"API reference","title":"FastIsostasy.matrify","text":"matrify(x, N)\n\nGenerate a vector of constant matrices from a vector of constants.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_r","page":"API reference","title":"FastIsostasy.get_r","text":"get_r(x::T, y::T) where {T<:Real}\n\nGet euclidean distance of point (x, y) to origin.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.meshgrid","page":"API reference","title":"FastIsostasy.meshgrid","text":"meshgrid(x, y)\n\nReturn a 2D meshgrid spanned by x, y.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.dist2angulardist","page":"API reference","title":"FastIsostasy.dist2angulardist","text":"dist2angulardist(r::Real)\n\nConvert Euclidean to angular distance along great circle.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.latlon2stereo","page":"API reference","title":"FastIsostasy.latlon2stereo","text":"latlon2stereo(lat, lon, lat0, lon0)\n\nCompute stereographic projection (x,y) for a given latitude lat longitude lon, reference latitude lat0 and reference longitude lon0. Optionally one can provide lat::AbstractMatrix and lon::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be provided in degrees! Reference: John P. Snyder (1987), p. 157, eq. (21-2), (21-3), (21-4).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.stereo2latlon","page":"API reference","title":"FastIsostasy.stereo2latlon","text":"stereo2latlon(x, y, lat0, lon0)\n\nCompute the inverse stereographic projection (lat, lon) based on Cartesian coordinates (x,y) and for a given reference latitude lat0 and reference longitude lon0. Optionally one can provide x::AbstractMatrix and y::AbstractMatrix if the projection is to be computed for the whole domain. Note: angles must be para elloprovided in degrees!\n\nConvert stereographic (x,y)-coordinates to latitude-longitude. Reference: John P. Snyder (1987), p. 159, eq. (20-14), (20-15), (20-18), (21-15).\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_rigidity","page":"API reference","title":"FastIsostasy.get_rigidity","text":"get_rigidity(t::T, E::T, nu::T) where {T<:AbstractFloat}\n\nCompute rigidity D based on thickness t, Young modulus E and Poisson ration nu.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.loginterp_viscosity","page":"API reference","title":"FastIsostasy.loginterp_viscosity","text":"loginterp_viscosity(tvec, layer_viscosities, layers_thickness, pseudodiff)\n\nCompute a log-interpolator of the equivalent viscosity from provided viscosity fields layer_viscosities at time stamps tvec.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_elasticgreen","page":"API reference","title":"FastIsostasy.get_elasticgreen","text":"get_elasticgreen(Omega, quad_support, quad_coeffs)\n\nIntegrate load response over field by using 2D quadrature with specified support points and associated coefficients.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.get_quad_coeffs","page":"API reference","title":"FastIsostasy.get_quad_coeffs","text":"get_quad_coeffs(T, n)\n\nReturn support points and associated coefficients with specified Type for Gauss-Legendre quadrature.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.quadrature1D","page":"API reference","title":"FastIsostasy.quadrature1D","text":"quadrature1D(f, n, x1, x2)\n\nCompute 1D Gauss-Legendre quadrature of f between x1 and x2 based on n support points.\n\n\n\n\n\n","category":"function"},{"location":"APIref.html#FastIsostasy.kernelpromote","page":"API reference","title":"FastIsostasy.kernelpromote","text":"kernelpromote(X, arraykernel)\n\nPromote X to the kernel (Array or CuArray) specified by arraykernel.\n\n\n\n\n\n","category":"function"},{"location":"index.html#Introduction","page":"Introduction","title":"Introduction","text":"","category":"section"},{"location":"index.html#Getting-started","page":"Introduction","title":"Getting started","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FastIsostasy.jl is work under devlopment and is not a registered julia package so far. To install it, please run:","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"using Pkg\nPkg.add(\"https://github.com/JanJereczek/FastIsostasy.jl\")","category":"page"},{"location":"index.html#FastIsostasy.jl-–-For-whom?","page":"Introduction","title":"FastIsostasy.jl – For whom?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"This package is mainly addressed to ice sheet modellers looking for a regional model of glacial isostatic adjustment (GIA) that (1) captures the 3D structure of solid-Earth parameters, (2) computes an approximation of the sea-level equation, (3) runs kiloyear simulations on high resolution within minutes (without the need of HPC hardware) and (4) comes with ready-to-use calibration tools. For GIA \"purists\", this package is likely to miss interesting processes but we belive that the ridiculous run-time of FastIsostasy.jl can help them to perform some fast prototypting of a problem they might then transfer to a more comprehensive model.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"tip: Star us on GitHub!\nIf you have found this library useful, please consider starring it on GitHub. This gives us a lower bound of the satisfied user count.","category":"page"},{"location":"index.html#How-to-read-the-docs?","page":"Introduction","title":"How to read the docs?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"If you already know about GIA, skip to Overview of GIA for ice-sheet simulation. If you are already familiar with the complexity range of GIA models, skip to Why FastIsostasy?. If you want to have a more thorough but still very accessbile introduction to GIA, we highly recommend reading Whitehouse et al. 2018. If you want to get started right away, feel free to directly go to the Examples. If you face any problem using the code or want to know more about the functionalities of the package, visit the API reference. If you face a problem you cannot solve, please open a GitHub issue with a minimal and reproduceable example.","category":"page"},{"location":"index.html#What-is-glacial-isostatic-adjustment?","page":"Introduction","title":"What is glacial isostatic adjustment?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"The evolution of cryosphere components leads to changes in the ice and liquid water column and therefore in the vertical load applied upon the solid Earth. Glacial isostatic adjustment (GIA) denotes the mechanical response of the solid Earth, which is characterized by its vertical and horizontal displacement. GIA models usually encompass related processes, such as the resulting changes in sea-surface height and sea level.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"The magnitude and time scale of GIA depends on the applied load and on solid-Earth parameters, here assumed to be the density, the viscosity and the lithospheric thickness. These parameters display a radial and sometimes also a lateral variability, further jointly denoted by parameter \"heterogeneity\". For further details, please refer to Wiens et al. 2021 and Ivins et al. 2023.","category":"page"},{"location":"index.html#Why-should-we-care?","page":"Introduction","title":"Why should we care?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"GIA is known to present many feedbacks on ice-sheet evolution. Their net effect is negative, meaning that GIA inhibits ice-sheet growth and retreat. In other words, it tends to stabilize a given state and is therefore particularly important in the context of paleo-climate and climate change.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"The speed and magnitude of anthropogenic warming is a potential threat to the Greenland and the West-Antarctic ice sheets. They both represent an ice volume that could lead to multi-meter sea-level rise. The effect of GIA in this context appears to be particularly relevant - not only from a theoretical but also from a practical perspective, as a large portion of human livelihoods are concentrated along coasts.","category":"page"},{"location":"index.html#Motivation","page":"Introduction","title":"Motivation","text":"","category":"section"},{"location":"index.html#Overview-of-GIA-for-ice-sheet-simulation","page":"Introduction","title":"Overview of GIA for ice-sheet simulation","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"GIA models present a wide range of complexity, which can only be briefly mentioned here. On the lower end, models such as the Elastic-Lithopshere/Viscous-Asthenopshere are (1) cheap to run and (2) easy to implement, which has made them popular within the ice-sheet modelling community. They present some acceptable limitations such as (3) regionally approximating a global problem and (4) lacking the radially layered structure of the solid Earth. However, some limitations have shown to be too important to be overlooked – mainly the fact that (5) the heterogeneity of the lithospheric thickness and upper-mantle viscosity cannot be represented.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"On the higher end of the complexity spectrum, we find the 3D GIA models which address all the limitations of low-complexity models but are (1) expensive to run, (2) more tedious to couple to an ice-sheet model and (3) generally lack a well-documented and open-source code base. Due to these drawbacks, they do not represent a standard tool within the ice-sheet modelling community. Nonetheless, they are becoming increasingly used, as for instance in Gomez et al. 2018 and Van Calcar et al. 2023.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"We here willingly omit to speak about 1D GIA models, as they lack the representation of heterogeneous solid-Earth parameters.","category":"page"},{"location":"index.html#Where-is-FastIsosatsy.jl-on-the-complexity-range?","page":"Introduction","title":"Where is FastIsosatsy.jl on the complexity range?","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"Although they are increasingly being coupled to ice-sheet models, we believe that the expense of 3D GIA models can be avoided while still addressing the aforementioned limitations of simplistic models. Models specifically designed for ice-sheet modelling, such as Bueler et al. 2007 and Coulon et al. 2021, have shown first improvements in closing the gap between simplistic and expensive models. FastIsostasy continues this work by generalizing both of these contributions into one, while benchmarking results against 1D and 3D GIA models.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FastIsostasy heavily relies on the Fast-Fourier Transform (FFT), as (1) its central PDE is solved by applying a Fourier collocation scheme and (2) important diagnostic fields are computed by matrix convolutions which can famously be accelerated by the use of FFT. FFT therefore inspired the name \"FastIsostasy\", along with a GitHub repository that eased the first steps of this package. The use of a performant language such as julia, as well as supporting performance-relevant computations on GPU allows FastIsostasy to live up to the expectations of low computation time.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"We believe that FastIsostasy drastically reduces the burdens associated with using a 3D GIA model while offering all the complexity needed for ice-sheet modelling. As targeted and efficient climate-change mitigation relies on a good representation of important mechanisms in numerical models, we believe that this can be a significant contribution for future research.","category":"page"},{"location":"index.html#Technical-details","page":"Introduction","title":"Technical details","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"In case you wonder, FastIsostasy.jl:","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"Takes all parameters in SI units. This might be made more flexible in future by the use of Unitful.jl.\nRelies on a regular, square grid as those typically used for finite-difference schemes.\nHas a hybrid approach to solving its underlying PDE: while some terms are evaluated by finite differences, the usual expense of such method is avoided by applying a Fourier collocation scheme.\nFor now only supports square domains with the number of points being a power of 2. This can accelerate computations of FFTs but will be made more flexible in future work.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FastIsostasy.jl largely relies on following packages:","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"FFTW.jl\nInterpolations.jl\nCUDA.jl\nDSP.jl\nKalmanEnsembleProcesses.jl","category":"page"},{"location":"index.html#References","page":"Introduction","title":"References","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"Whitehouse et al. 2018\nWiens et al. 2021\nIvins et al. 2023.\nGomez et al. 2018\nVan Calcar et al. 2023\nBueler et al. 2007\nCoulon et al. 2021","category":"page"}] } diff --git a/docs/src/APIref.md b/docs/src/APIref.md index 1eda468..971561c 100644 --- a/docs/src/APIref.md +++ b/docs/src/APIref.md @@ -52,8 +52,7 @@ update_slc_den! years2seconds seconds2years m_per_sec2mm_per_yr -matrify_vectorconstant -matrify_constant +matrify get_r meshgrid dist2angulardist diff --git a/src/FastIsostasy.jl b/src/FastIsostasy.jl index 55efa22..b262e52 100644 --- a/src/FastIsostasy.jl +++ b/src/FastIsostasy.jl @@ -1,14 +1,14 @@ module FastIsostasy using LinearAlgebra +using Statistics: mean using Distributions: MvNormal using Interpolations: linear_interpolation, Flat -using FFTW -using StatsBase -using FastGaussQuadrature +using FFTW: fft, plan_fft, plan_ifft +using AbstractFFTs: Plan, ScaledPlan +using FastGaussQuadrature: gausslegendre using DSP: conv -using CUDA -using CUDA: CuArray +using CUDA: CuArray, CUFFT.plan_fft, CUFFT.plan_ifft, allowscalar using OrdinaryDiffEq: ODEProblem, solve, OrdinaryDiffEqAlgorithm using EnsembleKalmanProcesses using EnsembleKalmanProcesses.Observations @@ -16,10 +16,10 @@ using EnsembleKalmanProcesses.ParameterDistributions using Reexport @reexport using Interpolations -@reexport using OrdinaryDiffEq: Euler, BS3, Tsit5, TanYam7, Vern9, VCABM, Rosenbrock23, QNDF, FBDF, ImplicitEuler +@reexport using OrdinaryDiffEq: Euler, Midpoint, Heun, Ralston, BS3, BS5, RK4, OwrenZen3, + OwrenZen4, OwrenZen5, Tsit5, DP5, RKO65, TanYam7, DP8, Feagin10, Feagin12, Feagin14, + TsitPap8, Vern6, Vern7, Vern8, Vern9, VCABM, Rosenbrock23, QNDF, FBDF, ImplicitEuler -# Euler, Midpoint, Heun, Ralston, RK4, OwrenZen3, OwrenZen4, OwrenZen5, DP5, RKO65, -# TanYam7, DP8, Feagin10, Feagin12, Feagin14, TsitPap8, BS5, Vern6, Vern7, Vern8, # KuttaPRK2p5(dt=), Trapezoid(autodiff = false), PDIRK44(autodiff = false) include("structs.jl") @@ -30,7 +30,6 @@ include("mechanics.jl") include("inversion.jl") # structs.jl -export XMatrix export ComputationDomain export PhysicalConstants export MultilayerEarth @@ -39,30 +38,24 @@ export GeoState export SuperStruct # utils.jl -export years2seconds -export seconds2years -export m_per_sec2mm_per_yr +export years2seconds, seconds2years, m_per_sec2mm_per_yr +export meshgrid, dist2angulardist, latlon2stereo, stereo2latlon +export kernelpromote, convert2Array, copystructs2cpu -export meshgrid -export dist2angulardist -export latlon2stereo -export stereo2latlon - -export kernelpromote -export convert2Array -export copystructs2cpu - -export matrify_vectorconstant +export matrify, matrify export loginterp_viscosity export get_rigidity export get_r export gauss_distr +export samesize_conv + # derivatives.jl export mixed_fdx export mixed_fdy export mixed_fdxx export mixed_fdyy +export get_differential_fourier # mechanics.jl export quadrature1D diff --git a/src/derivatives.jl b/src/derivatives.jl index 77cdfa2..ad3f18b 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -1,81 +1,81 @@ # FDM in x, 1st order -function central_fdx(M::XMatrix, h::T) where {T<:AbstractFloat} +function central_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n2 = size(M, 2) return (view(M, :, 3:n2) - view(M, :, 1:n2-2)) ./ (2*h) end -function forward_fdx(M::XMatrix, h::T) where {T<:AbstractFloat} +function forward_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return (view(M, :, 2) - view(M, :, 1)) ./ h end -function backward_fdx(M::XMatrix, h::T) where {T<:AbstractFloat} +function backward_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n2 = size(M, 2) return (view(M, :, n2) - view(M, :, n2-1)) ./ h end # FDM in y, 1st order -function mixed_fdx(M::XMatrix, h::T) where {T<:AbstractFloat} +function mixed_fdx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return cat( forward_fdx(M,h), central_fdx(M,h), backward_fdx(M,h), dims=2 ) end -function central_fdy(M::XMatrix, h::T) where {T<:AbstractFloat} +function central_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n1 = size(M, 1) return (view(M, 3:n1, :) - view(M, 1:n1-2, :)) ./ (2*h) end -function forward_fdy(M::XMatrix, h::T) where {T<:AbstractFloat} +function forward_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return (view(M, 2, :) - view(M, 1, :)) ./ h end -function backward_fdy(M::XMatrix, h::T) where {T<:AbstractFloat} +function backward_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n1 = size(M, 1) return (view(M, n1, :) - view(M, n1-1, :)) ./ h end -function mixed_fdy(M::XMatrix, h::T) where {T<:AbstractFloat} +function mixed_fdy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return cat( forward_fdy(M,h)', central_fdy(M,h), backward_fdy(M,h)', dims=1 ) end # FDM in x, 2nd order -function central_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat} +function central_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n2 = size(M, 2) return (view(M, :, 3:n2) - 2 .* view(M, :, 2:n2-1) + view(M, :, 1:n2-2)) ./ h^2 end -function forward_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat} +function forward_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return (view(M, :, 3) - 2 .* view(M, :, 2) + view(M, :, 1)) ./ h^2 end -function backward_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat} +function backward_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n2 = size(M, 2) return (view(M, :, n2) - 2 .* view(M, :, n2-1) + view(M, :, n2-2)) ./ h^2 end -function mixed_fdxx(M::XMatrix, h::T) where {T<:AbstractFloat} +function mixed_fdxx(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return cat( forward_fdxx(M,h), central_fdxx(M,h), backward_fdxx(M,h), dims=2 ) end # FDM in y, 2nd order -function central_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat} +function central_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n1 = size(M, 1) return (view(M, 3:n1, :) - 2 .* view(M, 2:n1-1, :) + view(M, 1:n1-2, :)) ./ h^2 end -function forward_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat} +function forward_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return (view(M, 3, :) - 2 .* view(M, 2, :) + view(M, 1, :)) ./ h^2 end -function backward_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat} +function backward_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} n1 = size(M, 1) return (view(M, n1, :) - 2 .* view(M, n1-1, :) + view(M, n1-2, :)) ./ h^2 end -function mixed_fdyy(M::XMatrix, h::T) where {T<:AbstractFloat} +function mixed_fdyy(M::AbstractMatrix{T}, h::T) where {T<:AbstractFloat} return cat( forward_fdyy(M,h)', central_fdyy(M,h), backward_fdyy(M,h)', dims=1 ) end -function mixed_fdxy(M::XMatrix, hx::T, hy::T) where {T<:AbstractFloat} +function mixed_fdxy(M::AbstractMatrix{T}, hx::T, hy::T) where {T<:AbstractFloat} return mixed_fdy(mixed_fdx(M, hx), hy) end @@ -91,13 +91,11 @@ end Compute the matrices representing the differential operators in the fourier space. """ -function get_differential_fourier( - W::T, - N2::Int, -) where {T<:Real} - mu = T(π / W) - raw_coeffs = mu .* T.( vcat(0:N2, N2-1:-1:1) ) - x_coeffs, y_coeffs = raw_coeffs, raw_coeffs +function get_differential_fourier(Wx::T, Wy::T, Nx::Int, Ny::Int) where {T<:Real} + mu_x = π / Wx + mu_y = π / Wy + x_coeffs = mu_x .* fftint(Nx) + y_coeffs = mu_y .* fftint(Ny) X_coeffs, Y_coeffs = meshgrid(x_coeffs, y_coeffs) harmonic_coeffs = X_coeffs .^ 2 + Y_coeffs .^ 2 pseudodiff_coeffs = sqrt.(harmonic_coeffs) @@ -105,24 +103,11 @@ function get_differential_fourier( return pseudodiff_coeffs, harmonic_coeffs, biharmonic_coeffs end -function precomp_fourier_dxdy( - M::XMatrix, - L1::T, - L2::T, -) where {T<:AbstractFloat} - n1, n2 = size(M) - k1 = 2 * π / L1 .* vcat(0:n1/2-1, 0, -n1/2+1:-1) - k2 = 2 * π / L2 .* vcat(0:n2/2-1, 0, -n2/2+1:-1) - p1, p2 = plan_fft(k1), plan_fft(k2) - ip1, ip2 = plan_ifft(k1), plan_ifft(k2) - return k1, k2, p1, p2, ip1, ip2 -end - -function fourier_dnx( - M::XMatrix, - k1::Vector{T}, - p1::AbstractFFTs.Plan, - ip1::AbstractFFTs.ScaledPlan, -) where {T<:AbstractFloat} - return vcat( [real.(ip1 * ( ( im .* k1 ) .^ n .* (p1 * M[i, :]) ))' for i in axes(M,1)]... ) +function fftint(N::Int) + N2 = N ÷ 2 + if iseven(N) + return vcat(0:N2, N2-1:-1:1) + else + return vcat(0:N2, N2:-1:1) + end end \ No newline at end of file diff --git a/src/geostate.jl b/src/geostate.jl index 4115f95..46ffcdf 100644 --- a/src/geostate.jl +++ b/src/geostate.jl @@ -5,11 +5,8 @@ Update the geoid of a `::GeoState` by convoluting the Green's function with the load anom. """ function update_geoid!(sstruct::SuperStruct{<:AbstractFloat}) - sstruct.geostate.geoid .= view( - conv( sstruct.tools.geoidgreen, loadanom_green(sstruct) ), - sstruct.Omega.N2:2*sstruct.Omega.N-1-sstruct.Omega.N2, - sstruct.Omega.N2:2*sstruct.Omega.N-1-sstruct.Omega.N2, - ) + sstruct.geostate.geoid .= samesize_conv(sstruct.tools.geoidgreen, + loadanom_green(sstruct), sstruct.Omega) return nothing end @@ -109,12 +106,12 @@ end """ - update_loadcolumns!(sstruct::SuperStruct, u::XMatrix, H_ice::XMatrix) + update_loadcolumns!(sstruct::SuperStruct, u::AbstractMatrix{T}, H_ice::AbstractMatrix{T}) Update the load columns of a `::GeoState`. """ -function update_loadcolumns!(sstruct::SuperStruct{<:AbstractFloat}, - u::XMatrix, H_ice::XMatrix) +function update_loadcolumns!(sstruct::SuperStruct{T}, u::AbstractMatrix{T}, + H_ice::AbstractMatrix{T}) where {T<:AbstractFloat} sstruct.geostate.b .= sstruct.refgeostate.b .+ u sstruct.geostate.H_ice .= H_ice diff --git a/src/mechanics.jl b/src/mechanics.jl index 5c32b58..872b7d5 100644 --- a/src/mechanics.jl +++ b/src/mechanics.jl @@ -19,7 +19,7 @@ function fastisostasy( ODEsolver::Any = "ExplicitEuler", dt::T = T(years2seconds(1.0)), t_eta_snapshots::Vector{T} = [t_out[1], t_out[end]], - eta_snapshots::Vector{<:XMatrix} = [p.effective_viscosity, p.effective_viscosity], + eta_snapshots::Vector{<:AbstractMatrix{T}} = [p.effective_viscosity, p.effective_viscosity], u_viscous_0::Matrix{T} = copy(Omega.null), u_elastic_0::Matrix{T} = copy(Omega.null), interactive_geostate::Bool = false, @@ -62,7 +62,7 @@ Forward-integrate the isostatic adjustment. function forward_isostasy( dt::T, t_out::Vector{T}, - u::XMatrix, + u::AbstractMatrix{T}, sstruct::SuperStruct{T}, ODEsolver::Any, verbose::Bool, @@ -78,7 +78,12 @@ function forward_isostasy( end if isa(ODEsolver, OrdinaryDiffEqAlgorithm) prob = ODEProblem(forwardstep_isostasy!, u, (t0, t_out[k]), sstruct) - sol = solve(prob, ODEsolver, reltol=1e-3) #, dtmin = years2seconds(0.1)), , dt=dt + + if ODEsolver == Euler() + sol = solve(prob, ODEsolver, reltol=1e-3, dt = dt) + else + sol = solve(prob, ODEsolver, reltol=1e-3) #, dtmin = years2seconds(0.1)) + end u .= sol(t_out[k], Val{0}) dudt .= sol(t_out[k], Val{1}) @@ -93,6 +98,8 @@ function forward_isostasy( dudt_out[k] .= copy(kernelpromote(dudt, Array)) geoid_out[k] .= copy(kernelpromote(sstruct.geostate.geoid, Array)) sealevel_out[k] .= copy(kernelpromote(sstruct.geostate.sealevel, Array)) + u_el_out[k] .= samesize_conv(columnanom_load(sstruct), + sstruct.tools.elasticgreen, sstruct.Omega) end return u_out, dudt_out, u_el_out, geoid_out, sealevel_out end @@ -122,12 +129,12 @@ Forward integrate the isostatic adjustment over a single time step by updating t displacement rate `dudt`, and the sea-level state contained within `sstruct::SuperStruct`. """ function forwardstep_isostasy!( - dudt::XMatrix, - u::XMatrix, + dudt::AbstractMatrix{T}, + u::AbstractMatrix{T}, sstruct::SuperStruct{T}, t::T, ) where {T<:AbstractFloat} - corner_bc!(u, sstruct.Omega.N) + corner_bc!(u, sstruct.Omega.Nx, sstruct.Omega.Ny) # Only update the geoid and sea level if geostate is interactive. # As integration requires smaller time steps than diagnostics, @@ -140,7 +147,20 @@ function forwardstep_isostasy!( # println("Updated GeoState at t=$(seconds2years(t))") end + coupled_viscoelastic = false + + if coupled_viscoelastic + sstruct.geostate.dtloadanom .= columnanom_load(sstruct) + end + update_loadcolumns!(sstruct, u, sstruct.Hice(t)) + + if coupled_viscoelastic + sstruct.geostate.dtloadanom -= columnanom_load(sstruct) + u += samesize_conv(sstruct.geostate.dtloadanom, sstruct.tools.elasticgreen, + sstruct.Omega) + end + dudt_isostasy!(dudt, u, sstruct, t) return nothing end @@ -152,8 +172,8 @@ end Update the displacement rate `dudt`. """ function dudt_isostasy!( - dudt::XMatrix, - u::XMatrix, + dudt::AbstractMatrix{T}, + u::AbstractMatrix{T}, sstruct::SuperStruct{T}, t::T, ) where {T<:AbstractFloat} @@ -178,7 +198,7 @@ function dudt_isostasy!( sstruct.Omega.pseudodiff)) ./ (2 .* sstruct.p.effective_viscosity) # dudt[:, :] .= real.(tools.pifft * ((tools.pfft * (rhs ./ (2 .* p.effective_viscosity)) ) ./ Omega.pseudodiff)) - corner_bc!(dudt, sstruct.Omega.N) + corner_bc!(dudt, sstruct.Omega.Nx, sstruct.Omega.Ny) return nothing end @@ -190,8 +210,8 @@ Update the state `u` by performing an explicit Euler integration of its derivati over a time step `dt`. """ function explicit_euler!( - u::XMatrix, - dudt::XMatrix, + u::AbstractMatrix{T}, + dudt::AbstractMatrix{T}, dt::T, ) where {T<:AbstractFloat} u .+= dudt .* dt @@ -209,14 +229,14 @@ end Apply boundary condition on Fourier collocation solution. Assume that mean deformation at corners of domain is 0. """ -function corner_bc(u::XMatrix, N) +function corner_bc(u::AbstractMatrix{<:AbstractFloat}, Nx::Int, Ny::Int) u_bc = copy(u) - return corner_bc!(u_bc, N) + return corner_bc!(u_bc, Nx, Ny) end -function corner_bc!(u::XMatrix, N) - CUDA.allowscalar() do - u .-= (view(u, 1, 1) + view(u, 1, N) + view(u, N, 1) + view(u, N, N)) / 4 +function corner_bc!(u::AbstractMatrix{<:AbstractFloat}, Nx::Int, Ny::Int) + allowscalar() do + u .-= ( view(u,1,1) + view(u,1,Nx) + view(u,Ny,1) + view(u,Ny,Nx) ) / 4 end end @@ -228,14 +248,15 @@ Apply boundary condition on Fourier collocation solution. Assume that mean deformation at borders of domain is 0. Same as Bueler et al. (2007). """ -function border_bc(u::XMatrix, N) +function border_bc(u::AbstractMatrix{<:AbstractFloat}, Nx::Int, Ny::Int) u_bc = copy(u) - return border_bc!(u_bc, N) + return border_bc!(u_bc, Nx, Ny) end -function border_bc!(u::XMatrix, N) - CUDA.allowscalar() do - u .-= sum(view(u, 1, :) + view(u, :, N) + view(u, N, :) + view(u, :, 1)) / (4*N) +function border_bc!(u::AbstractMatrix{<:AbstractFloat}, Nx::Int, Ny::Int) + allowscalar() do + u .-= sum( view(u,1,:) + view(u,:,Nx) + view(u,Ny,:) + view(u,:,1) ) / + (2*Nx + 2*Ny) end end @@ -251,9 +272,9 @@ by convoluting the `load` with the Green's function stored in the pre-computed ` (elements obtained from Farell 1972). """ function compute_elastic_response( - Omega::ComputationDomain, - tools::PrecomputedFastiso, - load::XMatrix, -) - return conv(load, tools.elasticgreen)[Omega.N2:end-Omega.N2, Omega.N2:end-Omega.N2] + Omega::ComputationDomain{T}, + tools::PrecomputedFastiso{T}, + load::AbstractMatrix{T}, +) where {T<:AbstractFloat} + return samesize_conv(load, tools.elasticgreen, Omega) end \ No newline at end of file diff --git a/src/structs.jl b/src/structs.jl index c9a2100..cbea469 100644 --- a/src/structs.jl +++ b/src/structs.jl @@ -1,4 +1,9 @@ -XMatrix = Union{Matrix{T}, CuArray{T, 2}} where {T<:Real} +KernelMatrix = Union{Matrix{T}, CuArray{T, 2}} where {T<:Real} +# FFTPlan = Union{cFFTWPlan{ComplexF64, -1, false, 2, UnitRange{Int64}}, +# CUFFT.cCuFFTPlan{ComplexF64, -1, false, 2}} +# IFFTPlan = Union{ +# AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, false, 2, UnitRange{Int64}}, Float64}, +# AbstractFFTs.ScaledPlan{ComplexF64, CUDA.CUFFT.cCuFFTPlan{ComplexF64, 1, false, 2}, Float64}} ######################################################### # Computation domain @@ -16,31 +21,41 @@ Omega = ComputationDomain(W, n) struct ComputationDomain{T<:AbstractFloat} Wx::T # Domain length in x (m) Wy::T # Domain length in y (m) - N::Int # Average number of grid points along one dimension - N2::Int # N/2 + Nx::Int # Number of grid points in x-dimension + Ny::Int # Number of grid points in y-dimension + Mx::Int # Nx/2 + My::Int # Ny/2 dx::T # Spatial discretization in x dy::T # Spatial discretization in y x::Vector{T} y::Vector{T} - X::XMatrix - Y::XMatrix - R::XMatrix - Theta::XMatrix - Lat::XMatrix - Lon::XMatrix - K::XMatrix + X::AbstractMatrix{T} + Y::AbstractMatrix{T} + R::AbstractMatrix{T} + Theta::AbstractMatrix{T} + Lat::AbstractMatrix{T} + Lon::AbstractMatrix{T} + K::AbstractMatrix{T} projection_correction::Bool - null::XMatrix # a zero matrix of size Nx x Ny - pseudodiff::XMatrix # pseudodiff operator - harmonic::XMatrix # harmonic operator - biharmonic::XMatrix # biharmonic operator + null::AbstractMatrix{T} # a zero matrix of size Nx x Ny + pseudodiff::AbstractMatrix{T} # pseudodiff operator + harmonic::AbstractMatrix{T} # harmonic operator + biharmonic::AbstractMatrix{T} # biharmonic operator use_cuda::Bool - arraykernel # Array or CuArray depending on chosen hardware + arraykernel # Array or CuArray depending on chosen hardware +end + +function ComputationDomain(W::T, n::Int; kwargs...,) where {T<:AbstractFloat} + Wx, Wy = W, W + Nx, Ny = 2^n, 2^n + return ComputationDomain(Wx, Wy, Nx, Ny; kwargs...) end function ComputationDomain( - W::T, - n::Int; + Wx::T, + Wy::T, + Nx::Int, + Ny::Int; use_cuda::Bool = false, lat0::T = T(-71.0), lon0::T = T(0.0), @@ -48,11 +63,9 @@ function ComputationDomain( ) where {T<:AbstractFloat} # Geometry - Wx, Wy = W, W - N = 2^n - N2 = Int(floor(N/2)) - dx = T(2*Wx) / N - dy = T(2*Wy) / N + Mx, My = Nx ÷ 2, Ny ÷ 2 + dx = 2*Wx / Nx + dy = 2*Wy / Ny x = collect(-Wx+dx:dx:Wx) y = collect(-Wy+dy:dy:Wy) X, Y = meshgrid(x, y) @@ -63,23 +76,19 @@ function ComputationDomain( arraykernel = use_cuda ? CuArray : Array K = kernelpromote(scalefactor(deg2rad.(Lat), deg2rad.(Lon), deg2rad(lat0), deg2rad(lon0)), arraykernel) - null = fill(T(0.0), N, N) + null = matrify(T(0), Nx, Ny) # Differential operators in Fourier space - pseudodiff, harmonic, biharmonic = get_differential_fourier(W, N2) + pseudodiff, harmonic, biharmonic = get_differential_fourier(Wx, Wy, Nx, Ny) pseudodiff[1, 1] = mean([pseudodiff[1,2], pseudodiff[2,1]]) pseudodiff, harmonic, biharmonic = kernelpromote( [pseudodiff, harmonic, biharmonic], arraykernel) # Avoid division by zero. Tolerance ϵ of the order of the neighboring terms. # Tests show that it does not lead to errors wrt analytical or benchmark solutions. - return ComputationDomain( - Wx, Wy, N, N2, - dx, dy, x, y, - X, Y, R, Theta, - Lat, Lon, K, projection_correction, null, - pseudodiff, harmonic, biharmonic, - use_cuda, arraykernel, + return ComputationDomain(Wx, Wy, Nx, Ny, Mx, My, dx, dy, x, y, X, Y, + R, Theta, Lat, Lon, K, projection_correction, null, + pseudodiff, harmonic, biharmonic, use_cuda, arraykernel, ) end @@ -120,9 +129,9 @@ its parameters. mutable struct MultilayerEarth{T<:AbstractFloat} mean_gravity::T mean_density::T - effective_viscosity::XMatrix - litho_thickness::XMatrix - litho_rigidity::XMatrix + effective_viscosity::AbstractMatrix{T} + litho_thickness::AbstractMatrix{T} + litho_rigidity::AbstractMatrix{T} litho_poissonratio::T layers_density::Vector{T} layer_viscosities::Array{T, 3} @@ -158,15 +167,15 @@ function MultilayerEarth( T<:AbstractFloat, A<:Union{Vector{T}, Array{T, 3}}, B<:Union{Vector{T}, Array{T, 3}}, - C<:Union{T, XMatrix}, - D<:Union{T, XMatrix}, + C<:Union{T, AbstractMatrix{T}}, + D<:Union{T, AbstractMatrix{T}}, } if layer_boundaries isa Vector{<:Real} - layer_boundaries = matrify_vectorconstant(layer_boundaries, Omega.N) + layer_boundaries = matrify(layer_boundaries, Omega.Nx, Omega.Ny) end if layer_viscosities isa Vector{<:Real} - layer_viscosities = matrify_vectorconstant(layer_viscosities, Omega.N) + layer_viscosities = matrify(layer_viscosities, Omega.Nx, Omega.Ny) end litho_thickness = layer_boundaries[:, :, 1] @@ -184,7 +193,7 @@ function MultilayerEarth( layers_thickness, ) - mean_density = fill(mean(layers_density), Omega.N, Omega.N) + mean_density = matrify(mean(layers_density), Omega.Nx, Omega.Ny) litho_rigidity, effective_viscosity, mean_density = kernelpromote( [litho_rigidity, effective_viscosity, mean_density], Omega.arraykernel) @@ -209,11 +218,11 @@ end Return a struct containing the reference geostate. We define the geostate to be all quantities related to sea-level. """ struct RefGeoState{T<:AbstractFloat} - H_ice::XMatrix # reference height of ice column - H_water::XMatrix # reference height of water column - b::XMatrix # reference bedrock position - z0::XMatrix # reference height to allow external sea-level forcing - sealevel::XMatrix # reference sealevel field + H_ice::AbstractMatrix{T} # reference height of ice column + H_water::AbstractMatrix{T} # reference height of water column + b::AbstractMatrix{T} # reference bedrock position + z0::AbstractMatrix{T} # reference height to allow external sea-level forcing + sealevel::AbstractMatrix{T} # reference sealevel field sle_af::T # reference sl-equivalent of ice volume above floatation V_pov::T # reference potential ocean volume V_den::T # reference potential ocean volume associated with V_den @@ -226,11 +235,11 @@ end Return a mutable struct containing the geostate which will be updated over the simulation. """ mutable struct GeoState{T<:AbstractFloat} - H_ice::XMatrix # current height of ice column - H_water::XMatrix # current height of water column - b::XMatrix # vertical bedrock position - geoid::XMatrix # current geoid displacement - sealevel::XMatrix # current sealevel field + H_ice::AbstractMatrix{T} # current height of ice column + H_water::AbstractMatrix{T} # current height of water column + b::AbstractMatrix{T} # vertical bedrock position + geoid::AbstractMatrix{T} # current geoid displacement + sealevel::AbstractMatrix{T} # current sealevel field V_af::T # ice volume above floatation sle_af::T # sl-equivalent of ice volume above floatation slc_af::T # sl-contribution of Vice above floatation @@ -241,6 +250,7 @@ mutable struct GeoState{T<:AbstractFloat} slc::T # total sealevel contribution countupdates::Int # count the updates of the geostate dt::T # update step + dtloadanom::AbstractMatrix{T} # load anomaly wrt previous time step end """ @@ -248,32 +258,32 @@ end PrecomputedFastiso(Omega::ComputationDomain, c::PhysicalConstants, p::MultilayerEarth) Return a `struct` containing pre-computed tools to perform forward-stepping of the model, namely: - - elasticgreen::XMatrix - - fourier_elasticgreen::XMatrix{Complex{T}} + - elasticgreen::AbstractMatrix{T} + - fourier_elasticgreen::AbstractMatrix{T}{Complex{T}} - pfft::AbstractFFTs.Plan - pifft::AbstractFFTs.ScaledPlan - - Dx::XMatrix - - Dy::XMatrix - - Dxx::XMatrix - - Dyy::XMatrix - - Dxy::XMatrix + - Dx::AbstractMatrix{T} + - Dy::AbstractMatrix{T} + - Dxx::AbstractMatrix{T} + - Dyy::AbstractMatrix{T} + - Dxy::AbstractMatrix{T} - negligible_gradD::Bool - rhog::T - - geoidgreen::XMatrix + - geoidgreen::AbstractMatrix{T} """ struct PrecomputedFastiso{T<:AbstractFloat} - elasticgreen::XMatrix - fourier_elasticgreen::XMatrix{Complex{T}} - pfft::AbstractFFTs.Plan - pifft::AbstractFFTs.ScaledPlan - Dx::XMatrix - Dy::XMatrix - Dxx::XMatrix - Dyy::XMatrix - Dxy::XMatrix + elasticgreen::AbstractMatrix{T} + fourier_elasticgreen::AbstractMatrix{Complex{T}} + pfft::Plan + pifft::ScaledPlan + Dx::AbstractMatrix{T} + Dy::AbstractMatrix{T} + Dxx::AbstractMatrix{T} + Dyy::AbstractMatrix{T} + Dxy::AbstractMatrix{T} negligible_gradD::Bool rhog::T - geoidgreen::XMatrix + geoidgreen::AbstractMatrix{T} end @@ -298,7 +308,7 @@ function PrecomputedFastiso( Dyy = mixed_fdyy(D, Omega.dy) Dxy = mixed_fdy( mixed_fdx(D, Omega.dx), Omega.dy ) - omega_zeros = fill(T(0.0), Omega.N, Omega.N) + omega_zeros = matrify(T(0), Omega.Nx, Omega.Ny) zero_tol = 1e-2 negligible_gradD = isapprox(Dx, omega_zeros, atol = zero_tol) & isapprox(Dy, omega_zeros, atol = zero_tol) & @@ -364,7 +374,7 @@ function SuperStruct( t_Hice_snapshots::Vector{T}, Hice_snapshots::Vector{Matrix{T}}, t_eta_snapshots::Vector{T}, - eta_snapshots::Vector{<:XMatrix}, + eta_snapshots::Vector{<:AbstractMatrix{T}}, interactive_geostate::Bool; geoid_0::Matrix{T} = copy(Omega.null), sealevel_0::Matrix{T} = copy(Omega.null), @@ -406,6 +416,7 @@ function SuperStruct( T(0.0), T(0.0), # V_den terms T(0.0), # total sl-contribution & conservation term 0, years2seconds(10.0), # countupdates, update step + copy(Omega.null), # dtloadanom ) return SuperStruct(Omega, c, p, tools, Hice, Hice_cpu, eta, eta_cpu, refgeostate, geostate, interactive_geostate) diff --git a/src/utils.jl b/src/utils.jl index 5451221..abe6c4d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -40,18 +40,44 @@ end """ - matrify_vectorconstant(x, N) + matrify(x, Nx, Ny) Generate a vector of constant matrices from a vector of constants. """ -function matrify_vectorconstant(x::Vector{T}, N::Int) where {T<:AbstractFloat} - X = zeros(T, N, N, length(x)) +function matrify(x::Vector{<:Real}, N::Int) + return matrify(x, N, N) +end + +function matrify(x::Vector{T}, Nx::Int, Ny::Int) where {T<:Real} + X = zeros(T, Ny, Nx, length(x)) @inbounds for i in eachindex(x) - X[:, :, i] = fill(x[i], N, N) + X[:, :, i] = matrify(x[i], Nx, Ny) end return X end +matrify(x::Real, Nx::Int, Ny::Int) = fill(x, Ny, Nx) + + +function samesize_conv(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, + Omega::ComputationDomain) where {T<:AbstractFloat} + if iseven(Omega.Ny) + i1 = Omega.My + else + i1 = Omega.My+1 + end + i2 = 2*Omega.Ny-1-Omega.My + + if iseven(Omega.Nx) + j1 = Omega.Mx + else + j1 = Omega.Mx+1 + end + j2 = 2*Omega.Nx-1-Omega.Mx + + return view( conv(X, Y), i1:i2, j1:j2 ) +end + ##################################################### # Domain and projection utils ##################################################### @@ -101,11 +127,11 @@ function scalefactor(lat::T, lon::T, lat0::T, lon0::T; k0::T = T(1)) where {T<:R return 2*k0 / (1 + sin(lat0)*sin(lat) + cos(lat0)*cos(lat)*cos(lon-lon0)) end -function scalefactor(lat::XMatrix, lon::XMatrix, lat0::T, lon0::T; +function scalefactor(lat::AbstractMatrix{T}, lon::AbstractMatrix{T}, lat0::T, lon0::T; kwargs... ) where {T<:Real} K = similar(lat) @inbounds for idx in CartesianIndices(lat) - K[idx] = scalefactor(lat[idx], lon[idx], lat0, lon0, kwargs...) + K[idx] = scalefactor(lat[idx], lon[idx], lat0, lon0; kwargs...) end return K end @@ -124,22 +150,22 @@ Reference: John P. Snyder (1987), p. 157, eq. (21-2), (21-3), (21-4). function latlon2stereo(lat::T, lon::T, lat0::T, lon0::T; R::T = T(6.371e6), kwargs...) where {T<:Real} lat, lon, lat0, lon0 = deg2rad.([lat, lon, lat0, lon0]) - k = scalefactor(lat, lon, lat0, lon0, kwargs...) + k = scalefactor(lat, lon, lat0, lon0; kwargs...) x = R * k * cos(lat) * sin(lon - lon0) y = R * k * (cos(lat0) * sin(lat) - sin(lat0) * cos(lat) * cos(lon-lon0)) return k, x, y end function latlon2stereo( - lat::XMatrix, - lon::XMatrix, + lat::AbstractMatrix{T}, + lon::AbstractMatrix{T}, lat0::T, lon0::T; kwargs..., ) where {T<:Real} K, X, Y = similar(lat), similar(lat), similar(lat) @inbounds for idx in CartesianIndices(lat) - K[idx], X[idx], Y[idx] = latlon2stereo(lat[idx], lon[idx], lat0, lon0, kwargs...) + K[idx], X[idx], Y[idx] = latlon2stereo(lat[idx], lon[idx], lat0, lon0; kwargs...) end return K, X, Y end @@ -167,16 +193,11 @@ function stereo2latlon(x::T, y::T, lat0::T, lon0::T; return rad2deg(lat), rad2deg(lon) end -function stereo2latlon( - x::XMatrix, - y::XMatrix, - lat0::T, - lon0::T; - kwargs..., -) where {T<:Real} +function stereo2latlon(x::AbstractMatrix{T}, y::AbstractMatrix{T}, lat0::T, lon0::T; + kwargs...) where {T<:Real} Lat, Lon = copy(x), copy(x) @inbounds for idx in CartesianIndices(x) - Lat[idx], Lon[idx] = stereo2latlon(x[idx], y[idx], lat0, lon0) + Lat[idx], Lon[idx] = stereo2latlon(x[idx], y[idx], lat0, lon0; kwargs...) end return Lat, Lon end @@ -188,13 +209,13 @@ end """ - gauss_distr(X::XMatrix, Y::XMatrix, mu::Vector{<:Real}, sigma::Matrix{<:Real}) + gauss_distr(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, mu::Vector{<:Real}, sigma::Matrix{<:Real}) Compute `Z = f(X,Y)` with `f` a Gaussian function parametrized by mean `mu` and covariance `sigma`. """ -function gauss_distr(X::XMatrix, Y::XMatrix, - mu::Vector{<:Real}, sigma::Matrix{<:Real}) +function gauss_distr(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, + mu::Vector{<:Real}, sigma::Matrix{<:Real}) where {T<:AbstractFloat} k = length(mu) G = similar(X) invsigma = inv(sigma) @@ -219,7 +240,7 @@ end """ get_effective_viscosity( - layer_viscosities::Vector{XMatrix}, + layer_viscosities::Vector{AbstractMatrix{T}}, layers_thickness::Vector{T}, Omega::ComputationDomain{T}, ) where {T<:AbstractFloat} @@ -231,7 +252,7 @@ function get_effective_viscosity( Omega::ComputationDomain{T}, layer_viscosities::Array{T, 3}, layers_thickness::Array{T, 3}, - # pseudodiff::XMatrix, + # pseudodiff::AbstractMatrix{T}, ) where {T<:AbstractFloat} # Recursion has to start with half space = n-th layer: @@ -296,10 +317,10 @@ function loginterp_viscosity( tvec::AbstractVector{T}, layer_viscosities::Array{T, 4}, layers_thickness::Array{T, 3}, - pseudodiff::XMatrix, + pseudodiff::AbstractMatrix{T}, ) where {T<:AbstractFloat} n1, n2, n3, nt = size(layer_viscosities) - log_eqviscosity = [fill(T(0.0), n1, n2) for k in 1:nt] + log_eqviscosity = matrify(zeros(nt), n1, n2) [log_eqviscosity[k] .= log10.(get_effective_viscosity( layer_viscosities[:, :, :, k], @@ -412,21 +433,20 @@ function get_elasticgreen( quad_coeffs::Vector{T}, ) where {T<:AbstractFloat} - h = Omega.dx - N = Omega.N + dx, dy = Omega.dx, Omega.dy elasticgreen = similar(Omega.X) - @inbounds for i = 1:N, j = 1:N - p = i - Omega.N2 - 1 - q = j - Omega.N2 - 1 - elasticgreen[i, j] = quadrature2D( + @inbounds for i = 1:Omega.Nx, j = 1:Omega.Ny + p = i - Omega.Mx - 1 + q = j - Omega.My - 1 + elasticgreen[j, i] = quadrature2D( greenintegrand_function, quad_support, quad_coeffs, - p*h, - p*h+h, - q*h, - q*h+h, + p*dx, + (p+1)*dx, + q*dy, + (q+1)*dy, ) end return elasticgreen @@ -560,8 +580,7 @@ function copystructs2cpu( p::MultilayerEarth{T}, ) where {T<:AbstractFloat} - n = Int( round( log2(Omega.N) ) ) - Omega_cpu = ComputationDomain(Omega.Wx, n, use_cuda = false) + Omega_cpu = ComputationDomain(Omega.Wx, Omega.Wy, Omega.Nx, Omega.Ny, use_cuda = false) p_cpu = MultilayerEarth( Omega_cpu, diff --git a/test/external_load_maps.jl b/test/external_load_maps.jl index 82283a7..42d85a9 100644 --- a/test/external_load_maps.jl +++ b/test/external_load_maps.jl @@ -16,7 +16,7 @@ function interpolated_glac1d_snapshots(Omega::ComputationDomain{T}) where {T<:Ab xl = range(xl[1], stop = xl[end], length = length(xl)) yl = range(yl[1], stop = yl[end], length = length(yl)) - H = fill(0.0, Omega.N, Omega.N, length(tvecl)) + H = matrify(zeros(length(tvecl)), Omega.Nx, Omega.Ny) for k in axes(H, 3) itp = extrapolate(scale(interpolate( Hl[:, :, k], diff --git a/test/helpers_compute.jl b/test/helpers_compute.jl index adf4163..a77e43f 100644 --- a/test/helpers_compute.jl +++ b/test/helpers_compute.jl @@ -1,10 +1,10 @@ using LinearAlgebra -function mask_disc(X::XMatrix, Y::XMatrix, R::T) where {T<:AbstractFloat} +function mask_disc(X::AbstractMatrix{T}, Y::AbstractMatrix{T}, R::T) where {T<:AbstractFloat} return mask_disc(sqrt.(X.^2 + Y.^2), R) end -function mask_disc(r::XMatrix, R::T) where {T<:AbstractFloat} +function mask_disc(r::AbstractMatrix{T}, R::T) where {T<:AbstractFloat} return T.(r .< R) end @@ -156,51 +156,19 @@ end # Generate binary parameter fields for test 3 ################################################ -function generate_binary_field( - Omega::ComputationDomain{T}, - x_lo::T, - x_hi::T, -) where {T<:AbstractFloat} - lo_half = fill(x_lo, Omega.N, Int(Omega.N/2)) - hi_half = fill(x_hi, Omega.N, Int(Omega.N/2)) - return cat(lo_half, hi_half, dims=2) -end - -function generate_sigmoid_field( - Omega::ComputationDomain{T}, - x_lo::T, - x_hi::T, -) where {T<:AbstractFloat} - - delta_x = x_hi - x_lo - scaled_sigmoid(x) = sigmoid(x, 1e-2) - return x_lo .+ delta_x .* scaled_sigmoid.(Omega.X) -end -sigmoid(x::Real, c::Real) = 1 / (1 + exp(-c*x)) - -function generate_window_field( - Omega::ComputationDomain{T}, - x_background::T, - x_lo::T, - x_hi::T, -) where {T<:AbstractFloat} - N = Omega.N - N2 = Int(N/2) - N4 = Int(N/4) - X = fill(x_background, N, N) - X[N4+1:end-N4, N4+1:N2] .= x_lo - X[N4+1:end-N4, N2+1:end-N4] .= x_hi - return X -end - function generate_gaussian_field( Omega::ComputationDomain{T}, z_background::T, xy_peak::Vector{T}, z_peak::T, - sigma::XMatrix, + sigma::AbstractMatrix{T}, ) where {T<:AbstractFloat} - N = Omega.N + if Omega.Nx == Omega.Ny + N = Omega.Nx + else + error("Automated generation of Gaussian parameter fields only supported for" * + "square domains.") + end G = gauss_distr( Omega.X, Omega.Y, xy_peak, sigma ) G = G ./ maximum(G) .* z_peak return fill(z_background, N, N) + G diff --git a/test/helpers_plot.jl b/test/helpers_plot.jl index 0181ee5..2efbc17 100644 --- a/test/helpers_plot.jl +++ b/test/helpers_plot.jl @@ -59,7 +59,7 @@ function plot_response( width = Relative(0.8), ) end - plotname = "plots/test1/2D/$(case)_N$(Omega.N)" + plotname = "plots/test1/2D/$(case)_Nx$(Omega.Nx)_Ny$(Omega.Ny)" save("$plotname.png", fig) save("$plotname.pdf", fig) return fig diff --git a/test/test1_compute.jl b/test/test1_compute.jl index db8921b..12719ba 100644 --- a/test/test1_compute.jl +++ b/test/test1_compute.jl @@ -16,7 +16,7 @@ function main( p = MultilayerEarth(Omega, c) kernel = use_cuda ? "gpu" : "cpu" - println("Computing on $kernel and $(Omega.N) x $(Omega.N) grid...") + println("Computing on $kernel and $(Omega.Ny) x $(Omega.Nx) grid...") R = T(1000e3) # ice disc radius (m) H = T(1000) # ice disc thickness (m) @@ -35,7 +35,13 @@ function main( end gs = active_gs ? "geostate" : "isostate" - filename = "$(solver)_N$(Omega.N)_$(kernel)_$(gs)" + if solver == BS3() + solvername = "BS3" + elseif solver == "ExplicitEuler" + solvername = "ExplicitEuler" + end + + filename = "$(solvername)_Nx$(Omega.Nx)_Ny$(Omega.Ny)_$(kernel)_$(gs)" jldsave( "data/test1/$filename.jld2", Omega = Omega, c = c, p = p, @@ -45,7 +51,6 @@ function main( ) end -# ["ExplicitEuler", BS3(), VCABM(), Rosenbrock23(autodiff=false)] for use_cuda in [false] # [false, true] for active_gs in [false] # [false, true] for n in 6:6 # 3:8 @@ -53,3 +58,23 @@ for use_cuda in [false] # [false, true] end end end + +#= +Slight speed up if using powers of 2: + +This file: +main(n, use_cuda = false, solver = BS3(), active_gs = false) +Took 0.6100420951843262 seconds! + +main(n, use_cuda = false, solver = "ExplicitEuler", active_gs = false) +Took 14.107969999313354 seconds! + +------------------------------------ + +test1_rectangle.jl: +main(63, 64, use_cuda = false, solver = BS3(), active_gs = false) +Took 0.6303250789642334 seconds! + +main(63, 64, use_cuda = false, solver = "ExplicitEuler", active_gs = false) +Took 14.486158847808838 seconds! +=# \ No newline at end of file diff --git a/test/test1_rectangle.jl b/test/test1_rectangle.jl new file mode 100644 index 0000000..826af65 --- /dev/null +++ b/test/test1_rectangle.jl @@ -0,0 +1,50 @@ +push!(LOAD_PATH, "../") +using FastIsostasy +using JLD2 +include("helpers_compute.jl") + +function main( + Nx::Int, + Ny::Int; + use_cuda::Bool = false, + solver::Any = "ExplicitEuler", + active_gs::Bool = true, +) + T = Float64 + W = T(3000e3) # half-length of the square domain (m) + Omega = ComputationDomain(W, W, Nx, Ny, use_cuda = use_cuda) + c = PhysicalConstants() + p = MultilayerEarth(Omega, c) + + kernel = use_cuda ? "gpu" : "cpu" + println("Computing on $kernel and $(Omega.Ny) x $(Omega.Nx) grid...") + + R = T(1000e3) # ice disc radius (m) + H = T(1000) # ice disc thickness (m) + Hcylinder = uniform_ice_cylinder(Omega, R, H) + t_out = years2seconds.([0.0, 100.0, 500.0, 1500.0, 5000.0, 10_000.0, 50_000.0]) + + t1 = time() + results = fastisostasy(t_out, Omega, c, p, Hcylinder, ODEsolver=solver, + interactive_geostate=active_gs) + t_fastiso = time() - t1 + println("Took $t_fastiso seconds!") + println("-------------------------------------") + + if use_cuda + Omega, p = copystructs2cpu(Omega, c, p) + end + + gs = active_gs ? "geostate" : "isostate" + filename = "$(solver)_Nx$(Omega.Nx)_Ny$(Omega.Ny)_$(kernel)_$(gs)" + jldsave( + "data/test1/$filename.jld2", + Omega = Omega, c = c, p = p, + results = results, + t_fastiso = t_fastiso, + R = R, H = H, + ) +end + +main(63, 64, use_cuda = false, solver = "ExplicitEuler", active_gs = false) + diff --git a/test/test2_compute.jl b/test/test2_compute.jl index 19ee6ab..e940f0e 100644 --- a/test/test2_compute.jl +++ b/test/test2_compute.jl @@ -28,7 +28,7 @@ function main( ) kernel = use_cuda ? "gpu" : "cpu" - println("Computing on $kernel and $(Omega.N) x $(Omega.N) grid...") + println("Computing on $kernel and $(Omega.Nx) x $(Omega.Ny) grid...") if occursin("disc", case) alpha = T(10) # max latitude (°) of uniform ice disc @@ -42,7 +42,7 @@ function main( end t_out = years2seconds.([0.0, 1.0, 1e3, 2e3, 5e3, 1e4, 1e5]) - sl0 = fill(-Inf, Omega.N, Omega.N) + sl0 = matrify(-Inf, Omega.Nx, Omega.Ny) t1 = time() results = fastisostasy(t_out, Omega, c, p, H_ice, sealevel_0=sl0, interactive_geostate=true) t_fastiso = time() - t1 @@ -53,7 +53,7 @@ function main( Omega, p = copystructs2cpu(Omega, c, p) end - filename = "$(case)_N$(Omega.N)_$(kernel)" + filename = "$(case)_Nx$(Omega.Nx)_Ny$(Omega.Ny)_$(kernel)" jldsave( "data/test2/$filename.jld2", Omega = Omega, c = c, p = p, diff --git a/test/test3_cases.jl b/test/test3_cases.jl index 50831ae..503e344 100644 --- a/test/test3_cases.jl +++ b/test/test3_cases.jl @@ -6,7 +6,7 @@ function choose_case(case::String, Omega::ComputationDomain, c::PhysicalConstant W = (Omega.Wx + Omega.Wy) / 2 sigma = diagm([(W/4)^2, (W/4)^2]) layer1_begin = generate_gaussian_field(Omega, 150e3, [0.0, 0.0], -100e3, sigma) - layer2_begin = fill(250e3, Omega.N, Omega.N) + layer2_begin = matrify(250e3, Omega.Nx, Omega.Ny) lb = cat(layer1_begin, layer2_begin, dims=3) p = MultilayerEarth(Omega, c, layer_boundaries = lb, layer_viscosities = [eta, eta], layers_density = [rho]) @@ -14,7 +14,7 @@ function choose_case(case::String, Omega::ComputationDomain, c::PhysicalConstant W = (Omega.Wx + Omega.Wy) / 2 sigma = diagm([(W/4)^2, (W/4)^2]) layer1_begin = generate_gaussian_field(Omega, 150e3, [0.0, 0.0], 100e3, sigma) - layer2_begin = fill(250e3, Omega.N, Omega.N) + layer2_begin = matrify(250e3, Omega.Nx, Omega.Ny) lb = cat(layer1_begin, layer2_begin, dims=3) p = MultilayerEarth(Omega, c, layer_boundaries = lb, layer_viscosities = [eta, eta], layers_density = [rho]) @@ -28,7 +28,7 @@ function choose_case(case::String, Omega::ComputationDomain, c::PhysicalConstant W = (Omega.Wx + Omega.Wy) / 2 sigma = diagm([(W/4)^2, (W/4)^2]) gauss_visc = 10.0 .^ generate_gaussian_field(Omega, 21.0, [0.0, 0.0], 1.0, sigma) - halfspace_viscosity = fill(1e21, Omega.N, Omega.N) + halfspace_viscosity = matrify(1e21, Omega.Nx, Omega.Ny) lv = cat(gauss_visc, halfspace_viscosity, dims=3) p = MultilayerEarth(Omega, c, layer_viscosities = lv) end diff --git a/test/test3_compute.jl b/test/test3_compute.jl index 3388566..0a00a76 100644 --- a/test/test3_compute.jl +++ b/test/test3_compute.jl @@ -18,7 +18,7 @@ function main( p = choose_case(case, Omega, c) kernel = use_cuda ? "gpu" : "cpu" - println("Computing on $kernel and $(Omega.N) x $(Omega.N) grid...") + println("Computing on $kernel and $(Omega.Nx) x $(Omega.Ny) grid...") R = T(1000e3) # ice disc radius (m) H = T(1000) # ice disc thickness (m) @@ -44,7 +44,7 @@ function main( Omega, p = copystructs2cpu(Omega, c, p) end - filename = "$(case)_$(kernel)_N$(Omega.N)_$densekey" + filename = "$(case)_$(kernel)_Nx$(Omega.Nx)_Ny$(Omega.Ny)_$densekey" jldsave( "data/test3/$filename.jld2", Omega = Omega, c = c, p = p, diff --git a/test/test4_compute.jl b/test/test4_compute.jl index d7e2d02..6992cf9 100644 --- a/test/test4_compute.jl +++ b/test/test4_compute.jl @@ -13,8 +13,8 @@ function main(n::Int, case::String; use_cuda::Bool = true, solver = "ExplicitEul Omega = ComputationDomain(W, n) # domain parameters c = PhysicalConstants() if occursin("homogeneous", case) - channel_viscosity = fill(1e20, Omega.N, Omega.N) - halfspace_viscosity = fill(1e21, Omega.N, Omega.N) + channel_viscosity = matrify(1e20, Omega.Nx, Omega.Ny) + halfspace_viscosity = matrify(1e21, Omega.Nx, Omega.Ny) lv = cat(channel_viscosity, halfspace_viscosity, dims=3) p = MultilayerEarth( Omega, @@ -24,7 +24,7 @@ function main(n::Int, case::String; use_cuda::Bool = true, solver = "ExplicitEul elseif occursin("meanviscosity", case) log10_eta_channel = interpolate_visc_wiens_on_grid(Omega.X, Omega.Y) channel_viscosity = 10 .^ (log10_eta_channel) - halfspace_viscosity = fill(1e21, Omega.N, Omega.N) + halfspace_viscosity = matrify(1e21, Omega.Nx, Omega.Ny) lv = cat(channel_viscosity, halfspace_viscosity, dims=3) p = MultilayerEarth( Omega, @@ -33,7 +33,7 @@ function main(n::Int, case::String; use_cuda::Bool = true, solver = "ExplicitEul ) elseif occursin("scaledviscosity", case) lb = [88e3, 180e3, 280e3, 400e3] - halfspace_logviscosity = fill(21.0, Omega.N, Omega.N) + halfspace_logviscosity = matrify(21.0, Omega.Nx, Omega.Ny) Eta, Eta_mean, z = load_wiens_2021(Omega.X, Omega.Y) eta_interpolators, eta_mean_interpolator = interpolate_viscosity_xy( @@ -53,7 +53,7 @@ function main(n::Int, case::String; use_cuda::Bool = true, solver = "ExplicitEul end kernel = use_cuda ? "gpu" : "cpu" - println("Computing on $kernel and $(Omega.N) x $(Omega.N) grid...") + println("Computing on $kernel and $(Omega.Nx) x $(Omega.Ny) grid...") R = T(2000e3) # ice disc radius (m) H = T(1000) # ice disc thickness (m) @@ -74,7 +74,7 @@ function main(n::Int, case::String; use_cuda::Bool = true, solver = "ExplicitEul end jldsave( - "data/test4/discload_nonlocal_$(case)_N$(Omega.N).jld2", + "data/test4/discload_nonlocal_$(case)_Nx$(Omega.Nx)_Ny$(Omega.Ny).jld2", Omega = Omega, c = c, p = p, results = results, t_fastiso = t_fastiso, diff --git a/test/test5_compute.jl b/test/test5_compute.jl index 6d03f62..3c28360 100644 --- a/test/test5_compute.jl +++ b/test/test5_compute.jl @@ -13,7 +13,7 @@ function main(n::Int, active_gs::Bool; use_cuda::Bool = false,solver = "Explicit c = PhysicalConstants() lb = [88e3, 180e3, 280e3, 400e3] - halfspace_logviscosity = fill(21.0, Omega.N, Omega.N) + halfspace_logviscosity = matrify(21.0, Omega.Nx, Omega.Ny) Eta, Eta_mean, z = load_wiens_2021(Omega.X, Omega.Y) eta_interpolators, eta_mean_interpolator = interpolate_viscosity_xy( Omega.X, Omega.Y, Eta, Eta_mean) @@ -30,7 +30,7 @@ function main(n::Int, active_gs::Bool; use_cuda::Bool = false,solver = "Explicit ) kernel = use_cuda ? "gpu" : "cpu" - println("Computing on $kernel and $(Omega.N) x $(Omega.N) grid...") + println("Computing on $kernel and $(Omega.Nx) x $(Omega.Ny) grid...") t_out, deltaH, H = interpolated_glac1d_snapshots(Omega) dH = [deltaH[:, :, k] for k in axes(deltaH, 3)] @@ -50,7 +50,7 @@ function main(n::Int, active_gs::Bool; use_cuda::Bool = false,solver = "Explicit case = active_gs ? "geostate" : "isostate" jldsave( - "data/test5/$(case)_N$(Omega.N).jld2", + "data/test5/$(case)_Nx$(Omega.Nx)_Ny$(Omega.Ny).jld2", Omega = Omega, c = c, p = p, results = results, t_fastiso = t_fastiso, @@ -60,5 +60,5 @@ end cases = [false, true] for active_gs in cases[1:1] - main(6, active_gs, use_cuda=true, solver="ExplicitEuler") + main(6, active_gs, use_cuda=false, solver="ExplicitEuler") end \ No newline at end of file diff --git a/test/test5_plot.jl b/test/test5_plot.jl index 2aed678..00fbb71 100644 --- a/test/test5_plot.jl +++ b/test/test5_plot.jl @@ -13,7 +13,8 @@ function main( ) N = 2^n - sol = load("data/test5/$(case)_N$(N).jld2") + filekey = "$(case)_Nx$(N)_Ny$(N)" + sol = load("data/test5/$filekey.jld2") results = sol["results"] t_out = results.t_out t_out_kyr = round.(seconds2years.(t_out) ./ 1e3, digits=1) @@ -177,7 +178,7 @@ function main( ) framerate = 24 - plotname = "plots/test5/loaduplift_$(case)_N$(N)" + plotname = "plots/test5/loaduplift_$filekey" record(fig, "$plotname.mp4", eachindex(u), framerate = framerate) do k i[] = k end @@ -188,7 +189,7 @@ end cases = ["isostate", "geostate"] for case in cases[1:1] - main(7, case) + main(6, case) end diff --git a/test/test6_compute.jl b/test/test6_compute.jl index 78e27b0..a3e3ada 100644 --- a/test/test6_compute.jl +++ b/test/test6_compute.jl @@ -4,7 +4,7 @@ include("helpers_compute.jl") include("external_viscosity_maps.jl") function get_wiens_layervisc(Omega) - halfspace_logviscosity = fill(21.0, Omega.N, Omega.N) + halfspace_logviscosity = matrify(21.0, Omega.Nx, Omega.Ny) Eta, Eta_mean, z = load_wiens_2021(Omega.X, Omega.Y) eta_interpolators, eta_mean_interpolator = interpolate_viscosity_xy( Omega.X, Omega.Y, Eta, Eta_mean) diff --git a/test/test_viscosityscaling.jl b/test/test_viscosityscaling.jl index e54c6af..a246448 100644 --- a/test/test_viscosityscaling.jl +++ b/test/test_viscosityscaling.jl @@ -15,8 +15,8 @@ W = T(3000e3) # half-length of the square domain (m) Omega = ComputationDomain(W, n) # domain parameters c = PhysicalConstants() if occursin("homogeneous", case) - channel_viscosity = fill(1e20, Omega.N, Omega.N) - halfspace_viscosity = fill(1e21, Omega.N, Omega.N) + channel_viscosity = matrify(1e20, Omega.Nx, Omega.Ny) + halfspace_viscosity = matrify(1e21, Omega.Nx, Omega.Ny) lv = cat(channel_viscosity, halfspace_viscosity, dims=3) p = MultilayerEarth( Omega, @@ -26,7 +26,7 @@ if occursin("homogeneous", case) elseif occursin("meanviscosity", case) log10_eta_channel = interpolate_visc_wiens_on_grid(Omega.X, Omega.Y) channel_viscosity = 10 .^ (log10_eta_channel) - halfspace_viscosity = fill(1e21, Omega.N, Omega.N) + halfspace_viscosity = matrify(1e21, Omega.Nx, Omega.Ny) lv = cat(channel_viscosity, halfspace_viscosity, dims=3) p = MultilayerEarth( Omega, @@ -35,7 +35,7 @@ elseif occursin("meanviscosity", case) ) elseif occursin("scaledviscosity", case) lb = [88e3, 180e3, 280e3, 400e3] - halfspace_logviscosity = fill(21.0, Omega.N, Omega.N) + halfspace_logviscosity = matrify(21.0, Omega.Nx, Omega.Ny) Eta, Eta_mean, z = load_wiens_2021(Omega.X, Omega.Y) eta_interpolators, eta_mean_interpolator = interpolate_viscosity_xy(