From 613da3dd0fe6cfa334ff1715048f7535d90b97b9 Mon Sep 17 00:00:00 2001 From: principe Date: Thu, 14 Mar 2024 12:23:22 +0100 Subject: [PATCH 1/6] Added old scripts with solid copling, they run with PA<3.0 --- drafts/solid_main.jl | 143 ++++++++++++++++++++++++++ drafts/solid_main_mpi.jl | 216 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 359 insertions(+) create mode 100644 drafts/solid_main.jl create mode 100644 drafts/solid_main_mpi.jl diff --git a/drafts/solid_main.jl b/drafts/solid_main.jl new file mode 100644 index 0000000..07e5ebd --- /dev/null +++ b/drafts/solid_main.jl @@ -0,0 +1,143 @@ +module Solid + +using Gridap +using Gridap.Helpers +using Gridap.Algebra +using Gridap.CellData +using Gridap.ReferenceFEs +using Gridap.Geometry +using GridapMHD: main +using PartitionedArrays + +ν=1.0 +ρ=1.0 +σ=1.0 +σ1=0.1 +σ2=10.0 +B=VectorValue(0.0,50.0,0.0) +f=VectorValue(0.0,0.0,1.0) +u0=1.0 +B0=norm(B) + +# Reduced quantities +L = 1.0 # Do not rescale by length +Re = u0*L/ν +Ha = B0*L*sqrt(σ/(ρ*ν)) +N = Ha^2/Re +f̄ = (L/(ρ*u0^2))*f +B̄ = (1/B0)*B +α = 1.0 +β = 1.0/Re +γ = N +σ̄1 = σ1/σ +σ̄2 = σ2/σ + +Lt = 2 +n = ceil(10*Lt) +const Lf = 1.8 +domain = (-Lt/2,Lt/2,-Lt/2,Lt/2,-0.1,0.1) +cells = (n,n,3) +layer(x,a) = sign(x)*abs(x)^(1/a) +xmap((x,y,z)) = VectorValue(layer(x,3),layer(y,3),z) +model = CartesianDiscreteModel(domain,cells,map=xmap,isperiodic=(false,false,true)) +Ω = Interior(model) + +labels = get_face_labeling(model) +cell_entity = labels.d_to_dface_to_entity[end] +solid_1 = maximum(cell_entity) + 1 +solid_2 = solid_1 + 1 +fluid = solid_2 + 1 + +function find_tag(xs) + tol = 1.0e-9 + if all(x->(x[1]>0.5*Lf-tol)||x[1]<-0.5*Lf+tol,xs) + solid_1 + elseif all(x->(x[2]>0.5*Lf-tol)||x[2]<-0.5*Lf+tol,xs) + solid_2 + else + fluid + end +end +cell_xs = get_cell_coordinates(Ω) +copyto!(cell_entity,map(find_tag,cell_xs)) +add_tag!(labels,"solid_1",[solid_1]) +add_tag!(labels,"solid_2",[solid_2]) +add_tag!(labels,"solid",[solid_1,solid_2]) +add_tag!(labels,"fluid",[fluid]) + +D = num_cell_dims(model) +topo = get_grid_topology(model) +wall = fluid + 1 +for d in 0:D-1 + dface_entity = labels.d_to_dface_to_entity[d+1] + dface_cells = get_faces(topo,d,D) + cache = array_cache(dface_cells) + for dface in 1:length(dface_cells) + cells = getindex!(cache,dface_cells,dface) + solid_found = false + fluid_found = false + for cell in cells + solid_found = solid_found || cell_entity[cell] in (solid_1,solid_2) + fluid_found = fluid_found || cell_entity[cell] == fluid + end + if solid_found && fluid_found + dface_entity[dface] = wall + end + end +end +add_tag!(labels,"wall",[wall]) + +insulating_tags = vcat(collect(1:(8+12)),collect((1:4).+(8+12+2))) +add_tag_from_tags!(labels,"insulating",insulating_tags) + +writevtk(model,"s_model") + +function entity_to_σ(entity) + if entity == solid_1 + σ̄1 + elseif entity == solid_2 + σ̄2 + else + 1.0 + end +end + +σ_Ω = CellField(map(entity_to_σ,cell_entity),Ω) + +params = Dict( + :ptimer => PTimer(get_part_ids(sequential,1),verbose=true), + :debug=>false, + :res_assemble=>false, + :jac_assemble=>false, + :solve=>true, + :model => model, + :solid => Dict(:domain=>"solid",:σ=>σ_Ω), + :fluid => Dict( + :domain=>"fluid", + :α=>α, + :β=>β, + :γ=>γ, + :B=>B̄, + :f=>f̄, + ), + :bcs=>Dict( + :u => Dict(:tags=>"wall"), + :j => Dict(:tags=>"insulating"), + ), +) + +xh, fullparams = main(params) + +ūh,p̄h,j̄h,φ̄h = xh +uh = u0*ūh +ph = (ρ*u0^2)*p̄h +jh = (σ*u0*B0)*j̄h +φh = (u0*B0*L)*φ̄h +dh = ∇⋅jh + +writevtk(model,"s_model") +writevtk(Ω,"s_Ω",order=2,cellfields=["uh"=>uh,"jh"=>jh,"dh"=>dh,"sigma"=>σ_Ω]) + + +end #module + diff --git a/drafts/solid_main_mpi.jl b/drafts/solid_main_mpi.jl new file mode 100644 index 0000000..4aa699e --- /dev/null +++ b/drafts/solid_main_mpi.jl @@ -0,0 +1,216 @@ +module Solid + +using Gridap +using Gridap.Helpers +using Gridap.Algebra +using Gridap.CellData +using Gridap.ReferenceFEs +using Gridap.Geometry +using GridapMHD: main +using GridapDistributed +using PartitionedArrays + +parts = (2,2,1) +ranks = get_part_ids(SequentialBackend(),parts) + +ν=1.0 +ρ=1.0 +σ=1.0 +σ1=0.1 +σ2=10.0 +B=VectorValue(0.0,50.0,0.0) +f=VectorValue(0.0,0.0,1.0) +u0=1.0 +B0=norm(B) + +# Reduced quantities +L = 1.0 # Do not rescale by length +Re = u0*L/ν +Ha = B0*L*sqrt(σ/(ρ*ν)) +N = Ha^2/Re +f̄ = (L/(ρ*u0^2))*f +B̄ = (1/B0)*B +α = 1.0 +β = 1.0/Re +γ = N +σ̄1 = σ1/σ +σ̄2 = σ2/σ + +Lt = 2 +n = ceil(10*Lt) +const Lf = 1.8 +domain = (-Lt/2,Lt/2,-Lt/2,Lt/2,-0.1,0.1) +cells = (n,n,3) +layer(x,a) = sign(x)*abs(x)^(1/a) +xmap((x,y,z)) = VectorValue(layer(x,3),layer(y,3),z) +model = CartesianDiscreteModel(ranks,domain,cells,map=xmap,isperiodic=(false,false,true)) +Ω = Interior(model) + +# This is not really needed, we know that the maximum tag is the +# same in all parts... +solid_tag = map_parts(model.models) do model + labels = get_face_labeling(model) + cell_entity = labels.d_to_dface_to_entity[end] + solid = maximum(cell_entity) + 1 +end +solid_tag = reduce_all(max,solid_tag,init=0) +solid_1 = get_main_part(solid_tag) +solid_2 = solid_1 + 1 +fluid = solid_2 + 1 +wall = fluid + 1 + +function find_tag(xs) + tol = 1.0e-9 + if all(x->(x[1]>0.5*Lf-tol)||x[1]<-0.5*Lf+tol,xs) + solid_1 + elseif all(x->(x[2]>0.5*Lf-tol)||x[2]<-0.5*Lf+tol,xs) + solid_2 + else + fluid + end +end + +D = num_cell_dims(model) + +# Dfaces = get_face_gids(model,D) +# +# cell_entity = map_parts(model.models,Dfaces.partition) do model,partition +# labels = get_face_labeling(model) +# cell_entity = labels.d_to_dface_to_entity[end] +# grid = get_grid(model) +# cell_xs = get_cell_coordinates(grid) +# copyto!(cell_entity,map(find_tag,cell_xs)) +# add_tag!(labels,"solid_1",[solid_1]) +# add_tag!(labels,"solid_2",[solid_2]) +# add_tag!(labels,"solid",[solid_1,solid_2]) +# add_tag!(labels,"fluid",[fluid]) +# # Do we need to make tags consistent? No, we know they already are... +# topo = get_grid_topology(model) +# for d in 0:D-1 +# dface_entity = labels.d_to_dface_to_entity[d+1] +# dface_cells = get_faces(topo,d,D) +# cache = array_cache(dface_cells) +# for dface in 1:length(dface_cells) +# cells = getindex!(cache,dface_cells,dface) +# solid_found = false +# fluid_found = false +# for cell in cells +# solid_found = solid_found || cell_entity[cell] in (solid_1,solid_2) +# fluid_found = fluid_found || cell_entity[cell] == fluid +# end +# if solid_found && fluid_found +# dface_entity[dface] = wall +# end +# end +# end +# add_tag!(labels,"wall",[wall]) +# insulating_tags = vcat(collect(1:(8+12)),collect((1:4).+(8+12+2))) +# add_tag_from_tags!(labels,"insulating",insulating_tags) +# cell_entity[partition.oid_to_lid] +# end + +# Separate previous code in two steps: +# 1) generating tags and +# 2) generating a cell_field to define a variable conductivity based on tags +# Only step 2 needs to be done if tags are generated with, e.g. gmsh + +# 1) Generate tags +# 1.1) Materials: 2 solids and a fluid +# 1.2) Boundaries: non-slip walls (at the interface between a solid and the fluid) +# and insulating walls (at exterior faces) +map_parts(model.models) do model + labels = get_face_labeling(model) + cell_entity = labels.d_to_dface_to_entity[end] + grid = get_grid(model) + cell_xs = get_cell_coordinates(grid) + copyto!(cell_entity,map(find_tag,cell_xs)) + add_tag!(labels,"solid_1",[solid_1]) + add_tag!(labels,"solid_2",[solid_2]) + add_tag!(labels,"solid",[solid_1,solid_2]) + add_tag!(labels,"fluid",[fluid]) + # Do we need to make tags consistent? No, we know they already are... + topo = get_grid_topology(model) + for d in 0:D-1 + dface_entity = labels.d_to_dface_to_entity[d+1] + dface_cells = get_faces(topo,d,D) + cache = array_cache(dface_cells) + for dface in 1:length(dface_cells) + cells = getindex!(cache,dface_cells,dface) + solid_found = false + fluid_found = false + for cell in cells + solid_found = solid_found || cell_entity[cell] in (solid_1,solid_2) + fluid_found = fluid_found || cell_entity[cell] == fluid + end + if solid_found && fluid_found + dface_entity[dface] = wall + end + end + end + add_tag!(labels,"wall",[wall]) + insulating_tags = vcat(collect(1:(8+12)),collect((1:4).+(8+12+2))) + add_tag_from_tags!(labels,"insulating",insulating_tags) + return nothing +end +writevtk(model,"s_model") + +# 2) Get cell entities to define σ +Dfaces = get_face_gids(model,D) +cell_entity = map_parts(model.models,Dfaces.partition) do model,partition + labels = get_face_labeling(model) + cell_entity = labels.d_to_dface_to_entity[end] + cell_entity[partition.oid_to_lid] +end + +function entity_to_σ(entity) + if entity == solid_1 + σ̄1 + elseif entity == solid_2 + σ̄2 + else + 1.0 + end +end + +fields = map_parts(Ω.trians,cell_entity) do trian,cell_entity + CellField(map(entity_to_σ,cell_entity),trian) +end +σ_Ω = GridapDistributed.DistributedCellField(fields) + +# Generate Dict to call GridapMHD +params = Dict( + :ptimer => PTimer(get_part_ids(sequential,1),verbose=true), + :debug=>false, + :res_assemble=>false, + :jac_assemble=>false, + :solve=>true, + :model => model, + :solid => Dict(:domain=>"solid",:σ=>σ_Ω), + :fluid => Dict( + :domain=>"fluid", + :α=>α, + :β=>β, + :γ=>γ, + :B=>B̄, + :f=>f̄, + ), + :bcs=>Dict( + :u => Dict(:tags=>"wall"), + :j => Dict(:tags=>"insulating"), + ), +) + +xh, fullparams = main(params) + +ūh,p̄h,j̄h,φ̄h = xh +uh = u0*ūh +ph = (ρ*u0^2)*p̄h +jh = (σ*u0*B0)*j̄h +φh = (u0*B0*L)*φ̄h +dh = ∇⋅jh + +writevtk(Ω,"s_Ω",order=2,cellfields=["uh"=>uh,"jh"=>jh,"dh"=>dh,"sigma"=>σ_Ω]) + + +end #module + From b151e26e0cd44f653c04b6a2c3d4d9dfb44303d7 Mon Sep 17 00:00:00 2001 From: principe Date: Thu, 14 Mar 2024 12:24:40 +0100 Subject: [PATCH 2/6] Implemented solid coupling for Hunt problem, still in debuging phase --- src/Applications/hunt.jl | 68 +++++++++++++++++++++++++++++++++----- src/Meshers/hunt_mesher.jl | 64 +++++++++++++++++++++++++++++------ src/Meshers/meshers.jl | 2 ++ 3 files changed, 114 insertions(+), 20 deletions(-) diff --git a/src/Applications/hunt.jl b/src/Applications/hunt.jl index eaafafa..7d863fd 100644 --- a/src/Applications/hunt.jl +++ b/src/Applications/hunt.jl @@ -51,6 +51,9 @@ function _hunt(; L=1.0, u0=1.0, B0=norm(VectorValue(B)), + σw1=0.1, + σw2=10.0, + tw=0.0, nsums = 10, vtk=true, title = "test", @@ -105,10 +108,12 @@ function _hunt(; α = 1.0 β = 1.0/Re γ = N + σ̄1 = σw1/σ + σ̄2 = σw2/σ # DiscreteModel in terms of reduced quantities - model = hunt_mesh(parts,params,nc,rank_partition,L,Ha,kmap_x,kmap_y,BL_adapted,ranks_per_level) + model = hunt_mesh(parts,params,nc,rank_partition,L,tw,Ha,kmap_x,kmap_y,BL_adapted,ranks_per_level) Ω = Interior(model) if debug && vtk writevtk(model,"data/hunt_model") @@ -124,6 +129,11 @@ function _hunt(; :ζ=>ζ, ) + if tw > 0.0 + σ_Ω = solid_conductivity(model) + params[:solid] = Dict(:domain=>"solid",:σ=>σ_Ω) + end + # Boundary conditions params[:bcs] = Dict( @@ -188,11 +198,19 @@ function _hunt(; toc!(t,"post_process") if vtk - writevtk(Ω_phys,joinpath(path,title), - order=2, - cellfields=[ - "uh"=>uh,"ph"=>ph,"jh"=>jh,"phi"=>φh, - "u"=>u,"j"=>j,"u_ref"=>u_ref,"j_ref"=>j_ref]) + if tw > 0.0 + writevtk(Ω_phys,joinpath(path,title), + order=2, + cellfields=[ + "uh"=>uh,"ph"=>ph,"jh"=>jh,"phi"=>φh, + "u"=>u,"j"=>j,"u_ref"=>u_ref,"j_ref"=>j_ref,"σ"=>σ_Ω]) + else + writevtk(Ω_phys,joinpath(path,title), + order=2, + cellfields=[ + "uh"=>uh,"ph"=>ph,"jh"=>jh,"phi"=>φh, + "u"=>u,"j"=>j,"u_ref"=>u_ref,"j_ref"=>j_ref]) + end toc!(t,"vtk") end if verbose @@ -227,13 +245,13 @@ end function hunt_mesh( parts,params, - nc::Tuple,np::Tuple,L::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool, + nc::Tuple,np::Tuple,L::Real,tw::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool, ranks_per_level) if isnothing(ranks_per_level) # Single grid - model = Meshers.hunt_generate_base_mesh(parts,np,nc,L,Ha,kmap_x,kmap_y,BL_adapted) + model = Meshers.hunt_generate_base_mesh(parts,np,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted) params[:model] = model else # Multigrid - base_model = Meshers.hunt_generate_base_mesh(nc,L,Ha,kmap_x,kmap_y,BL_adapted) + base_model = Meshers.hunt_generate_base_mesh(nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted) mh = Meshers.generate_mesh_hierarchy(parts,base_model,0,ranks_per_level) params[:multigrid] = Dict{Symbol,Any}( :mh => mh, @@ -246,6 +264,38 @@ function hunt_mesh( return model end +function solid_conductivity(model::GridapDistributed.DistributedDiscreteModel) + D = num_cell_dims(model) + cells = get_face_gids(model,D) + Ω = Interior(model) + fields = map(local_views(model),local_views(cells),local_views(Ω)) do model,partition,trian + labels = get_face_labeling(model) + cell_entity = labels.d_to_dface_to_entity[end] + σ_field(labels,trian,cell_entity[partition.own_to_local]) + end + GridapDistributed.DistributedCellField(fields) +end + +function solid_conductivity(model) + labels = get_face_labeling(model) + σ_field(labels,Interior(model),get_cell_entity(labels)) +end + +function σ_field(labels,Ω,cell_entity) + solid_1 = get_tag_entities(labels,"solid_1") + solid_2 = get_tag_entities(labels,"solid_2") + function entity_to_σ(entity) + if entity == solid_1 + σ̄1 + elseif entity == solid_2 + σ̄2 + else + 1.0 + end + end + σ_field = CellField(map(entity_to_σ,cell_entity),Ω) +end + # This is not very elegant. This needs to be solved by Gridap and GridapDistributed function _warp(model::DiscreteModel,Ω::Triangulation,L) grid_phys = UnstructuredGrid(get_grid(model)) diff --git a/src/Meshers/hunt_mesher.jl b/src/Meshers/hunt_mesher.jl index 7751a63..64562f5 100644 --- a/src/Meshers/hunt_mesher.jl +++ b/src/Meshers/hunt_mesher.jl @@ -45,12 +45,52 @@ function hunt_stretch_map( return coord_map end -function hunt_add_tags!(model) +function hunt_add_tags!(model::GridapDistributed.DistributedDiscreteModel,L::Real,tw::Real) + map(local_views(model)) do model + hunt_add_tags!(model,L,tw) + end +end + +function hunt_add_tags!(model,L::Real,tw::Real) labels = get_face_labeling(model) - tags_u = append!(collect(1:20),[23,24,25,26]) - tags_j = append!(collect(1:20),[25,26]) - add_tag_from_tags!(labels,"noslip",tags_u) - add_tag_from_tags!(labels,"insulating",tags_j) + if tw > 0.0 ## add solid tags + # When model is part of a distributed model we are using + # that the maximum entity (9) is the same in all parts, + # which is true for a GenericDistributedDiscreteModel + # A reduction would be needed in general + # cell_entity = labels.d_to_dface_to_entity[end] + cell_entity = get_cell_entity(labels) + solid_1 = maximum(cell_entity) + 1 + solid_2 = solid_1 + 1 + fluid = solid_2 + 1 + noslip = fluid + 1 + function set_entities(xs) + tol = 1.0e-9 + if all(x->(x[1]>0.5*L-tol)||x[1]<-0.5*L+tol,xs) + solid_1 + elseif all(x->(x[2]>0.5*L-tol)||x[2]<-0.5*L+tol,xs) + solid_2 + else + fluid + end + end + grid = get_grid(model) + cell_coords = get_cell_coordinates(grid) + copyto!(cell_entity,map(set_entities,cell_coords)) + add_tag!(labels,"solid_1",[solid_1]) + add_tag!(labels,"solid_2",[solid_2]) + add_tag!(labels,"solid",[solid_1,solid_2]) + add_tag!(labels,"fluid",[fluid]) + tags_j = vcat(collect(1:(8+12)),collect((1:4).+(8+12+2))) + add_tag_from_tags!(labels,"insulating",tags_j) + add_non_slip_at_solid_entity!(model,[solid_1,solid_2],fluid,noslip) + add_tag!(labels,"noslip",[noslip]) + else + tags_u = append!(collect(1:20),[23,24,25,26]) + tags_j = append!(collect(1:20),[25,26]) + add_tag_from_tags!(labels,"noslip",tags_u) + add_tag_from_tags!(labels,"insulating",tags_j) + end end """ @@ -60,14 +100,15 @@ end """ function hunt_generate_base_mesh( ranks,np::Tuple, - nc::Tuple,L::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool + nc::Tuple,L::Real,tw::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool ) - coord_map = hunt_stretch_map(L,Ha,kmap_x,kmap_y,BL_adapted) + Lt = L+tw + coord_map = hunt_stretch_map(Lt,Ha,kmap_x,kmap_y,BL_adapted) _nc = (nc[1],nc[2],3) _np = (np[1],np[2],1) domain = (-1.0,1.0,-1.0,1.0,0.0,0.1) model = CartesianDiscreteModel(ranks,_np,domain,_nc;isperiodic=(false,false,true),map=coord_map) - hunt_add_tags!(model) + hunt_add_tags!(model,L,tw) return model end @@ -77,12 +118,13 @@ end Generate a serial mesh for the Hunt problem. """ function hunt_generate_base_mesh( - nc::Tuple,L::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool + nc::Tuple,L::Real,tw::Real,Ha::Real,kmap_x::Number,kmap_y::Number,BL_adapted::Bool ) - coord_map = hunt_stretch_map(L,Ha,kmap_x,kmap_y,BL_adapted) + Lt = L+tw + coord_map = hunt_stretch_map(Lt,Ha,kmap_x,kmap_y,BL_adapted) _nc = (nc[1],nc[2],3) domain = (-1.0,1.0,-1.0,1.0,0.0,0.1) model = CartesianDiscreteModel(domain,_nc;isperiodic=(false,false,true),map=coord_map) - hunt_add_tags!(model) + hunt_add_tags!(model,L,tw) return model end diff --git a/src/Meshers/meshers.jl b/src/Meshers/meshers.jl index 209ae23..6224b50 100644 --- a/src/Meshers/meshers.jl +++ b/src/Meshers/meshers.jl @@ -1,10 +1,12 @@ module Meshers using DrWatson using PartitionedArrays, MPI + using GridapDistributed using Gridap, GridapP4est, GridapSolvers using Gridap.Arrays, Gridap.Geometry, Gridap.ReferenceFEs, Gridap.Adaptivity + include("tools.jl") include("p4est.jl") export generate_refined_mesh From b98c9f3057d704c3021e4577fa9a79cb0353d41b Mon Sep 17 00:00:00 2001 From: principe Date: Thu, 14 Mar 2024 14:14:27 +0100 Subject: [PATCH 3/6] Minor fix --- src/Applications/hunt.jl | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/Applications/hunt.jl b/src/Applications/hunt.jl index 7d863fd..613f179 100644 --- a/src/Applications/hunt.jl +++ b/src/Applications/hunt.jl @@ -130,7 +130,7 @@ function _hunt(; ) if tw > 0.0 - σ_Ω = solid_conductivity(model) + σ_Ω = solid_conductivity(Ω,get_cell_gids(model),get_face_labeling(model)) params[:solid] = Dict(:domain=>"solid",:σ=>σ_Ω) end @@ -264,23 +264,15 @@ function hunt_mesh( return model end -function solid_conductivity(model::GridapDistributed.DistributedDiscreteModel) - D = num_cell_dims(model) - cells = get_face_gids(model,D) - Ω = Interior(model) - fields = map(local_views(model),local_views(cells),local_views(Ω)) do model,partition,trian - labels = get_face_labeling(model) +function solid_conductivity(Ω::GridapDistributed.DistributedTriangulation,cells,labels) + D = num_cell_dims(Ω) + fields = map(local_views(labels),local_views(cells),local_views(Ω)) do labels,partition,trian cell_entity = labels.d_to_dface_to_entity[end] σ_field(labels,trian,cell_entity[partition.own_to_local]) end GridapDistributed.DistributedCellField(fields) end -function solid_conductivity(model) - labels = get_face_labeling(model) - σ_field(labels,Interior(model),get_cell_entity(labels)) -end - function σ_field(labels,Ω,cell_entity) solid_1 = get_tag_entities(labels,"solid_1") solid_2 = get_tag_entities(labels,"solid_2") From 68eee08d0d63a57a3864727fcf6968e1c8141b19 Mon Sep 17 00:00:00 2001 From: principe Date: Thu, 14 Mar 2024 14:29:46 +0100 Subject: [PATCH 4/6] Added missing file --- src/Meshers/tools.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 src/Meshers/tools.jl diff --git a/src/Meshers/tools.jl b/src/Meshers/tools.jl new file mode 100644 index 0000000..31fc7ab --- /dev/null +++ b/src/Meshers/tools.jl @@ -0,0 +1,23 @@ +function add_non_slip_at_solid_entity!(model,solid_entities,fluid_entity,name) + D = num_cell_dims(model) + labels = get_face_labeling(model) + topo = get_grid_topology(model) + cell_entity = get_cell_entity(labels) + for d in 0:D-1 + dface_entity = labels.d_to_dface_to_entity[d+1] + dface_cells = get_faces(topo,d,D) + cache = array_cache(dface_cells) + for dface in 1:length(dface_cells) + cells = getindex!(cache,dface_cells,dface) + solid_found = false + fluid_found = false + for cell in cells + solid_found = solid_found || cell_entity[cell] in solid_entities + fluid_found = fluid_found || cell_entity[cell] == fluid_entity + end + if solid_found && fluid_found + dface_entity[dface] = name + end + end + end +end From 1ed47d57c8076897ae96838dbd0f00043282dbb0 Mon Sep 17 00:00:00 2001 From: principe Date: Wed, 27 Mar 2024 13:07:47 +0100 Subject: [PATCH 5/6] Solid coupling working for Hunt problem --- meshes/Expansion_31k.geo | 130 +++++++++++++++++++++++++++++++++++++ src/Applications/hunt.jl | 51 ++++++++++----- src/Meshers/hunt_mesher.jl | 6 +- test/seq/hunt_tests.jl | 16 +++++ 4 files changed, 183 insertions(+), 20 deletions(-) create mode 100644 meshes/Expansion_31k.geo diff --git a/meshes/Expansion_31k.geo b/meshes/Expansion_31k.geo new file mode 100644 index 0000000..a796514 --- /dev/null +++ b/meshes/Expansion_31k.geo @@ -0,0 +1,130 @@ +// Gmsh project created on Wed Jan 19 14:55:13 2022 + +//Puntos de la figura 2D en z=-1.0 para luego extruir + +//+ +Point(1) = {-8, 0, -1.0, 1.0}; +//+ +Point(2) = {-8, 0.25, -1.0, 1.0}; +//+ +Point(3) = {0, 0.25, -1.0, 1.0}; +//+ +Point(4) = {0, 0.625, -1.0, 1.0}; +//+ +Point(5) = {0, 1.0, -1.0, 1.0}; +//+ +Point(6) = {8, 1.0, -1.0, 1.0}; +//+ +Point(7) = {8, 0.625, -1.0, 1.0}; +//+ +Point(8) = {8, 0.25, -1.0, 1.0}; +//+ +Point(9) = {8, 0, -1.0, 1.0}; +//+ +Point(10) = {8, -0.25, -1.0, 1.0}; +//+ +Point(11) = {8, -0.625, -1.0, 1.0}; +//+ +Point(12) = {8, -1.0, -1.0, 1.0}; +//+ +Point(13) = {0, -1.0, -1.0, 1.0}; +//+ +Point(14) = {0, -0.625, -1.0, 1.0}; +//+ +Point(15) = {0, -0.25, -1.0, 1.0}; +//+ +Point(16) = {-8, -0.25, -1.0, 1.0}; + +//Punto extra en el centro +//+ +Point(17) = {0, 0, -1.0, 1.0}; + +//Lineas que unen los puntos + +//+ +Line(1) = {1, 2}; +//+ +Line(2) = {2,3}; +//+ +Line(3) = {3, 4}; +//+ +Line(4) = {4, 5}; +//+ +Line(5) = {5, 6}; +//+ +Line(6) = {6, 7}; +//+ +Line(7) = {7, 8}; +//+ +Line(8) = {8, 9}; +//+ +Line(9) = {9, 10}; +//+ +Line(10) = {10, 11}; +//+ +Line(11) = {11, 12}; +//+ +Line(12) = {12, 13}; +//+ +Line(13) = {13, 14}; +//+ +Line(14) = {14, 15}; +//+ +Line(15) = {15, 16}; +//+ +Line(16) = {16, 1}; +//+ + +//Dos líneas extra para poder definir luego las superficies con sus hexahedros +//+ +Line(17) = {3, 17}; +//+ +Line(18) = {17, 15}; + + +//Defino los nodos de las líneas + +//+ +Transfinite Curve {-2, 15, 5, -12} = 20 Using Progression 1.2; //Nodos axiales +//+ +Transfinite Curve {-1, 17, 8, 16, -18, -9} = 10 Using Progression 1.35; //Nodos dirección Hartmann en el canal pequeño +//+ +Transfinite Curve {3, -7, 13, -11, -14, 10, -4, 6} = 10 Using Progression 1.35; //Nodos dirección Hartmann en el canal pequeño + +//Defino la superfice + +//+ +Curve Loop(1) = {1, 2, 17, 18, 15, 16}; +//+ +Plane Surface(1) = {1}; +//+ +Curve Loop(2) = {-17, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, -18}; +//+ +Plane Surface(2) = {2}; +//+ +Transfinite Surface {1} = {2, 3, 15, 16}; +//+ +Transfinite Surface {2} = {5, 6, 12, 13}; +//+ +Recombine Surface {1,2}; //Esto hace desaparecer los triangulos +//+ +Extrude {0.0,0.0,2.0} { + Surface{1,2}; Layers{{4,10,4}, {0.1,0.9,1}}; Recombine; //Pongo 8 nodos en las capas limites side para Ha = 1000 +} + +//Defino las tags para las condiciones de contorno (los números de las superficies los genera gmsh automáticamente al extruir) + +//+ +Physical Surface("inlet") = {29, 49}; +//+ +Physical Surface("outlet") = {85, 89, 93, 97, 101, 105}; +//+ +Physical Surface("wall") = {1, 2, 33, 45, 109, 81, 50, 122, 77, 73, 117, 113}; +//+ +Physical Volume("PbLi") = {1, 2}; +//+ +Physical Line("inlet") = {27}; +Physical Line("outlet") = {84,88,92,96,100}; +Physical Line("wall") = {1,16,20,25,28,44,2,15,21,24,3,4,17,18,13,14,53,54,22,23,63,64,76,72,32,40,112,108,5,12,55,62,6,7,8,9,10,11,56,57,58,59,60,61,80,104}; +//+ +Physical Point("wall") = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,23,27,31,35,41,45,49,53,57,61,65,69,73,77,81}; diff --git a/src/Applications/hunt.jl b/src/Applications/hunt.jl index 613f179..681dec2 100644 --- a/src/Applications/hunt.jl +++ b/src/Applications/hunt.jl @@ -116,22 +116,31 @@ function _hunt(; model = hunt_mesh(parts,params,nc,rank_partition,L,tw,Ha,kmap_x,kmap_y,BL_adapted,ranks_per_level) Ω = Interior(model) if debug && vtk - writevtk(model,"data/hunt_model") + writevtk(model,"hunt_model") end - params[:fluid] = Dict( - :domain=>nothing, - :α=>α, - :β=>β, - :γ=>γ, - :f=>f̄, - :B=>B̄, - :ζ=>ζ, - ) - if tw > 0.0 - σ_Ω = solid_conductivity(Ω,get_cell_gids(model),get_face_labeling(model)) + σ_Ω = solid_conductivity(σ̄1,σ̄2,Ω,get_cell_gids(model),get_face_labeling(model)) params[:solid] = Dict(:domain=>"solid",:σ=>σ_Ω) + params[:fluid] = Dict( + :domain=>"fluid", + :α=>α, + :β=>β, + :γ=>γ, + :B=>B̄, + :f=>f̄, + :ζ=>ζ, + ) + else + params[:fluid] = Dict( + :domain=>nothing, + :α=>α, + :β=>β, + :γ=>γ, + :f=>f̄, + :B=>B̄, + :ζ=>ζ, + ) end # Boundary conditions @@ -264,18 +273,18 @@ function hunt_mesh( return model end -function solid_conductivity(Ω::GridapDistributed.DistributedTriangulation,cells,labels) +function solid_conductivity(σ̄1,σ̄2,Ω::GridapDistributed.DistributedTriangulation,cells,labels) D = num_cell_dims(Ω) fields = map(local_views(labels),local_views(cells),local_views(Ω)) do labels,partition,trian cell_entity = labels.d_to_dface_to_entity[end] - σ_field(labels,trian,cell_entity[partition.own_to_local]) + σ_field(σ̄1,σ̄2,labels,trian,cell_entity[partition.own_to_local]) end GridapDistributed.DistributedCellField(fields) end -function σ_field(labels,Ω,cell_entity) - solid_1 = get_tag_entities(labels,"solid_1") - solid_2 = get_tag_entities(labels,"solid_2") +function σ_field(σ̄1,σ̄2,labels,Ω,cell_entity) + solid_1 = get_tag_entities(labels,"solid_1")[1] + solid_2 = get_tag_entities(labels,"solid_2")[1] function entity_to_σ(entity) if entity == solid_1 σ̄1 @@ -322,6 +331,10 @@ function analytical_hunt_u( ξ = x[1]/a η = x[2]/a + if ξ > 1.0 || ξ < -1.0 || η > 1.0 || η < -1.0# think about a≠b... + return VectorValue(0.0,0.0,0.0) + end + V = 0.0; V0=0.0; for k in 0:n α_k = (k + 0.5)*π/l @@ -360,6 +373,10 @@ function analytical_hunt_j( ξ = x[1]/a η = x[2]/a + if ξ > 1.0 || ξ < -1.0 || η > 1.0 || η < -1.0# think about a≠b... + return VectorValue(0.0,0.0,0.0) + end + H_dx = 0.0; H_dy = 0.0 for k in 0:n α_k = (k + 0.5)*π/l diff --git a/src/Meshers/hunt_mesher.jl b/src/Meshers/hunt_mesher.jl index 64562f5..0338524 100644 --- a/src/Meshers/hunt_mesher.jl +++ b/src/Meshers/hunt_mesher.jl @@ -38,7 +38,7 @@ function hunt_stretch_map( return z end function map2(x) - layer(x,a) = sign(x)*abs(x)^(1/a) + layer(x,a) = sign(x)*abs(L*x)^(1/a) return VectorValue(layer(x[1],kmap_x),layer(x[2],kmap_y),x[3]) end coord_map = BL_adapted ? map1 : map2 @@ -66,9 +66,9 @@ function hunt_add_tags!(model,L::Real,tw::Real) noslip = fluid + 1 function set_entities(xs) tol = 1.0e-9 - if all(x->(x[1]>0.5*L-tol)||x[1]<-0.5*L+tol,xs) + if all(x->(x[1]>L-tol)||x[1]<-L+tol,xs) solid_1 - elseif all(x->(x[2]>0.5*L-tol)||x[2]<-0.5*L+tol,xs) + elseif all(x->(x[2]>L-tol)||x[2]<-L+tol,xs) solid_2 else fluid diff --git a/test/seq/hunt_tests.jl b/test/seq/hunt_tests.jl index 10a1212..349b45c 100644 --- a/test/seq/hunt_tests.jl +++ b/test/seq/hunt_tests.jl @@ -47,4 +47,20 @@ hunt( solver=:julia, ) +hunt( + nc=(12,12), + backend=:sequential, + np=(2,2), + L=1.0, + B=(0.,50.,0.), + tw=0.2, + debug=false, + vtk=true, + title="hunt-solid", + solver=:julia, + BL_adapted=false, + kmap_x = 3, + kmap_y = 3 +) + end # module From 9ea67293cd2810d923d3b60efbca4e240be4fadb Mon Sep 17 00:00:00 2001 From: Pere Antoni Martorell Date: Fri, 5 Jul 2024 18:32:01 +0200 Subject: [PATCH 6/6] fix function call --- src/Applications/hunt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Applications/hunt.jl b/src/Applications/hunt.jl index 15ab60d..a270750 100644 --- a/src/Applications/hunt.jl +++ b/src/Applications/hunt.jl @@ -277,7 +277,7 @@ function solid_conductivity(σ̄1,σ̄2,Ω::GridapDistributed.DistributedTriangu cell_entity = labels.d_to_dface_to_entity[end] σ_field(σ̄1,σ̄2,labels,trian,cell_entity[partition.own_to_local]) end - GridapDistributed.DistributedCellField(fields) + GridapDistributed.DistributedCellField(fields,Ω) end function σ_field(σ̄1,σ̄2,labels,Ω,cell_entity)