From fc56021c6686022063180031c0885a0c6f4df0bb Mon Sep 17 00:00:00 2001 From: Martin Hermann Date: Thu, 22 Jul 2021 13:59:24 +0200 Subject: [PATCH 1/6] fix small bug in adaptiveTOCollocationStiffnessMastrix --- src/TO.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/TO.jl b/src/TO.jl index 5bcc1f5e..30a511a9 100644 --- a/src/TO.jl +++ b/src/TO.jl @@ -170,6 +170,7 @@ function adaptiveTOCollocationStiffnessMatrix(ctx::GridContext{2}, flow_maps, ti new_ctx,new_bdata,new_density_bcdofvals = adaptiveTOFutureGrid(ctx, flow_map_t; on_torus=on_torus, on_cylinder=on_cylinder, + quadrature_order=quadrature_order, LL_future=LL_future, UR_future=UR_future, bdata=bdata, From 334552b8a35ce4aca442655c9db351a14a729e6d Mon Sep 17 00:00:00 2001 From: Martin Hermann Date: Thu, 22 Jul 2021 14:14:02 +0200 Subject: [PATCH 2/6] add basic functionality for dynamic isoperimetry --- src/dynamicIsoperimetry.jl | 241 +++++++++++++++++++++++++++++++++++++ src/exports.jl | 8 +- 2 files changed, 248 insertions(+), 1 deletion(-) create mode 100644 src/dynamicIsoperimetry.jl diff --git a/src/dynamicIsoperimetry.jl b/src/dynamicIsoperimetry.jl new file mode 100644 index 00000000..e79a4b6d --- /dev/null +++ b/src/dynamicIsoperimetry.jl @@ -0,0 +1,241 @@ +# functions are all for 2D only, as well as not really working with periodic boundaries + +function apply2curve(flowmap, curve::Curve2{T}) where {T} + moved_points = [flowmap(x) for x in zip(coordinates(curve)...)] + moved_points = [flowmap(x) for x in zip(coordinates(curve)...)] + return Curve2(moved_points) +end + +function get_length(curve, tensorfield=(x,y) -> 1) + xs, ys = coordinates(curve) + n = length(xs) + result = 0 + for j in 1:(n-1) + current_x, current_y = xs[j], ys[j] + next_x, next_y = xs[j+1], ys[j+1] + approxTangentVector = [next_x - current_x, next_y - current_y] + result += sqrt(approxTangentVector ⋅ (tensorfield(curx, cury) * approxTangentVector)) + end + return result +end + + +# If the first point is not the last point, we might get vol(Ω)-area instead +function get_euclidean_area(ctx,curve;tolerance=0.0) + + curve = close_curve(ctx,curve,tolerance=tolerance) + + xs, ys = coordinates(curve) + center_x, center_y = mean(xs), mean(ys) + n = length(xs) + result = 0.0 + for j in 1:(n-1) + current_x, current_y = xs[j], ys[j] + next_x, next_y = xs[j+1], ys[j+1] + + x1 = current_x - center_x + x2 = next_x - center_x + y1 = current_y - center_y + y2 = next_y - center_y + + area = 0.5*det([x2 y2; x1 y1]) + result +=area + end + return abs(result) +end + + +# from contour.jl we only get the guarantee that either the first and last points are equal, or the contour starts and +# ends at the boundary of Ω. in the last case, we will have to construct a closed curve e.g. in order to determine the +# enclosed area +function close_curve(ctx,curve::Curve2{T};tolerance=0.0,orientation=:clockwise) where {T} + xs, ys = coordinates(curve) + start_x, start_y, end_x, end_y = xs[1], ys[1], xs[end], ys[end] + + if (start_x == end_x && start_y == end_y) + return curve + end + + points = [SA[x,y] for (x,y) in zip(xs,ys)] + + LL, UR = ctx.spatialBounds + #Curve2 wants StaticArrays + LL, UR = SA[LL[1], LL[2]], SA[UR[1], UR[2]] + LR, UL = SA[UR[1], LL[2]], SA[LL[1], UR[2]] + + # We have a rectangle, and order its sides like + # 3 3 + # -- -- + # 2| |4 for clockwise, 4| |2 for anticlockwise orientation. + # -- -- + # 1 1 + # + # corners[i] then is the next corner in the respective orientation. + + if orientation==:clockwise + corners = [LL, UL, UR, LR] + elseif orientation ==:anticlockwise + corners = [LR, UR, UL, LL] + else + error("Unknown orientation '",orientation,"'!") + end + + if abs(end_y-corners[1][2])<=tolerance + side = 1 + elseif abs(end_x-corners[2][1])<=tolerance + side = 2 + elseif abs(end_y-corners[3][2])<=tolerance + side = 3 + elseif abs(end_x-corners[4][1])<=tolerance + side = 4 + else # point does not lie on the boundary + error("A curve has either to be closed or start and end at the boundary!") + end + + i = 0 + while(abs(start_x-end_x)>tolerance && abs(start_y-end_y)>tolerance) + push!(points,corners[side]) + end_x = corners[side][1] + end_y = corners[side][2] + side = (side%4)+1 + i += 1 + if i>4 # might happen if the first point is not on the border, prevent infinite loops + error("A curve has either to be closed or start and end at the boundary!") + end + end + push!(points,points[1]) # close curve + return Curve2(points) +end + + +function get_function_values(ctx, u; x_resolution=nothing, y_resolution=nothing, bdata=BoundaryData()) + + if isnothing(x_resolution) + x_resolution = ctx.numberOfPointsInEachDirection[1] + end + if isnothing(y_resolution) + y_resolution = ctx.numberOfPointsInEachDirection[2] + end + + xs = range(ctx.spatialBounds[1][1], stop=ctx.spatialBounds[2][1], length=x_resolution) + ys = range(ctx.spatialBounds[1][2], stop=ctx.spatialBounds[2][2], length=y_resolution) + + u_dofvals = undoBCS(ctx, u, bdata) + u_nodevals = u_dofvals[ctx.node_to_dof] + + fs = [evaluate_function_from_node_or_cellvals(ctx, u_nodevals, Vec{2}((x,y))) + for x in xs, y in ys] + + return xs,ys,fs +end + +""" + get_levelset(ctx, u, c; + x_resolution=nothing, y_resolution=nothing, bdata=BoundaryData()) + +Returns the levelset of u corresponding to the function value c +""" +function get_levelset(ctx, u, c; + x_resolution=nothing, y_resolution=nothing, bdata=BoundaryData()) + + xs, ys, fs = get_function_values(ctx,u,x_resolution=x_resolution,y_resolution=y_resolution,bdata=bdata) + + return Contour.contour(xs, ys, fs, c) +end + + + +""" + get_minimal_levelset(ctx, u, objective_function; + x_resolution=nothing, y_resolution=nothing, min=minimum(u), max=maximum(u), + n_candidates=100, bdata=BoundaryData()) + +Returns the levelset of u that achieves the minimal value of objective_function, considering levelsets +of values between min and max. +""" +function get_minimal_levelset(ctx, u, objective_function; + x_resolution=nothing, y_resolution=nothing, min=minimum(u), max=maximum(u), + n_candidates=100, bdata=BoundaryData()) + + xs, ys, fs = get_function_values(ctx,u,x_resolution=x_resolution,y_resolution=y_resolution,bdata=bdata) + + currentmin = Inf + result = nothing + + for cl in levels(contours(xs, ys, fs, range(min,stop=max,length=n_candidates))) + curves = lines(cl) + if length(curves) == 0 # this can happen e.g. for the levelset of max(u) + continue + end + # TODO + #if length(curves) != 1 + # @warn "Currently only connected levelsets are allowed! Levelset: ", level(cl) + #end + value = objective_function(curves[1]) + if value < currentmin + currentmin = value + result = cl + end + end + return result, currentmin +end + +""" + dynamic_cheeger_value(ctx, curve, flowmap; tolerance=0.0) + +Calculate the dynamic cheeger value of a curve, given a flowmap T(x). The gridcontext ctx is needed +to close curves that intersect with the boundary, the tolerance paramter determines how close to +the boundary a curve can end to still be considered to have a point on it. +""" +function dynamic_cheeger_value(ctx, curve, flowmap; tolerance=0.0) + LL, UR = ctx.spatialBounds + volume_Ω = (UR[1] - LL[1]) * (UR[2] - LL[2]) + curve_closed_cw = close_curve(ctx,curve,tolerance=tolerance,orientation=:clockwise) + image_curve_cw = apply2curve(flowmap,curve_closed_cw) + curve_closed_acw = close_curve(ctx,curve,tolerance=tolerance,orientation=:anticlockwise) + image_curve_acw = apply2curve(flowmap,curve_closed_acw) + combined_length = min(get_length(curve_closed_cw) + get_length(image_curve_cw), + get_length(curve_closed_acw) + get_length(image_curve_acw)) + volume_curve = get_euclidean_area(ctx,curve,tolerance=tolerance) + return 0.5 * combined_length / min(volume_curve, (volume_Ω - volume_curve)) +end + + +# simple finite differences for visualization only +# shouldn't be too hard to use the grid instead and do something like evaluate_function_from_node_or_cellvals +function gradient_from_values(xs,ys,fs) + dx = diff(fs,dims=1)./diff(xs) + dy = (diff(fs,dims=2)'./diff(ys))' + #simply extend by constant + dx = [dx ; dx[end,:]'] + dy = [dy dy[:,end]] + return dx, dy +end + +function get_gradient_field(ctx, u; + x_resolution=nothing, y_resolution=nothing, bdata=BoundaryData()) + xs, ys, fs = get_function_values(ctx,u,x_resolution=x_resolution,y_resolution=y_resolution,bdata=bdata) + dx, dy = gradient_from_values(xs,ys,fs) + + return xs, ys, dx, dy +end + +""" + get_levelset_evolution(ctx, u, u_dot; + x_resolution=nothing, y_resolution=nothing, bdata=BoundaryData()) + +Returns a vectorfield describing how levelsets of u change according to the derivative u_dot. +Mainly to be used for plotting. +""" +function get_levelset_evolution(ctx, u, u_dot; + x_resolution=nothing, y_resolution=nothing, bdata=BoundaryData()) + xs, ys, dx, dy = get_gradient_field(ctx, u, x_resolution=x_resolution, y_resolution=y_resolution,bdata=bdata) + xs_dot, ys_dot, fs_dot = get_function_values(ctx, u_dot, x_resolution=x_resolution, y_resolution=y_resolution,bdata=bdata) + + @assert xs==xs_dot + @assert ys==ys_dot + + norm_gradient = sqrt.(dx.*dx + dy.*dy) + + return xs, ys, -fs_dot.*dx./norm_gradient, -fs_dot.*dy./norm_gradient +end diff --git a/src/exports.jl b/src/exports.jl index fe11060f..cc035211 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -150,4 +150,10 @@ export #odesolvers.jl LinearImplicitEuler, - LinearMEBDF2 + LinearMEBDF2, + + #dynamicIsoperimetry.jl + dynamic_cheeger_value, + get_minimal_levelset, + get_levelset_evolution, + get_levelset From d9872d3bd5bec6e6c33a79e6acbd1dfee7ff7fd2 Mon Sep 17 00:00:00 2001 From: Martin Hermann Date: Thu, 22 Jul 2021 14:16:01 +0200 Subject: [PATCH 3/6] add basic linear response functionality --- src/exports.jl | 9 +- src/linearResponse.jl | 265 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 273 insertions(+), 1 deletion(-) create mode 100644 src/linearResponse.jl diff --git a/src/exports.jl b/src/exports.jl index fe11060f..ade40ce0 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -150,4 +150,11 @@ export #odesolvers.jl LinearImplicitEuler, - LinearMEBDF2 + LinearMEBDF2, + + #linearResponse.jl + linear_response_tensor, + linearized_flow_autodiff, + adaptiveTOCollocationLinearResponseMatrix, + nonadaptiveTOCollocationLinearResponseMatrix, + get_linear_response diff --git a/src/linearResponse.jl b/src/linearResponse.jl new file mode 100644 index 00000000..3cdf5b53 --- /dev/null +++ b/src/linearResponse.jl @@ -0,0 +1,265 @@ +# functions to be used in the context of the methods described in https://arxiv.org/abs/1907.10852 + + +""" + linear_response_tensor(parametrized_flow_map, u, p) + +Calculate the tensor dot_A from the linear response paper, when given a flow map T(x,p) +""" +function linear_response_tensor(parametrized_flow_map, u, p) + return linear_response_tensor( + x->parametrized_flow_map(x,p), + x->flow_map_parameter_derivative(parametrized_flow_map, x, p), + u, + p) +end + +""" + linear_response_tensor(flow_map, parameter_derivative, u, p) + +Calculate the tensor dot_A from the linear response paper, when given a flow map T(x) and its +parameter derivative dotT(x) +""" +function linear_response_tensor(flow_map, parameter_derivative, u, p) + DT = linearized_flow_autodiff(flow_map, u) + DTdot = linearized_flow_autodiff(parameter_derivative, u) + DTinv = inv(DT) + return -Tensors.symmetric(DTinv ⋅ DTdot ⋅ Tensors.dott(DTinv)) +end + +""" + linearized_flow_autodiff(flow_map, x) + +Calculate derivative of flow map by automatic differentiation + +Return the time-resolved base trajectory and its associated linearized flow maps. +""" +linearized_flow_autodiff(flow_map, x) = Tensor{2,2}(ForwardDiff.jacobian(flow_map, x)) + + +flow_map_parameter_derivative(flow_map, u, p) = ForwardDiff.jacobian(x -> flow_map(x[1:end - 1], x[end]), vcat(u, p))[:,end] + +""" + adaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, parametrized_flow_map, p; + quadrature_order=default_quadrature_order, + on_torus::Bool=false, + on_cylinder::Bool=false, + LL_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[1], + UR_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[2], + bdata::BoundaryData=BoundaryData(), + flow_map_mode=0 + ) + +Return the linear response matrix L computed by the adaptive TO method, using a parametrized flow map T(x,p) + +Also note that, as with all functions of the TO method, BoundaryData makes sense only for periodic +boundary data. If we have Dirichlet boundary conditions, it is best to simply ignore them here +and apply them to L afterwards. +""" +function adaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, parametrized_flow_map, p; + quadrature_order=default_quadrature_order, + on_torus::Bool=false, + on_cylinder::Bool=false, + LL_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[1], + UR_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[2], + bdata::BoundaryData=BoundaryData(), + flow_map_mode=0 + ) + return adaptiveTOCollocationLinearResponseMatrix( + ctx, + x->parametrized_flow_map(x,p), + x->flow_map_parameter_derivative(parametrized_flow_map, x, p), + p, + quadrature_order=quadrature_order,on_torus=on_torus,on_cylinder=on_cylinder,LL_future=LL_future, + UR_future=UR_future,bdata=bdata,flow_map_mode=flow_map_mode) +end + +""" + adaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, flow_map, parameter_derivative, p; + quadrature_order=default_quadrature_order, + on_torus::Bool=false, + on_cylinder::Bool=false, + LL_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[1], + UR_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[2], + bdata::BoundaryData=BoundaryData(), + flow_map_mode=0 + ) + +Return the linear response matrix L computed by the adaptive TO method, using a flow map T(x) and its +parameter derivative dotT(x) + +Also note that, as with all functions of the TO method, BoundaryData makes sense only for periodic +boundary data. If we have Dirichlet boundary conditions, it is best to simply ignore them here +and apply them to L afterwards. +""" +function adaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, flow_map, parameter_derivative, p; + quadrature_order=default_quadrature_order, + on_torus::Bool=false, + on_cylinder::Bool=false, + LL_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[1], + UR_future::Tuple{<:Real,<:Real}=ctx.spatialBounds[2], + bdata::BoundaryData=BoundaryData(), + flow_map_mode=0 + ) + + # adapted from adaptiveTOCollocationStiffnessMatrix + # functionality is slightly reduced from there, as we assume only the volume preserving case + # and a single time in the future + + flow_map_images = zeros(Vec{2}, ctx.n) + for i in 1:ctx.n + if flow_map_mode == 0 + flow_map_images[i] = Vec{2}(flow_map(ctx.grid.nodes[i].x)) + else + flow_map_images[i] = Vec{2}(flow_map(i)) + end + end + + flow_map_t(j) = flow_map_images[j] + + # we only support volume_preserving here for now + new_ctx,new_bdata,_ = adaptiveTOFutureGrid(ctx, flow_map_t; + on_torus=on_torus, + on_cylinder=on_cylinder, + LL_future=LL_future, + UR_future=UR_future, + bdata=bdata, + quadrature_order=quadrature_order, + flow_map_mode=1) + + # find the derivatives on the nodes of the old grid, then reorder this so it is in dof order for the new grid + derivative_on_old_grid = zeros(2,ctx.n) + for i in 1:ctx.n + derivative_on_old_grid[:,i] = parameter_derivative(ctx.grid.nodes[i].x) + end + # keep in mind: node order for new_ctx is bcdof order for ctx, which allows us to link the two grids + # unfortunately we need some intermediate steps. new_ctx potentially has more nodes than ctx (we need + # to duplicate some if we are e.g. on a torus), so new_ctx.dof_to_node might fail + translation_table_bcdof_new_to_node_old = bcdof_to_node(ctx,bdata)[bcdof_to_node(new_ctx,new_bdata)] + derivative_in_bcdof_order_new_grid = derivative_on_old_grid[:,translation_table_bcdof_new_to_node_old] + # we need the derivatives in dof order instead of in bcdof order + W = mapslices(x->undoBCS(new_ctx, x, new_bdata),derivative_in_bcdof_order_new_grid,dims=2) + + L = assembleLinearResponseMatrix(new_ctx, W, bdata=new_bdata) + + # change L from bcdof order for ctx_new to node order for ctx_new. This is by design the same + # as bcdof order for ctx, which is the order we need it in (this can be slightly confusing) + translation_table_bcdof_new_to_bcdof_old = bcdof_to_node(new_ctx,new_bdata) + I, J, V = Main.CoherentStructures.findnz(L) + I .= translation_table_bcdof_new_to_bcdof_old[I] + J .= translation_table_bcdof_new_to_bcdof_old[J] + n = ctx.n - length(bdata.periodic_dofs_from) + L = sparse(I,J,V,n,n) + + return 0.5(L+L') +end + + + +""" +function nonadaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, parametrized_flow_map, A, p; + bdata::BoundaryData=BoundaryData()) + +Return the linear response matrix L computed by the nonadaptive TO method, using a parametrized flow +map T(x,p) as well as a representation matrix A computed by nonAdaptiveTOCollocation + +Also note that, as with all functions of the TO method, BoundaryData makes sense only for periodic +boundary data. If we have Dirichlet boundary conditions, it is best to simply ignore them here +and apply them to L afterwards. +""" +function nonadaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, parametrized_flow_map, A, p; + bdata::BoundaryData=BoundaryData()) + return nonadaptiveTOCollocationLinearResponseMatrix( + ctx, + x->parametrized_flow_map(x,p), + x->flow_map_parameter_derivative(parametrized_flow_map, x, p), + A, + p; + bdata=bdata) +end + + +""" +function nonadaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, parametrized_flow_map, A, p; + bdata::BoundaryData=BoundaryData()) + +Return the linear response matrix L computed by the nonadaptive TO method, using a flow map T(x) and +its parameter derivative dotT(x) as well as a representation matrix A computed by nonAdaptiveTOCollocation + +Also note that, as with all functions of the TO method, BoundaryData makes sense only for periodic +boundary data. If we have Dirichlet boundary conditions, it is best to simply ignore them here +and apply them to L afterwards. +""" +function nonadaptiveTOCollocationLinearResponseMatrix(ctx::GridContext{2}, flow_map, parameter_derivative, A, p; + bdata::BoundaryData=BoundaryData()) + W = zeros(2,ctx.n) + for i in 1:ctx.n + W[:,i] = parameter_derivative(ctx.grid.nodes[i].x) + end + # from node order to bcdof order + W = W[:,bcdof_to_node(ctx, bdata)]; + WA = W*A' + # and now back to dof order + WA = mapslices(x->undoBCS(ctx, x, bdata),WA,dims=2) + + L = assembleLinearResponseMatrix(ctx, WA, bdata=bdata) + + return 0.5A'*(L+L')*A +end + + +# adapted from _assembleStiffnessMatrix +# WA is a matrix containing all the derivatives at the dofs, possibly multiplied by a representation +# matrix A in the nonadaptive case +function assembleLinearResponseMatrix(ctx, WA; bdata::BoundaryData=BoundaryData()) + cv = FEM.CellScalarValues(ctx.qr, ctx.ip, ctx.ip_geom) + dh = ctx.dh + K = FEM.create_sparsity_pattern(dh) + a_K = FEM.start_assemble(K) + dofs = zeros(Int, FEM.ndofs_per_cell(dh)) + n = FEM.getnbasefunctions(cv) # number of basis functions + Ke = zeros(n, n) + + index = 1 # quadrature point counter + + @inbounds for cell in FEM.CellIterator(dh) + fill!(Ke, 0) + FEM.reinit!(cv, cell) + for q in 1:FEM.getnquadpoints(cv) + dΩ = FEM.getdetJdV(cv, q) * ctx.mass_weights[index] + for i in 1:n + ∇φᵢ = FEM.shape_gradient(cv, q, i) + for j in 1:n + ∇φⱼ = FEM.shape_gradient(cv, q, j) + for k in 1:n + ∇φₖ = FEM.shape_gradient(cv, q, k) + Ke[i,j] += (∇φᵢ ⋅ WA[:,FEM.celldofs(cell)[k]]) * (∇φₖ ⋅ ∇φⱼ) * dΩ + end + end + end + index += 1 + end + FEM.celldofs!(dofs, cell) + FEM.assemble!(a_K, dofs, Ke) + end + return applyBCS(ctx, K, bdata) +end + +""" + get_linear_response(u₀,λ₀,M,K,L) + +Solve the linear system defined in the linear response paper, given an eigenpair (u₀,λ₀), mass matrix +M, stiffness matrix K and linear response matrix L +""" +function get_linear_response(u₀,λ₀,M,K,L) + lhs = [K-λ₀*M -M*u₀ + u₀'*M 0] + rhs = [-L*u₀ ; 0] + + sol = lhs \ rhs + + u_dot = sol[1:end-1] + λ_dot = sol[end] + + return u_dot, λ_dot +end From 8fd3206101f16efa11d13907d913ad787bfd2175 Mon Sep 17 00:00:00 2001 From: Martin Hermann Date: Tue, 27 Jul 2021 14:15:20 +0200 Subject: [PATCH 4/6] add minor updates (exports, documentation) --- docs/src/basics.md | 1 + docs/src/fem.md | 4 +++- src/CoherentStructures.jl | 8 ++++++++ 3 files changed, 12 insertions(+), 1 deletion(-) diff --git a/docs/src/basics.md b/docs/src/basics.md index 40e5c4f9..a9cb88ca 100644 --- a/docs/src/basics.md +++ b/docs/src/basics.md @@ -61,6 +61,7 @@ pullback_tensors pullback_metric_tensor pullback_diffusion_tensor pullback_SDE_diffusion_tensor +linear_response_tensor ``` A second-order symmetric two-dimensional tensor (field) may be diagonalized diff --git a/docs/src/fem.md b/docs/src/fem.md index cf06ebe4..374fb8ad 100644 --- a/docs/src/fem.md +++ b/docs/src/fem.md @@ -2,7 +2,7 @@ These methods rely on the theory outlined by Froyland's [*Dynamical Laplacian*](http://arxiv.org/pdf/1411.7186v4.pdf) -and the [*Geometric Heat Flow*](https://www.researchgate.net/publication/306291640_A_geometric_heat-flow_theory_of_Lagrangian_coherent_structures) of Karrasch & Keller. +and the [*Geometric Heat Flow*](https://www.researchgate.net/publication/306291640_A_geometric_heat-flow_theory_of_Lagrangian_coherent_structures) of Karrasch & Keller. For the [*Linear Response*](https://arxiv.org/abs/1907.10852), the main source is Antown et al. The Laplace-like operators are best discretized by finite-element-based methods, see this [paper](https://arxiv.org/pdf/1705.03640.pdf) by Froyland & Junge. @@ -188,6 +188,8 @@ CurrentModule = CoherentStructures assembleStiffnessMatrix assembleMassMatrix adaptiveTOCollocationStiffnessMatrix +adaptiveTOCollocationLinearResponseMatrix +nonadaptiveTOCollocationLinearResponseMatrix ``` ### Constructing Grids diff --git a/src/CoherentStructures.jl b/src/CoherentStructures.jl index 2aa6f0d8..8da21f26 100644 --- a/src/CoherentStructures.jl +++ b/src/CoherentStructures.jl @@ -46,6 +46,8 @@ const FEM = Ferrite #Other miscallaneous packages using RecipesBase import ColorTypes +using ForwardDiff +using Contour # contains a list of exported functions include("exports.jl") @@ -113,4 +115,10 @@ include("ulam.jl") # plotting functionality include("plotting.jl") +# linear response related methods +include("linearResponse.jl") + +# functions for dynamic isoperimetry +include("dynamicIsoperimetry.jl") + end From a073a3f56194bb69e95ae729be895a4c0b06f191 Mon Sep 17 00:00:00 2001 From: Martin Hermann Date: Tue, 27 Jul 2021 14:53:16 +0200 Subject: [PATCH 5/6] add example of linear response --- docs/examples/rot_double_gyre.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/docs/examples/rot_double_gyre.jl b/docs/examples/rot_double_gyre.jl index 0eccfec4..0a739943 100644 --- a/docs/examples/rot_double_gyre.jl +++ b/docs/examples/rot_double_gyre.jl @@ -90,6 +90,30 @@ res = Plots.plot([plot_u(ctx2, u[:,i], 200, 200, c=:viridis, cbar=:none) for i i DISPLAY_PLOT(res, rot_double_gyre_fem) +# In order to compute the sensitivity towards a parameter (e.g. the stoptime of the system), we can also use FEM methods: + +T(x,p) = flow(rot_double_gyre, x, [0.0, 1.0 + p], + tolerance = 1e-5)[end]; +Adot = x -> linear_response_tensor(T, x, 0) +L = assembleStiffnesMatrix(ctx, Adot) +v_dot, λ_dot = get_linear_response(v[:,1],λ[i],M,K,L) + +# Now `v_dot` is a linear approximation to the change in `v` when we increase or decrease the stoptime of the system. To illustrate this, let us compute the eigenvalues at a different time. + +ϵ = -0.2 +Aϵ = x -> mean_diff_tensor(rot_double_gyre, x, [0.0, 1.0 + ϵ], 1.e-10, tolerance= 1.e-4) +Kϵ = assembleStiffnessMatrix(ctx, Aϵ) +λϵ, vϵ = eigs(-Kϵ, M, which=:SM); + +# We can now compare `vϵ[:,1]` to `v[:,1] - ϵ*v_dot`: + +res = Plots.plot([plot_u(ctx, v[:,1], 100, 100, colorbar=:none, clim=(-3,3)), + plot_u(ctx, v_dot, 100, 100, colorbar=:none, clim=(-3,3)), + plot_u(ctx, vϵ[:,1], 100, 100, colorbar=:none, clim=(-3,3)), + plot_u(ctx, v[:,1] + ϵ*v_dot, 100, 100, colorbar=:none, clim=(-3,3))]) + +DISPLAY_PLOT(res, rot_double_gyre_linear_response) + # ## Geodesic vortices # Here, we demonstrate how to calculate black-hole vortices, see From 41b4802c7e4f09c2c19672712c26c410c1c1bebd Mon Sep 17 00:00:00 2001 From: Martin Hermann Date: Tue, 27 Jul 2021 15:14:36 +0200 Subject: [PATCH 6/6] adjust to the fact that the first ev is trivial --- docs/examples/rot_double_gyre.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/examples/rot_double_gyre.jl b/docs/examples/rot_double_gyre.jl index 0a739943..5d838558 100644 --- a/docs/examples/rot_double_gyre.jl +++ b/docs/examples/rot_double_gyre.jl @@ -96,7 +96,7 @@ T(x,p) = flow(rot_double_gyre, x, [0.0, 1.0 + p], tolerance = 1e-5)[end]; Adot = x -> linear_response_tensor(T, x, 0) L = assembleStiffnesMatrix(ctx, Adot) -v_dot, λ_dot = get_linear_response(v[:,1],λ[i],M,K,L) +v_dot, λ_dot = get_linear_response(v[:,2],λ[i],M,K,L) # Now `v_dot` is a linear approximation to the change in `v` when we increase or decrease the stoptime of the system. To illustrate this, let us compute the eigenvalues at a different time. @@ -105,12 +105,12 @@ Aϵ = x -> mean_diff_tensor(rot_double_gyre, x, [0.0, 1.0 + ϵ], 1.e-10, tolera Kϵ = assembleStiffnessMatrix(ctx, Aϵ) λϵ, vϵ = eigs(-Kϵ, M, which=:SM); -# We can now compare `vϵ[:,1]` to `v[:,1] - ϵ*v_dot`: +# We can now compare `vϵ[:,2]` to `v[:,2] - ϵ*v_dot`: -res = Plots.plot([plot_u(ctx, v[:,1], 100, 100, colorbar=:none, clim=(-3,3)), +res = Plots.plot([plot_u(ctx, v[:,2], 100, 100, colorbar=:none, clim=(-3,3)), plot_u(ctx, v_dot, 100, 100, colorbar=:none, clim=(-3,3)), - plot_u(ctx, vϵ[:,1], 100, 100, colorbar=:none, clim=(-3,3)), - plot_u(ctx, v[:,1] + ϵ*v_dot, 100, 100, colorbar=:none, clim=(-3,3))]) + plot_u(ctx, vϵ[:,2], 100, 100, colorbar=:none, clim=(-3,3)), + plot_u(ctx, v[:,2] + ϵ*v_dot, 100, 100, colorbar=:none, clim=(-3,3))]) DISPLAY_PLOT(res, rot_double_gyre_linear_response)