Update two region example

examples/lid_two_regions/.gitignore

+inputs/*

examples/lid_two_regions/input.yml

@@ -53,7 +53,7 @@ stabilization:
 # - body_force: Body force for mode==stress||r13
 nsd: 2
 mode: r13
-heat_source: 0
+heat_source: 1
 mass_source: 0
 body_force: [0,0]
 
@@ -82,7 +82,7 @@ bcs:
   3000: # upper slip wall
     chi_tilde: 1.0
     theta_w: 1
-    u_t_w: -1
+    u_t_w: 0
     u_n_w: 0
     p_w: 0
     epsilon_w: 0
@@ -116,8 +116,8 @@ convergence_study:
 # - write_vecs: Write all solution fields as vectors
 # - massflow: List of BC IDs to compute massflow J=int_bc dot(u,n) ds
 postprocessing:
-  write_pdfs: True
-  write_vecs: False
+  write_pdfs: False
+  write_vecs: True
   massflow: []
 
 # Parameter Study
@@ -126,6 +126,6 @@ postprocessing:
 # - parameter_key: Key as list, e.g. ["elemenets", "p", "degree"]
 # - parameter_values: List of value for parameter, e.g. [0.01,0.1,1,10]
 parameter_study:
-  enable: False
-  parameter_key: []
-  parameter_values: []
+  enable: True
+  parameter_key: ["regs", 4100, "kn"]
+  parameter_values: [0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]

examples/lid_two_regions/process-study.py

+import numpy
+
+dofs = 1425
+kn_range = [(i + 1) / 10 for i in range(10)]
+numknudsen = len(kn_range)
+fields = [
+    "theta", "s_x", "s_y", "p", "u_x", "u_x", "sigma_xx", "sigma_xy", "sigma_yx"
+]
+numfields = len(fields)
+
+V = numpy.zeros(shape=(numknudsen**2, numfields * dofs))
+
+for iKn, kn1 in enumerate(kn_range):
+    for jKn, kn2 in enumerate(kn_range):
+        # row = numpy.zeros(shape=(numfields * dofs))
+
+        field_data = numpy.zeros(shape=(numfields, dofs))
+        for i, field in enumerate(fields):
+            data = numpy.loadtxt("{}_{}/{}_0.mat".format(kn1, kn2, field))
+            field_data[i] = data
+            # print(field_data)
+            # print(data)
+            # print(len(data))
+            # print(data.shape)
+            # print("next")
+        row = field_data.reshape(numfields * dofs)
+        # print(row)
+        print(iKn * numknudsen + jKn)
+        V[iKn * numknudsen + jKn] = row
+
+numpy.savetxt("V_{}.mat".format(numknudsen**2), V)
+print(V)

examples/lid_two_regions/run-study.py

+import yaml
+import os
+import subprocess
+
+
+# setup
+input_file = "input.yml"
+output_folder = "inputs/"
+
+# setup parameters
+kn_range = [(i + 1) / 10 for i in range(10)]
+
+# read template file
+try:
+    with open(input_file, "r") as stream:
+        data = yaml.safe_load(stream)
+except FileNotFoundError:
+    print(f"[Error] The input file {input_file} is not found.\n")
+    raise
+
+# create output folder if not exists
+os.makedirs(os.path.dirname(output_folder), exist_ok=True)
+
+# write all input_files
+for kn1 in kn_range:
+    data["regs"][4000]["kn"] = kn1
+    data["output_folder"] = "{}_".format(kn1)
+
+    output_path = output_folder + "{}.yml".format(kn1)
+    with open(output_path, "w") as stream:
+        yaml.dump(data, stream)
+
+# run study
+for kn1 in kn_range:
+    subprocess.check_call([
+        "fenicsR13", output_folder + "{}.yml".format(kn1)
+    ], cwd=".")

examples/lid_two_regions/svd-experiment.jl

+using DelimitedFiles
+using LinearAlgebra
+using Plots
+using SymPy
+
+kn_range = 0.1:0.1:1.0
+
+V100 = readdlm("V_100.mat", ' ', Float64, '\n');
+F100 = svd(V100);
+F100SN = normalize(F100.S, Inf)
+
+# low rank approximation with rank R
+R = 20
+V100LR = F100.U[:,1:R] * diagm(F100.S[1:R]) * F100.V[:,1:R]'
+
+# sanity check: Frobenius norm should be S[R+1] is not the case but still low
+# ! FIXME
+println(norm(V100 - V100LR ,2))
+
+L = F100.U[:,1:R]
+OPCompl = I(100) - L * L'  # projector to complement of LL'
+
+# ansatz
+# p(t1, t2, deg) = t1^deg  # monomials
+# p(t1, t2, deg) = (1/t1)^deg  # 1/x
+# p(t1, t2, deg) = exp(t1*deg)  # exponentials
+# p(t1, t2, deg) = log(t1)^deg  # logrithms
+# p(t1, t2, deg) = (t1*t2)^deg  # logrithms
+# p(t1, t2, deg) = sin(deg * pi * ((t2-0.1)*100+t1*10) / 100)  # logrithms
+# p(t1, t2, deg) = (t2-0.1)*100+t1*10  # logrithms
+# p(t1, t2, deg) = 1/(t1+t2)
+# p(t1, t2, deg) = t2 + 1/(t1)^deg
+p(t1, t2, deg) = 1/(t2)^deg + 1/(t1)^deg
+# function p(t1, t2, deg)
+#   x = Sym("x")
+#   f(x) = x^float(deg-1)
+#   G(y) = subs(integrate(f(x), x), x, y)
+#   return N(G(float(t1)))
+# end
+
+degrange = -3:3
+P = vcat(
+  [[p(kn1, kn2, d) for d in degrange]' for kn1 in kn_range, kn2 in kn_range]...,
+)
+Q = Matrix(qr(P).Q)
+objective = OPCompl * Q  # --> MIN!!!
+
+display(plot([svd(OPCompl).S, svd(objective).S], yaxis=:log, title="Singular Values, max ansatz degree, R=$R", xlabel="index", ylabel="value", label=["(I - LL')" "(I - LL') * Q"], marker = (:dot)))
+
+
+
+if false # comparison
+  V361 = readdlm("V_361.mat", ' ', Float64, '\n');
+  F361 = svd(V361);
+  F361SN = normalize(F361.S, Inf)
+
+  V400 = readdlm("V_400.mat", ' ', Float64, '\n');
+  F400 = svd(V400);
+  F400SN = normalize(F400.S, Inf)
+
+  max = 100
+  plot([F100SN, F400SN[1:max]], yaxis=:log, title="Normlized Singular Values", xlabel="index", ylabel="value", label=["h=0.1" "h=0.05"], marker = (:dot))
+end
+
+

examples/lid_two_regions/svd-experiment_2.jl

+using DelimitedFiles
+using LinearAlgebra
+using Plots
+
+kn_range = 0.1:0.1:1.0
+
+V100 = readdlm("V_100.mat", ' ', Float64, '\n');
+F100 = svd(V100);
+F100SN = normalize(F100.S, Inf)
+
+svArray = []
+
+for R=10:10:100
+  # low rank approximation with rank R
+  # R = 50
+  V100LR = F100.U[:,1:R] * diagm(F100.S[1:R]) * F100.V[:,1:R]'
+
+  # sanity check: Frobenius norm should be S[R+1] is not the case but still low
+  # ! FIXME
+  println(norm(V100 - V100LR ,2))
+
+  L = F100.U[:,1:R]
+  OPCompl = I(100) - L * L'  # projector to complement of LL'
+
+  # ansatz
+  # p(t1, t2, deg) = t1^deg  # monomials
+  p(t1, t2, deg) = (1/t1)^deg  # monomials
+  maxdeg = 1
+  P = vcat(
+    [[p(kn1, kn2, d) for d=1:1]' for kn1 in kn_range, kn2 in kn_range]...
+  )
+  Q = Matrix(qr(P).Q)
+  objective = OPCompl * Q  # --> MIN!!!
+
+  push!(svArray, svd(objective).S...)
+end
+
+plot([svd(OPCompl).S, svArray], yaxis=:log, title="Singular Values, max ansatz degree $maxdeg, R=$R", xlabel="index", ylabel="value", label=["(I - LL')" "rank=10,20,...,100"], marker = (:dot))
+
+
+
+if false # comparison
+  V361 = readdlm("V_361.mat", ' ', Float64, '\n');
+  F361 = svd(V361);
+  F361SN = normalize(F361.S, Inf)
+
+  V400 = readdlm("V_400.mat", ' ', Float64, '\n');
+  F400 = svd(V400);
+  F400SN = normalize(F400.S, Inf)
+
+  max = 100
+  plot([F100SN, F400SN[1:max]], yaxis=:log, title="Normlized Singular Values", xlabel="index", ylabel="value", label=["h=0.1" "h=0.05"], marker = (:dot))
+end
+
+