more experiments

jehicken · Dec 6, 2023 · 04fee94 · 04fee94
1 parent 4efc1cd
commit 04fee94
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 9 deletions.
diff --git a/src/norm.jl b/src/norm.jl
@@ -427,12 +427,21 @@ function penalty(root::Cell{Data, Dim, T, L}, xc, xc_init, dist_ref, mu, degree
         for d = 1:Dim 
             dist += (xc[d,i] - xc_init[d,i])^2
         end 
-        phi += dist/(dist_ref[i]^2)
+        phi += 0.5*mu*dist/(dist_ref[i]^2)
     end
 
     # compute the diagonal norm based on x 
     H = diagonal_norm(root, xc, degree)
 
+    # # compute the discrete KS function
+    # rho = 100.0
+    # minH = minimum(real.(H))     
+    # sum_exp = 0.0
+    # for i = 1:num_nodes 
+    #     sum_exp += exp(rho*(minH - H[i]))
+    # end 
+    # phi += -minH + log(sum_exp/num_nodes)/rho
+
     # add the penalties 
     tol = 1e-5
     for i = 1:num_nodes 
@@ -443,7 +452,6 @@ function penalty(root::Cell{Data, Dim, T, L}, xc, xc_init, dist_ref, mu, degree
         end 
     end
 
-    phi *= 0.5
     return phi
 end
 
@@ -454,11 +462,31 @@ function penalty_grad!(g::AbstractVector{T}, root::Cell{Data, Dim, T, L},
     # need the diagonal norm for the reverse sweep 
     H = diagonal_norm(root, xc, degree)
 
+    # rho = 100.0
+    # minH = minimum(real.(H))
+    # sum_exp = 0.0
+    # for i = 1:num_nodes 
+    #     sum_exp += exp(rho*(minH - H[i]))
+    # end 
+
     # start the reverse sweep 
     fill!(g, zero(T))
+
     # return phi
     # phi *= 0.5 
-    phi_bar = 0.5 
+    # phi_bar = 0.5 
+
+    # return phi 
+    phi_bar = 1.0 
+
+    # # phi += -minH + log(sum_exp/num_nodes)/rho
+    # sum_exp_bar = phi_bar/(sum_exp * rho)
+    # H_bar = zero(H)
+    # for i = 1:num_nodes 
+    #     # sum_exp += exp(rho*(minH - H[i]))
+    #     H_bar[i] -= sum_exp_bar*exp(rho*(minH - H[i]))*rho
+    # end
+
     # add the penalties
     tol = 1e-5
     H_bar = zero(H)
@@ -478,8 +506,8 @@ function penalty_grad!(g::AbstractVector{T}, root::Cell{Data, Dim, T, L},
 
     # compute the norm part of the penalty
     for i = 1:num_nodes 
-        # phi += dist/(dist_ref[i]^2)
-        dist_bar = phi_bar/(dist_ref[i]^2)
+        # phi += 0.5*mu*dist/(dist_ref[i]^2)
+        dist_bar = 0.5*mu*phi_bar/(dist_ref[i]^2)
         for d = 1:Dim 
             # dist += (xc[d,i] - xc_init[d,i])^2
             xc_bar[d,i] += dist_bar*2.0*(xc[d,i] - xc_init[d,i])

diff --git a/test/opt-norm.jl b/test/opt-norm.jl
@@ -14,7 +14,7 @@ using SparseArrays
 #Random.seed!(42)
 
 Dim = 2
-degree = 5
+degree = 2
 
 # use a unit HyperRectangle 
 root = Cell(SVector(ntuple(i -> 0.0, Dim)),
@@ -24,7 +24,7 @@ root = Cell(SVector(ntuple(i -> 0.0, Dim)),
 # DGD dof locations
 num_basis = binomial(Dim + degree, Dim)
 
-num_nodes = 10*num_basis
+num_nodes = 100*num_basis
 points = rand(Dim, num_nodes)
 
 # refine mesh, build sentencil, and evaluate particular quad and nullspace
@@ -46,7 +46,7 @@ CutDGD.build_nn_stencils!(root, points, degree)
 #end
 #println(dist_ref)
 dist_ref = ones(num_nodes)
-mu = 100.0
+mu = 0.0001
 
 g = zeros(num_nodes*Dim)
 
@@ -71,7 +71,7 @@ for n = 1:max_iter
 
     obj0 = obj
     dxc = reshape(g, (Dim, num_nodes))
-    alpha = 1.0
+    alpha = 10.0
     for k = 1:max_line
         global points -= alpha*dxc
         obj = CutDGD.penalty(root, points, points_init, dist_ref, mu, degree)