gridap / Gridap.jl

Grid-based approximation of partial differential equations in Julia
MIT License
665 stars 94 forks source link

Visualization of the solution of a multifield time dependant problem #1001

Open Mahdimasmoudi0 opened 3 months ago

Mahdimasmoudi0 commented 3 months ago

hi, I am solving an electrophysiology problem using Gridap. I have solved the problem, the issue is that while trying to write my solution in vtu file, I am getting this : ERROR: This function belongs to an interface definition and cannot be used. The issue is within these lines : createpvd("Cardiac_Electrophysiology_transient_solution") do pvd for (uₕ,t) in uₕₜ pvd[t] = createvtk(Ω,"Cardiac_Electrophysiology_transient_solution_$t"*".vtu",cellfields=["uh"=>uₕ[1],"ph"=>uₕ[2]]) end end

This is the complete code : `using Gridap

mesh generation

n = 100 domain = (0,1,0,1) partition = (n,n) model = CartesianDiscreteModel(domain,partition)

Define test and trial spaces

order = 1 reffeᵤ = ReferenceFE(lagrangian,Float64,order) V = TestFESpace(model,reffeᵤ,conformity=:H1) reffeₚ = ReferenceFE(lagrangian,Float64,order) Q = TestFESpace(model,reffeₚ,conformity=:L2) U = TransientTrialFESpace(V) P = TransientTrialFESpace(Q) Y = MultiFieldFESpace([V, Q]) X = TransientMultiFieldFESpace([U, P])

Triangulation and integration quadrature

degree = order Ωₕ = Triangulation(model) dΩ = Measure(Ωₕ,degree)

weak formulation

α=0.01 c1= 8 c2=8 γ=0.002 μ1=0.2 μ2=0.3 bcnst=0.15 D = 0.1 f(ϕ,r)=c1ϕ(ϕ−α)(1-ϕ)-c2rϕ g(ϕ,r)=( γ + μ1r(1/(μ2+ϕ)))(-r-c2ϕ(ϕ-bcnst-1))
res(t,(ϕ,r),(ν,η)) = ∫(∂t(ϕ)ν +∂t(r)η+ D( ∇(ν)⋅(∇(ϕ))) - νf(ϕ,r) - ηg(ϕ,r))dΩ op_AD = TransientFEOperator(res,X,Y)

Transient solver

linear_solver = LUSolver() Δt = 0.05 θ = 0.5 ode_solver = ThetaMethod(linear_solver,Δt,θ) t₀ = 0.0 T = 10.0 Uh0 = interpolate_everywhere(0.0, U(t₀)) # not like tutorial Ph0 = interpolate_everywhere(0.0, P(t₀)) # not like tutorial xh0 = interpolate_everywhere([Uh0, Ph0], X(t₀)) xhs0 = (xh0,) uₕₜ = solve(ode_solver,op_AD,xhs0,t₀,T) createpvd("Cardiac_Electrophysiology_transient_solution") do pvd for (uₕ,t) in uₕₜ pvd[t] = createvtk(Ω,"Cardiac_Electrophysiology_transientsolution$t"*".vtu",cellfields=["uh"=>uₕ[1],"ph"=>uₕ[2]]) end end`