Adding Fully Coupled Non-Linear Case

Author: carlodev

examples/TGVNlin.jl

@@ -0,0 +1,61 @@
+using PartitionedArrays
+using SegregatedVMSSolver
+using SegregatedVMSSolver.ParametersDef
+using SegregatedVMSSolver.SolverOptions
+
+using MPI
+
+
+
+
+
+t0 = 0.0
+dt = 0.01
+tF = 2.5
+vortex_diameter = 1.0
+N = 16
+Re = 1600
+D = 2
+
+backend = with_debug
+
+rank_partition = (2,2)
+
+
+
+
+
+solver_options =  options = "-snes_type newtonls -snes_linesearch_type basic -snes_linesearch_damping 1.0 -snes_rtol 1.0e-8 -snes_atol 0 -snes_monitor  -snes_max_it 20 \
+-pc_type asm -ksp_type fgmres -ksp_converged_reason -ksp_max_it 50 -ksp_rtol 1e-8 -ksp_atol 0.0"
+
+sprob = StabilizedProblem(VMS(1))
+timep = TimeParameters(t0=t0, dt=dt, tF=tF)
+
+physicalp = PhysicalParameters(Re=Re, c=vortex_diameter)
+solverp = SolverParameters(matrix_freq_update=1, Number_Skip_Expansion=10e6, M=40,
+petsc_options=solver_options, linear=false, θ=1.0)
+exportp = ExportParameters(printinitial=false, printmodel=false, 
+vtu_export = ["uh","ph","uh_analytic", "ph_analytic"], extra_export=["VelocityError","PressureError"])
+
+
+
+
+meshp = MeshParameters(rank_partition, D; N=N, L=0.5 * vortex_diameter)
+simparams = SimulationParameters(timep, physicalp, solverp, exportp)
+bc_tgv = Periodic(meshp, physicalp)
+
+
+
+mcase = TaylorGreen(bc_tgv, meshp, simparams, sprob)
+
+
+
+# Create folder and file
+# mkdir("Log")
+# open("Log/PrintSim.txt", "w") do file
+# end
+
+
+
+SegregatedVMSSolver.solve(mcase, backend)
+

src/Commons/Equations/CoupledEquations.jl

@@ -0,0 +1,79 @@
+function coupled_equations(params::Dict{Symbol,Any},simcase::SimulationCase)
+    @sunpack skew, ν,dt, θ, D = simcase
+    
+    sprob = simcase.sprob
+    @unpack dΩ = params
+    @unpack skew = sprob
+
+    #Conservations equations
+    Rc(u) = ∇ ⋅ u
+    dRc(du) = ∇ ⋅ du
+    Rm(t, (u, p)) = ∂t(u) + u ⋅ ∇(u) + ∇(p)  #- ν*Δ(u)
+    dRm((u, p), (du, dp), (v, q)) = du ⋅ ∇(u) + u ⋅ ∇(du) + ∇(dp) #- ν*Δ(du)
+    stab_coeff = compute_stab_coeff(simcase,params)
+
+    Tm(u) = momentum_stabilization(u, stab_coeff, simcase)
+    Tc(u) = continuity_stabilization(u, stab_coeff, simcase)
+  
+
+    θvp = 1.0
+
+    #VMS equations
+    #Function used in VMS terms
+    TRm(t, (u, p)) = Tm(u.cellfield) * Rm(t, (u, p))
+  
+    #Variational equations
+    Bᴳ(t, (u, p), (v, q)) = ∫(∂t(u) ⋅ v)dΩ + ∫((u ⋅ ∇(u)) ⋅ v)dΩ - ∫((∇ ⋅ v) * p)dΩ + ∫((q * (∇ ⋅ u)))dΩ + ν * ∫(∇(v) ⊙ ∇(u))dΩ 
+  
+    #SUPG terms 
+    B_SUPG(t, (u, p), (v, q)) = ∫((u ⋅ ∇(v) + ∇(q)) ⊙ TRm(t, (u, p)))dΩ + ∫((∇ ⋅ v) ⊙ (Tc(u.cellfield)  * Rc(u)))dΩ
+  
+    #First VMS term
+    B_VMS1(t, (u, p), (v, q)) = ∫((u ⋅ (∇(v))') ⊙ TRm(t, (u, p)))dΩ
+  
+    #Second VMS term
+    B_VMS2(t, (u, p), (v, q)) = -1 * ∫((∇(v)) ⊙ (outer(TRm(t, (u, p)), TRm(t, (u, p)))))dΩ
+  
+    #Adding all the contributions
+    Bᴹ(t, (u, p), (v, q)) = Bᴳ(t, (u, p), (v, q)) + B_SUPG(t, (u, p), (v, q)) + B_VMS1(t, (u, p), (v, q)) + B_VMS2(t, (u, p), (v, q))
+  
+    res(t, (u, p), (v, q)) = Bᴹ(t, (u, p), (v, q))
+  
+    #VMS Jacobian
+    #Function used in the VMS terms -  derivative
+    dTRm((u, p), (du, dp), (v, q)) =  Tm(u.cellfield)* dRm((u, p), (du, dp), (v, q)) #Derivative of the Function used in the second VMS term
+  
+    #Variational equations - derivative
+    dBᴳ(t, (u, p), (du, dp), (v, q)) = ∫(((du ⋅ ∇(u)) ⋅ v) + ((u ⋅ ∇(du)) ⋅ v) + (∇(dp) ⋅ v) + (q * (∇ ⋅ du)))dΩ + ν * ∫(∇(v) ⊙ ∇(du))dΩ
+  
+    #SUPG terms derivative
+    dB_SUPG(t, (u, p), (du, dp), (v, q)) = ∫(((du ⋅ ∇(v)) ⋅ TRm(t, (u, p))) + ((u ⋅ ∇(v) + ∇(q)) ⋅ dTRm((u, p), (du, dp), (v, q))) + ((∇ ⋅ v) ⋅ ( Tc(u.cellfield)  .* dRc(du))))dΩ
+  
+    #First VMS term derivative
+    dB_VMS1(t, (u, p), (du, dp), (v, q)) = ∫(((du ⋅ ∇(v)') ⊙ TRm(t, (u, p))) + ((u ⋅ ∇(v)') ⊙ dTRm((u, p), (du, dp), (v, q))))dΩ
+  
+  
+    #Second VMS term derivative
+    dB_VMS2(t, (u, p), (du, dp), (v, q)) = ∫(∇(v) ⊙ (outer(dTRm((u, p), (du, dp), (v, q)), TRm(t, (u, p)))))dΩ + ∫(∇(v) ⊙ (outer(TRm(t, (u, p)), dTRm((u, p), (du, dp), (v, q)))))dΩ
+  
+  
+    #Adding all the contributions
+    jac(t, (u, p), (du, dp), (v, q)) = dBᴳ(t, (u, p), (du, dp), (v, q)) + dB_SUPG(t, (u, p), (du, dp), (v, q)) + dB_VMS1(t, (u, p), (du, dp), (v, q)) - dB_VMS2(t, (u, p), (du, dp), (v, q))
+  
+    #VMS time-Jacobian
+    dtTRm(t, (u, p), (dut, dpt), (v, q)) = ( Tm(u.cellfield)) ⋅ dut #Derivative of the Function used in the second VMS term
+  
+    dtBᴳ(t, (u, p), (dut, dpt), (v, q)) = ∫(dut ⋅ v)dΩ 
+    dtB_SUPG(t, (u, p), (dut, dpt), (v, q)) =  ∫((u ⋅ ∇(v) + (θvp)*∇(q)) ⊙ dtTRm(t, (u, p), (dut, dpt), (v, q)))dΩ 
+    dtB_VMS1(t, (u, p), (dut, dpt), (v, q)) = ∫((u ⋅ (∇(v))') ⊙ dtTRm(t, (u, p), (dut, dpt), (v, q)))dΩ
+    dtB_VMS2(t, (u, p), (dut, dpt), (v, q)) =  ∫((∇(v)) ⊙ (outer(dtTRm(t, (u, p), (dut, dpt), (v, q)), TRm(t, (u, p)))))dΩ + ∫((∇(v)) ⊙ (outer(TRm(t, (u, p)), dtTRm(t, (u, p), (dut, dpt), (v, q)))))dΩ        
+    
+    jac_t(t, (u, p), (dut, dpt), (v, q)) = dtBᴳ(t, (u, p), (dut, dpt), (v, q)) +  dtB_SUPG(t, (u, p), (dut, dpt), (v, q)) + dtB_VMS1(t, (u, p), (dut, dpt), (v, q)) - dtB_VMS2(t, (u, p), (dut, dpt), (v, q))
+  
+    res, jac, jac_t
+  
+  end
+  
+      
+  
+  

src/Commons/Equations/Equations.jl

@@ -14,11 +14,12 @@ using GridapDistributed: Algebra
 using SegregatedVMSSolver.ParametersDef
 
 export segregated_equations
-
+export coupled_equations
 
 include("StabilizationParameters.jl")
 include("StabilizationOperations.jl")
 include("StabilizedEquations.jl")
+include("CoupledEquations.jl")
 
 
 end

src/Commons/Equations/StabilizationOperations.jl

@@ -45,7 +45,7 @@ function momentum_stabilization(uu, stab_coeff::TensorStabilization,simcase::Sim
         τ₂ = uu_new ⋅ G ⋅ uu_new
         return (τ₁ .+ τ₂ .+ τ₃) .^ (-1 / 2)
     end
-
+ 
     return τm ∘ (uu, G, GG)
 
 

src/Commons/Equations/StabilizedEquations.jl

@@ -62,5 +62,4 @@ function is_VMS(method::SUPG)
   return 0
 end
 
-    
 

src/Commons/ParametersDef/Params.jl

@@ -81,6 +81,7 @@ end
     a_err_threshold::Int64=200 ### threshold to be satisfied:: norm(accelation_initial)/norm(acceleration_final) > a_err_threshold
     M::Int64=20 ## maximum number of internal iterations
     Number_Skip_Expansion::Int64=100 #number of initial time steps where is not used the Taylor-Expansion to compute the velocity field at the next time step, but simply the velcity field at the previous one.
+    linear::Bool = true
 end
 
 

src/Commons/SolveNlin.jl

@@ -0,0 +1,80 @@
+module SolveNlin
+using Gridap
+using GridapDistributed
+using GridapPETSc
+using SparseArrays
+using PartitionedArrays
+using Parameters
+using Gridap.FESpaces
+
+using SegregatedVMSSolver.ParametersDef
+using SegregatedVMSSolver.SolverOptions
+using SegregatedVMSSolver.CreateProblem
+using SegregatedVMSSolver.MatrixCreation
+using SegregatedVMSSolver.VectorsOperations
+using SegregatedVMSSolver.ExportUtility
+using SegregatedVMSSolver.Interfaces
+using SegregatedVMSSolver.Equations
+
+export solve_nlin
+
+function intialize_fields(U,P,V,Q)
+    Y = MultiFieldFESpace([V, Q])
+    X = TransientMultiFieldFESpace([U, P])
+    return X,Y
+end
+
+function initialize_solve(simcase::SimulationCase,params::Dict{Symbol,Any})
+        @unpack trials,tests = params
+   
+        @sunpack t0,dt,save_sim_dir = simcase
+        
+        U,P = trials
+        V,Q = tests
+
+        X,Y = intialize_fields(U,P,V,Q)
+
+        Ut0 = U(t0)
+        Pt0 = P(t0)
+      
+        Ut0_1 = U(t0+dt)
+        Pt0_1 = P(t0+dt)
+      
+        merge!(params, Dict(:Utn => Ut0, :Ptn => Pt0, :Utn1 => Ut0_1, :Ptn1 => Pt0_1))
+      
+      
+        uh0, ph0 = create_initial_conditions(simcase,params)
+        xh0 = interpolate_everywhere([uh0, ph0], X(t0))
+
+
+    return xh0, X,Y
+
+end
+
+function solve_nlin(simcase::SimulationCase,params::Dict{Symbol,Any})
+    xh0, X,Y = initialize_solve(simcase,params)
+    
+    res, jac, jac_t = coupled_equations(params,simcase)
+
+    op = TransientFEOperator(res, jac, jac_t,X,Y)
+    @sunpack t0,dt,tF,petsc_options,θ = simcase
+
+        GridapPETSc.with(args=split(petsc_options)) do
+
+        nls = PETScNonlinearSolver()
+
+
+        ode_solver = ThetaMethod(nls,dt,θ)
+
+        nl_solution = solve(ode_solver, op,xh0, t0, tF)
+        
+        for ((uh_tn,ph_tn), tn) in nl_solution
+            println("Solution at time $(tn) computed")
+    
+        end
+
+    end #end GridapPETSc
+
+end
+
+end #end module

src/Main.jl

@@ -2,6 +2,7 @@ using MPI
 using SegregatedVMSSolver.ParametersDef
 using SegregatedVMSSolver.CreateProblem
 using SegregatedVMSSolver.SolveProblem
+using SegregatedVMSSolver.SolveNlin
 using PartitionedArrays, Gridap, GridapDistributed
 
 function solve(simcase::SimulationCase,backend::Function)
@@ -20,7 +21,12 @@ function solve(simcase::SimulationCase,backend::Function)
 
         printstructure(simcase)
         params = setup_case(simcase,distribute)
-        solve_case(params,simcase)
+        if simcase.simulationp.solverp.linear
+            solve_case(params,simcase)
+        else
+            solve_nlin(simcase,params)
+        end
+
     end
 
 return true

src/SegregatedVMSSolver.jl

@@ -17,6 +17,7 @@ include(joinpath("Commons","MatrixCreation.jl"))
 
 
 include(joinpath("Commons","SolveProblem.jl"))
+include(joinpath("Commons","SolveNlin.jl"))
 
 include("Main.jl")