update projection options

Author: carlodev

examples/2D/TaylorGreen2D.jl

@@ -12,7 +12,7 @@ t0 = 0.0
 dt = 0.01
 tF = 2.5
 vortex_diameter = 1.0
-N = 16
+N = 32
 Re = 1600
 D = 2
 
@@ -22,8 +22,6 @@ rank_partition = (2,2)
 
 
 
-
-
 solver_options = petsc_options(; vel_ksp="gmres", vel_pc="gamg", pres_ksp="cg", pres_pc="gamg")
 
 sprob = StabilizedProblem(VMS(3))
@@ -33,7 +31,7 @@ physicalp = PhysicalParameters(Re=Re, c=vortex_diameter)
 solverp = SolverParameters(matrix_freq_update=1, Number_Skip_Expansion=10e6, M=40,
 petsc_options=solver_options)
 exportp = ExportParameters(printinitial=true, printmodel=true, 
-vtu_export = ["uh","ph","uh_analytic", "ph_analytic"], extra_export=["VelocityError","PressureError"])
+vtu_export = ["uh","ph","uh_analytic", "ph_analytic"], extra_export=["VelocityError","PressureError"], projection_timesteps=[dt, 2*dt, 10*dt, 20*dt])
 
 
 

src/Commons/Equations/StabilizedEquations.jl

@@ -50,10 +50,10 @@ function segregated_equations(u_adv,params::Dict{Symbol,Any},simcase::Simulation
 
     rhs(v) = 0.0
 
-    if typeof(simcase) <: TaylorGreen
-      writevtk(Ω, "Stab_$dt", nsubcells=order, cellfields=["tm"=>Tm, "tc"=>Tc] )
-      @info "Exporting Stabilization Parameters"
-    end
+    # if typeof(simcase) <: TaylorGreen
+    #   writevtk(Ω, "Stab_$dt", nsubcells=order, cellfields=["tm"=>Tm, "tc"=>Tc] )
+    #   @info "Exporting Stabilization Parameters"
+    # end
 
     return Tuu,Tpu,Auu,Aup,Apu,App,ML,S,rhs
 

src/Commons/ParametersDef/Params.jl

@@ -138,6 +138,8 @@ const ALLOWED_VTU_EXPORTS = ["uh", "ph", "uh_analytic", "ph_analytic","uh_avg",
     vtu_export::Vector{String} = ["uh","ph"]
 
     extra_export::Vector{String} = String[]  
+    projection_timesteps::Vector{Float64} = Float64[]
+
     @assert all(x -> x in ALLOWED_EXTRA_EXPORTS, extra_export) "Invalid value in extra_export. Allowed values:$ALLOWED_EXTRA_EXPORTS"
 
 

src/Commons/Projection.jl

@@ -1,14 +1,18 @@
 module Projection
 using Gridap
+using GridapPETSc
+
 using GridapDistributed
 using CSV, DelimitedFiles
 using DataFrames
 using Parameters
 using PartitionedArrays
 using SegregatedVMSSolver
 using SegregatedVMSSolver.ParametersDef
+using SegregatedVMSSolver.SolverOptions
 using SegregatedVMSSolver.Interfaces
 using SegregatedVMSSolver.CreateProblem
+using SegregatedVMSSolver.MatrixCreation
 using SegregatedVMSSolver.ExportUtility: write_to_csv
 
 
@@ -20,8 +24,8 @@ export compute_VMS2_error
 
 function compute_VMS2_error(uh_fine, simcase::TaylorGreen{Periodic},params::Dict{Symbol,Any}, tn::Real)
     @unpack U, dΩ, degree,parts = params
-    @sunpack D = simcase
-    if D == 2
+    @sunpack D, projection_timesteps = simcase
+    if D == 2 && !isempty(    intersect(projection_timesteps, tn))
     ubar, uprime = project_solution(uh_fine, simcase, params, tn)
     norm_cross, norm_re, eps_cross, eps_re = compute_stresses(ubar, uprime, dΩ) 
     write_apriori_analysis(tn, D, norm_cross, norm_re, eps_cross, eps_re, parts)
@@ -37,10 +41,12 @@ function create_coarse_spaces(params,simcase,order::Int64)
         @unpack  model = params
         simcase_coarse = deepcopy(simcase)
         # simcase_coarse.meshp.meshinfo.N .= ones(Int64, D) .* Int(ceil(N[1] / 2))
-        simcase_coarse.sprob.method.order = order -1
+        simcase_coarse.sprob.method.order = 1
         boundary_conditions = create_boundary_conditions(simcase) 
         Vc, Uc, _, _ = creation_fe_spaces(simcase_coarse, model, boundary_conditions)
         merge!(params,Dict(:Uc=>Uc, :Vc=>Vc))
+        mkpath("Projections")
+
     end
 
     @unpack Vc, Uc = params
@@ -63,16 +69,22 @@ function project_solution(uh_fine, simcase::TaylorGreen{Periodic}, params::Dict{
 
     Vc, Uc =   create_coarse_spaces(params,simcase,order)
 
+    
     #L2 projection of uh_fine solution on lower order dimensional space (same mesh)
     # Build weak form
     a(u,v) = ∫( u ⋅ v )dΩ
     l(v)   = ∫( uh_fine ⋅ v )dΩ
 
     # Assemble and solve
-    op = AffineFEOperator(a,l,Uc(tn),Vc(tn))
-    ubar = solve(op)
+    op_proj = AffineFEOperator(a,l,Uc(tn),Vc(tn))
+    solver_proj =    PETScLinearSolver(pres_kspsetup)
+
+    ubar = solve(solver_proj,op_proj)
+    
+    @info "Projection ubar computed"
+    save_path= joinpath("Projections", "TGV_$(D)D_$(tn)_projections")
 
-    println("ubar - solved")
+    writevtk(Ω, save_path, nsubcells=order, cellfields=["uh_fine"=>uh_fine, "uh_projected"=>ubar])
 
     uprime = uh_fine - ubar