Skip to content
Snippets Groups Projects
Select Git revision
  • 30a7932b63df1abf895e74c5fdcaefeeea482fe5
  • master default protected
  • releases
  • releases/3.0.3
4 results

check-package.m4

Blame
  • run.jl 3.77 KiB
    export myrun
    
    using LinearAlgebra: norm
    
    
    
    function myrun()
        name = "test"
        mesh = init_grid(50)
        d = 2
        m = 1
    
        Vg = FeSpace(mesh, DP0(), (1,))
        Vu = FeSpace(mesh, P1(), (m,))
        Vp1 = FeSpace(mesh, DP0(), (1,))
        Vp2 = FeSpace(mesh, DP1(), (m, d))
    
        g = FeFunction(Vg, name="g")
        mask = FeFunction(Vg, name="mask")
        u = FeFunction(Vu, name="u")
        p1 = FeFunction(Vp1, name="p1")
        p2 = FeFunction(Vp2, name="p2")
        du = FeFunction(Vu)
        dp1 = FeFunction(Vp1)
        dp2 = FeFunction(Vp2)
        nablau = nabla(u)
        nabladu = nabla(du)
    
        g.data .= 0
        u.data .= 0
        p1.data .= 0
        p2.data .= 0
        du.data .= 0
        dp1.data .= 0
        dp2.data .= 0
    
        alpha1 = 0.
        alpha2 = 10.
        beta = 0.
        lambda = 0.01
        gamma1 = 1e-3
        gamma2 = 1e-3
    
        interpolate!(g, x -> norm(x - SA[0.5, 0.5]) < 0.3)
        interpolate!(mask, x -> abs(x[2] - 0.5) > 0.1)
    
        #T(data, u) = u
        #tdata = u # FIXME
    
        @inline T(data, u) = iszero(data) ? zero(u) : u
        tdata = mask
    
        #T(data, u) = dot(data, u)
        #tdata = nabla(fw)
    
        S(u) = u
    
        function dp1_update(x; g, u, p1, du, tdata)
    	m1 = max(gamma1, norm(T(tdata, u) - g))
    	cond = norm(T(tdata, u) - g) > gamma1 ?
    	    dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
    	    zeros(size(p1))
    	return -p1 + alpha1 / m1 * (T(tdata, u) + T(tdata, du) - g) - cond
        end
    
    
        function dp2_update(x; u, nablau, p2, du, nabladu)
    	m2 = max(gamma2, norm(nablau))
    	cond = norm(nablau) > gamma2 ?
    	    dot(nablau, nabladu) / norm(nablau)^2 * p2 :
    	    zeros(size(p2))
    	return -p2 + lambda / m2 * (nablau + nabladu) - cond
        end
    
        function p1_project!()
    	p1.space.element::DP0
    	p1.data .= clamp.(p1.data, -alpha1, alpha1)
        end
    
        function p2_project!()
    	p2.space.element::DP1
    	p2d = reshape(p2.data, m * d * 3, :) # no copy
    	for i in axes(p2d, 2)
    	    p2in = norm(p2d[:, i])
    	    if p2in > lambda
    		p2d[:, i] .*= lambda ./ p2in
    	    end
    	end
        end
    
        @inline function du_a(x, du, nabladu, phi, nablaphi; g, u, nablau, p1, p2, tdata)
    	m1 = max(gamma1, norm(T(tdata, u) - g))
    	cond1 = norm(T(tdata, u) - g) > gamma1 ?
    	    dot(T(tdata, u) - g, T(tdata, du)) / norm(T(tdata, u) - g)^2 * p1 :
    	    zeros(size(p1))
    	a1 = alpha1 / m1 * dot(T(tdata, du), T(tdata, phi)) -
    	    dot(cond1, T(tdata, phi))
    
    	m2 = max(gamma2, norm(nablau))
    	cond2 = norm(nablau) > gamma2 ?
    	    dot(nablau, nabladu) / norm(nablau)^2 * p2 :
    	    zeros(size(p2))
    	a2 = lambda / m2 * dot(nabladu, nablaphi) -
    	    dot(cond2, nablaphi)
    
    	aB = alpha2 * dot(T(tdata, du), T(tdata, phi)) +
    	    beta * dot(S(du), S(phi))
    
    	return a1 + a2 + aB
        end
    
        @inline function du_l(x, phi, nablaphi; g, u, nablau, tdata)
    	aB = alpha2 * dot(T(tdata, u), T(tdata, phi)) +
    	    beta * dot(S(u), S(phi))
    	m1 = max(gamma1, norm(T(tdata, u) - g))
    	p1part = alpha1 / m1 * dot(T(tdata, u) - g, T(tdata, phi))
    	m2 = max(gamma2, norm(nablau))
    	p2part = lambda / m2 * dot(nablau, nablaphi)
    	gpart = alpha2 * dot(g, T(tdata, phi))
    
    	return -aB - p1part - p2part + gpart
        end
    
        function save_step!(pvd, i)
    	vtk = vtk_mesh("$(name)_$(lpad(i, 5, '0')).vtu", mesh)
    	vtk_append!(vtk, g, u, p1, p2)
    	vtk_save(vtk)
    	pvd[i] = vtk
        end
    
        pvd = paraview_collection("$(name).pvd")
        save_step!(pvd, 0)
        for i = 1:20
    	print("[$(lpad(i, 3, '0'))] ")
    	# solve du
    	print("assemble ... ")
    	A = assemble(Vu, du_a; g, u, nablau, p1, p2, tdata)
    	b = assemble_rhs(Vu, du_l; g, u, nablau, tdata)
    	print("solve ... ")
    	du.data .= A \ b
    
    	# solve dp1, dp2
    	interpolate!(dp1, dp1_update; g, u, p1, du, tdata)
    	interpolate!(dp2, dp2_update; u, nablau, p2, du, nabladu)
    
    	# newton update
    	u.data .+= du.data
    	p1.data .+= dp1.data
    	p2.data .+= dp2.data
    
    	# reproject p1, p2
    	p1_project!()
    	p2_project!()
    
    	print("save ... ")
    	save_step!(pvd, i)
    	println()
        end
        vtk_save(pvd)
        return
    end