diff --git a/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl b/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl index de09c7d2..bd01282e 100644 --- a/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl +++ b/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl @@ -36,13 +36,12 @@ end function main() - solve = false γ = 2.1e7 # File path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/" - filename = string(path, "Stokes_01cpu_step1429_iter00.gzip.h5") + filename = string(path, "Stokes_01cpu_step01_iter01.gzip.h5") # path ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiverTom/" # filename = string(path, "Stokes_01cpu_step229_iter00.gzip.h5") # Matrix block A @@ -75,13 +74,21 @@ function main() model = ExtractData(filename,"/model/params"); nx, nz, Δx, Δz = model + display(reshape(equ, Int64(nx), Int64(nz+1))[:,end:-1:1]') + display(reshape(eqv, Int64(nx+1), Int64(nz))[:,end:-1:1]') + # Spies Msc = MA - MB*(MD*MC); @show norm(MA-MA') @show norm(MB+MC') @show norm(Msc-Msc') - p = spy(MA-MA') - display(p) + if norm(Msc-Msc')>0 + p = spy(MA-MA') + # p = spy(MA) + @show MA[11,101] + @show MA[101,11] + display(p) + end # Solve diff --git a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl index 7a445de6..4ae67b92 100644 --- a/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl +++ b/JuliaVisualisation/Main_Visualisation_Makie_MD7.jl @@ -2,7 +2,7 @@ import Pkg Pkg.activate(normpath(joinpath(@__DIR__, "."))) using HDF5, Printf, Colors, ColorSchemes, MathTeXEngine, LinearAlgebra, FFMPEG, Statistics using CairoMakie#, GLMakie -const Mak = CairoMakie +Mak = CairoMakie Makie.update_theme!(fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic))) const y = 365*24*3600 @@ -47,7 +47,7 @@ function main() path ="/home/larafriedrichs/repositories/MDOODZ7.0/MDLIB/" #path=raw"C:\Users\49176\OneDrive\Desktop\Test_c_code\\" path="/home/larafriedrichs/repositories/MDOODZ7.0/runs/firstmodel/" - #path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/" + path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/" # path ="/Users/tduretz/Downloads/" # path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/DoubleSubduction_OMP16/" @@ -60,9 +60,9 @@ function main() # path ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/1_NR09/" # File numbers - file_start = 100 + file_start = 1 file_step = 10 - file_end = 100 + file_end = 1 # Select field to visualise # field = :Phases @@ -72,10 +72,10 @@ function main() # field = :PlasticStrainrate # field = :Stress # field = :StrainRate - field = :Pressure + # field = :Pressure # field = :Divergence # field = :Temperature - # field = :Velocity_x + field = :Velocity_x # field = :Velocity_z # field = :Velocity # field = :GrainSize @@ -97,7 +97,7 @@ function main() fabric = false, # add fabric quiver (normal to director) topo = false, σ1_axis = false, - vel_vec = true, + vel_vec = false, ) α_heatmap = 1.0 #0.85 # transparency of heatmap vel_arrow = 5 diff --git a/JuliaVisualisation/Main_Visualisation_Makie_MD7_Residuals.jl b/JuliaVisualisation/Main_Visualisation_Makie_MD7_Residuals.jl new file mode 100644 index 00000000..8be98d6b --- /dev/null +++ b/JuliaVisualisation/Main_Visualisation_Makie_MD7_Residuals.jl @@ -0,0 +1,161 @@ +import Pkg +Pkg.activate(normpath(joinpath(@__DIR__, "."))) +using HDF5, Printf, Colors, ColorSchemes, MathTeXEngine, LinearAlgebra, FFMPEG, Statistics +using CairoMakie, GLMakie +# const Mak = GLMakie +Makie.update_theme!(fonts = (regular = texfont(), bold = texfont(:bold), italic = texfont(:italic))) + +const y = 365*24*3600 +const My = 1e6*y +const cm_y = y*100. + +function ExtractField(filename, field, size, mask_air, mask) + field = try (Float64.(reshape(ExtractData( filename, field), size...))) + catch + @warn "$field not found" + end + mask_air ? field[mask] .= NaN : nothing + return field +end +`d` +function main() + + # Set the path to your files + path ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/" + + # File numbers + istep = 1 + + # Select field to visualise + field = :Residual_x + # field = :Divergence + + # Switches + printfig = false # print figures to disk + printvid = false + framerate = 3 + PlotOnTop = ( + ph_contours = false, # add phase contours + T_contours = true, # add temperature contours + fabric = false, # add fabric quiver (normal to director) + topo = false, + σ1_axis = false, + vel_vec = false, + ) + α_heatmap = 1.0 #0.85 # transparency of heatmap + vel_arrow = 5 + vel_scale = 300000 + nap = 0.1 # pause for animation + resol = 500 + mov_name = "$(path)/_$(field)/$(field)" # Name of the movie + Lx, Lz = 1.0, 1.0 + + # Scaling + # Lc = 1000. + # tc = My + # Vc = 1e-9 + + Lc = 1.0 + tc = My + Vc = 1.0 + + probe = (ϕeff = Float64.([]), t = Float64.([])) + + # Time loop + f = Figure(size = (Lx/Lz*resol*1.2, resol), fontsize=25) + + name = @sprintf("Output%05d.gzip.h5", istep) + filename = string(path, name) + model = ExtractData( filename, "/Model/Params") + xc = ExtractData( filename, "/Model/xc_coord") + zc = ExtractData( filename, "/Model/zc_coord") + xv = ExtractData( filename, "/Model/xg_coord") + zv = ExtractData( filename, "/Model/zg_coord") + xvz = ExtractData( filename, "/Model/xvz_coord") + zvx = ExtractData( filename, "/Model/zvx_coord") + t = model[1] + tMy = round(t/tc, digits=6) + nvx = Int(model[4]) + nvz = Int(model[5]) + ncx, ncz = nvx-1, nvz-1 + centroids, vertices = (ncx, ncz), (nvx, nvz) + + name = @sprintf("Residuals%05d.gzip.h5", istep) + @info "Reading $(name)" + filename = string(path, name) + + Vx = Float64.(reshape(ExtractData( filename, "/VxNodes/ru"), (ncx+1), (ncz+2))) + Vz = Float64.(reshape(ExtractData( filename, "/VzNodes/rv"), (ncx+2), (ncz+1))) + divu = ExtractField(filename, "/Centers/rp", centroids, false, 0) + + # ##################################### + empty!(f) + f = Figure(fontsize=25) #size = (Lx/Lz*resol*1.2, resol), + + if field==:Residual_x + ax1 = Axis(f[1, 1], title = L"$Vx$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + hm = heatmap!(ax1, xv, zvx, Vx[:,:], colormap = (:jet, α_heatmap))#, colorrange=(0., 0.6) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[1, 2], hm, label = L"$Vx$ [cm.yr$^{-1}$]", width = 20, labelsize = 25, ticklabelsize = 14 ) + # Mak.colgap!(f.layout, 20) + if printfig Print2Disk( f, path, string(field), istep) end + end + + if field==:Residual_z + ax2 = Axis(f[1, 1], title = L"$Vz$ at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + hm = heatmap!(ax2, xvz, zv, Vz, colormap = (:jet, α_heatmap))#, colorrange=(0., 0.6) + # colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[1, 2], hm, label = L"$Vz$ [cm.yr$^{-1}$]", width = 20, labelsize = 25, ticklabelsize = 14 ) + # Mak.colgap!(f.layout, 20) + if printfig Print2Disk( f, path, string(field), istep) end + end + + if field==:Divergence + ax3 = Axis(f[1, 1], title = L"∇⋅V at $t$ = %$(tMy) Ma", xlabel = L"$x$ [km]", ylabel = L"$y$ [km]") + hm = heatmap!(ax3, xc./Lc, zc./Lc, divu, colormap = (:turbo, α_heatmap)) + colsize!(f.layout, 1, Aspect(1, Lx/Lz)) + Mak.Colorbar(f[1, 2], hm, label = L"∇⋅V [s$^{-1}$]", width = 20, labelsize = 25, ticklabelsize = 14 ) + Mak.colgap!(f.layout, 20) + if printfig Print2Disk( f, path, string(field), istep) end + end + + DataInspector(f) + display(f) + sleep(nap) + + yscale = Lz/Lx + + if printfig && printvid + FFMPEG.ffmpeg_exe(`-framerate $(framerate) -f image2 -pattern_type glob -i $(path)_$(field)/'*'.png -vf "scale=1080:1080*$(yscale)" -c:v libx264 -pix_fmt yuv420p -y "$(mov_name).mov"`) + end + +end + +function PrincipalStress(τxx, τzz, τxz, P) + σ1 = (x=zeros(size(τxx)), z=zeros(size(τxx)) ) + τxzc = 0.25*(τxz[1:end-1,1:end-1] .+ τxz[2:end-0,1:end-1] .+ τxz[1:end-1,2:end-0] .+ τxz[2:end-0,2:end-0]) + for i in eachindex(τxzc) + if P[i]>1e-13 + σ = [-P[i]+τxx[i] τxzc[i]; τxzc[i] -P[i]+τzz[i]] + v = eigvecs(σ) + σ1.x[i] = v[1,1] + σ1.z[i] = v[2,1] + end + end + return σ1 +end + +function ExtractData( file_path, data_path) + data = h5open(file_path, "r") do file + read(file, data_path) + end + return data +end + +function Print2Disk( f, path, field, istep; res=4) + path1 = path*"/_$field/" + mkpath(path1) + save(path1*"$field"*@sprintf("%05d", istep)*".png", f, px_per_unit = res) +end + +main() \ No newline at end of file diff --git a/MDLIB/Main_DOODZ.c b/MDLIB/Main_DOODZ.c index 93ef219c..c46d71d4 100755 --- a/MDLIB/Main_DOODZ.c +++ b/MDLIB/Main_DOODZ.c @@ -779,6 +779,9 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) { if (input.model.writer_debug == 1 ) WriteResiduals( mesh, input.model, Nmodel, input.scaling ); + // Print2DArrayChar( mesh.BCu.type, input.model.Nx, input.model.Nz+1, 1 ); + // Print2DArrayDouble( mesh.ru, input.model.Nx, input.model.Nz+1, 1.0 ); + // if pass --> clear matrix break if ( (Nmodel.resx < Nmodel.nonlin_abs_mom || Nmodel.resx/Nmodel.resx0 < Nmodel.nonlin_rel_mom) && (Nmodel.resz < Nmodel.nonlin_abs_mom || Nmodel.resz/Nmodel.resz0 < Nmodel.nonlin_rel_mom) && (Nmodel.resp < Nmodel.nonlin_abs_div || Nmodel.resp/Nmodel.resp0 < Nmodel.nonlin_rel_div) ) { printf( "Non-linear solver converged to nonlin_abs_mom = %2.2e nonlin_abs_div = %2.2e nonlin_rel_mom = %2.2e nonlin_rel_div = %2.2e\n", Nmodel.nonlin_abs_mom, Nmodel.nonlin_abs_div, Nmodel.nonlin_rel_mom, Nmodel.nonlin_rel_div ); diff --git a/MDLIB/RheologyDensity.c b/MDLIB/RheologyDensity.c index a0da6e7b..90de8c5b 100644 --- a/MDLIB/RheologyDensity.c +++ b/MDLIB/RheologyDensity.c @@ -1712,8 +1712,9 @@ void StrainRateComponents( grid* mesh, scale scaling, params* model ) { mesh->div_u[c0] = dvxdx + dvzdz + dvydy; // Normal strain rates - mesh->exxd[c0] = dvxdx - 1.0/3.0*mesh->div_u[c0]; - mesh->ezzd[c0] = dvzdz - 1.0/3.0*mesh->div_u[c0]; + mesh->exxd[c0] = dvxdx - model->compressible*1.0/3.0*mesh->div_u[c0]; + mesh->ezzd[c0] = dvzdz - model->compressible*1.0/3.0*mesh->div_u[c0]; + } } diff --git a/MDLIB/Solvers.c b/MDLIB/Solvers.c index 2b06d485..bf6ebd81 100755 --- a/MDLIB/Solvers.c +++ b/MDLIB/Solvers.c @@ -1168,8 +1168,8 @@ void KillerSolver( SparseMat *matA, SparseMat *matB, SparseMat *matC, SparseM L = cs_di_multiply( Dc, Cc ); // //----- test - in case C' != B (assume B is deficient) - // B1 = cs_di_transpose( Cc, 1); - // + // cs_di *B1 = cs_di_transpose( Cc, 1); + // // minus sign: B = -C' // for (k=0; knzmax; k++) B1->x[k] *= -1.0; // could be implicitly included in the next lines diff --git a/MDLIB/StokesAssemblyDecoupled.c b/MDLIB/StokesAssemblyDecoupled.c index 5a648467..98ac8668 100755 --- a/MDLIB/StokesAssemblyDecoupled.c +++ b/MDLIB/StokesAssemblyDecoupled.c @@ -890,6 +890,8 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars double out_of_plane = 1.0; if ( model.out_of_plane==1 ) out_of_plane = 3.0/2.0; double uC_corr = 0.0; + + comp = 0; int iVxC = c1; int iVxW = iVxC-1; @@ -913,15 +915,11 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars iVzNW = c3+nxvz-2; } - // double D11E = mesh->D11_n[iPrE]; - // double D11W = mesh->D11_n[iPrW]; - // double D33N = mesh->D33_s[ixyN]; - // double D33S = mesh->D33_s[ixyS]; - double D11E = 2*mesh->eta_n[iPrE]; - double D11W = 2*mesh->eta_n[iPrW]; - double D33N = mesh->eta_s[ixyN]; - double D33S = mesh->eta_s[ixyS]; - comp = 1.0; + double D11E = mesh->D11_n[iPrE]; + double D11W = mesh->D11_n[iPrW]; + double D33N = mesh->D33_s[ixyN]; + double D33S = mesh->D33_s[ixyS]; + comp = 1; double uS=0.0, uN=0.0, uW=0.0, uE=0.0, uC=0.0, vSW=0.0, vSE=0.0, vNW=0.0, vNE=0.0, pE=0.0, pW=0.0; @@ -930,22 +928,27 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars if (mesh->BCp.type[iPrW] == -1) inW = 1.0; if (mesh->BCp.type[iPrE] == -1) inE = 1.0; - double inS=0.0, inN = 0.0; - if (mesh->BCg.type[ixyS] != 30 && mesh->BCu.type[iVxS] != 13) inS = 1.0; - if (mesh->BCg.type[ixyN] != 30 && mesh->BCu.type[iVxN] != 13) inN = 1.0; + double inS=1.0, inN = 1.0, NeuS = 0.0, NeuN = 0.0, DirS = 0.0, DirN = 0.0; + if (mesh->BCg.type[ixyS] == 30 || mesh->BCu.type[iVxS] == 13) {NeuS = 1.0; inS = 0.0;} + if (mesh->BCg.type[ixyN] == 30 || mesh->BCu.type[iVxN] == 13) {NeuN = 1.0; inN = 0.0;} + if ( mesh->BCu.type[iVxS] == 11) {DirS = 1.0; inS = 0.0;} + if ( mesh->BCu.type[iVxN] == 11) {DirN = 1.0; inN = 0.0;} // Simpler - uW = (1.0/3.0)*D11W*inW*(comp*out_of_plane - 3)/pow(dx, 2); - uC = (1.0/3.0)*(3*pow(dx, 2)*(D33N*inN + D33S*inS) - pow(dz, 2)*(D11E*inE + D11W*inW)*(comp*out_of_plane - 3))/(pow(dx, 2)*pow(dz, 2)); - uE = (1.0/3.0)*D11E*inE*(comp*out_of_plane - 3)/pow(dx, 2); + uW = D11W*inW*(0.33333333333333331*inW*out_of_plane - 1)/pow(dx, 2); + uC = (pow(dx, 2)*(D33N*(2*DirS + inN) + D33S*(2*DirN + inS)) - pow(dz, 2)*(D11E*inE*(0.33333333333333331*inE*out_of_plane - 1) + D11W*inW*(0.33333333333333331*inW*out_of_plane - 1)))/(pow(dx, 2)*pow(dz, 2)); + uE = D11E*inE*(0.33333333333333331*inE*out_of_plane - 1)/pow(dx, 2); uS = -D33S*inS/pow(dz, 2); uN = -D33N*inN/pow(dz, 2); - vSW = ((1.0/3.0)*D11W*comp*inW*out_of_plane - D33S*inS)/(dx*dz); - vSE = (-1.0/3.0*D11E*comp*inE*out_of_plane + D33S*inS)/(dx*dz); - vNW = (-1.0/3.0*D11W*comp*inW*out_of_plane + D33N*inN)/(dx*dz); - vNE = ((1.0/3.0)*D11E*comp*inE*out_of_plane - D33N*inN)/(dx*dz); - pW = -inW*one_dx; - pE = inE*one_dx; + vSW = (0.33333333333333331*D11W*pow(inW, 2)*out_of_plane - D33S)/(dx*dz); + vSE = (-0.33333333333333331*D11E*pow(inE, 2)*out_of_plane + D33S)/(dx*dz); + vNW = (-0.33333333333333331*D11W*pow(inW, 2)*out_of_plane + D33N)/(dx*dz); + vNE = (0.33333333333333331*D11E*pow(inE, 2)*out_of_plane - D33N)/(dx*dz); + pW = -one_dx; + pE = one_dx; + + // if (eqn== 1) printf("uN = %1.4f\n", uN); + // if (eqn==10) printf("uS = %1.4f\n", uS); // Stabilisation with density gradients if ( stab==1 ) { @@ -956,10 +959,6 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars uC += uC_corr; } - // Add contribution from non-conforming Dirichlets - if ( mesh->BCu.type[iVxS] == 11 ) uC -= uS; - if ( mesh->BCu.type[iVxN] == 11 ) uC -= uN; - if ( Assemble == 1 ) { StokesA->b[eqn] *= celvol; StokesB->b[eqn] *= celvol; @@ -968,8 +967,8 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars if (mesh->BCu.type[iVxS] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxS], &(nnzc2A[ith]), uS*celvol, mesh->BCu.type[iVxS], mesh->BCu.val[iVxS], StokesA->bbc); if (mesh->BCu.type[iVxW] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxW], &(nnzc2A[ith]), uW*celvol, mesh->BCu.type[iVxW], mesh->BCu.val[iVxW], StokesA->bbc); AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), uC*celvol, mesh->BCu.type[iVxC], mesh->BCu.val[iVxC], StokesA->bbc); - if (mesh->BCu.type[iVxE] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxE], &(nnzc2A[ith]), uE*celvol, mesh->BCu.type[iVxE], mesh->BCu.val[iVxE], StokesA->bbc); - if (mesh->BCu.type[c1+nx] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxN], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[iVxN], mesh->BCu.val[iVxN], StokesA->bbc); + if (mesh->BCu.type[iVxE] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxE], &(nnzc2A[ith]), uE*celvol, mesh->BCu.type[iVxE], mesh->BCu.val[iVxE], StokesA->bbc); + if (mesh->BCu.type[iVxN] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxN], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[iVxN], mesh->BCu.val[iVxN], StokesA->bbc); //-------------------- if ( mesh->BCv.type[iVzSW] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSW], &(nnzc2A[ith]), vSW*celvol, mesh->BCv.type[iVzSW], mesh->BCv.val[iVzSW], StokesA->bbc); if ( mesh->BCv.type[iVzSE] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSE], &(nnzc2A[ith]), vSE*celvol, mesh->BCv.type[iVzSE], mesh->BCv.val[iVzSE], StokesA->bbc); @@ -996,17 +995,19 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars } } -#if 0 -// MD6 -void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int lev, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int ncx, int nxvz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { +/*--------------------------------------------------------------------------------------------------------------------*/ +/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/ +/*--------------------------------------------------------------------------------------------------------------------*/ + +void Xmomentum_WestNeumannDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int lev, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int ncx, int nxvz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { - double out_of_plane=1.0; - if (model.out_of_plane==1) { - out_of_plane = 3.0/2.0; - } + double dx = mesh->dx; + double dz = mesh->dz; + double out_of_plane = 1.0; + if ( model.out_of_plane==1 ) out_of_plane = 3.0/2.0; + double uC_corr = 0.0; int iVxC = c1; - int iVxW = iVxC-1; int iVxE = iVxC+1; int iVxS = iVxC-nx; int iVxN = iVxC+nx; @@ -1014,65 +1015,46 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars int iVzSE = c3-nxvz+1; int iVzNW = c3; int iVzNE = c3+1; - int iPrW = c2; int iPrE = c2+1; int ixyN = c1; int ixyS = c1-nx; - // Periodic stencil - if ( mesh->BCu.type[iVxC]==-2 ) { - iVxW = c1+nx-2; - iPrW = c2+ncx; - iVzSW = c3-2; - iVzNW = c3+nxvz-2; - } - - double D11E = mesh->D11_n[iPrE]; - double D11W = mesh->D11_n[iPrW]; - double D33N = mesh->D33_s[ixyN]; - double D33S = mesh->D33_s[ixyS]; + double D11E = 2*mesh->eta_n[iPrE]; + double D33N = mesh->eta_s[ixyN]; + double D33S = mesh->eta_s[ixyS]; + comp = 1.0; double uS=0.0, uN=0.0, uW=0.0, uE=0.0, uC=0.0, vSW=0.0, vSE=0.0, vNW=0.0, vNE=0.0, pE=0.0, pW=0.0; - double inE=0.0, inW=0.0, inSu=0.0, inSv=0.0, inNu=0.0, inNv=0.0; - if (mesh->BCp.type[iPrW] == -1) inW = 1.0; + + // Flags + double inE=0.0, inW=0.0; if (mesh->BCp.type[iPrE] == -1) inE = 1.0; - double inS=0.0, inN = 0.0; - if (mesh->BCg.type[ixyS] != 30 && mesh->BCu.type[ixyS ] != 13) inS = 1.0; - if (mesh->BCg.type[c1 ] != 30 && mesh->BCu.type[c1+nx] != 13) inN = 1.0; /// ????????????????????????? - - double dx = mesh->dx; - double dz = mesh->dz; - - comp = 1.0; - - // simpler - uW = (1.0/3.0)*D11W*inW*(comp*out_of_plane - 3)/pow(dx, 2); - uC = (1.0/3.0)*(3*pow(dx, 2)*(D33N*pow(inN, 2) + D33S*pow(inS, 2)) - pow(dz, 2)*(D11E*inE + D11W*inW)*(comp*out_of_plane - 3))/(pow(dx, 2)*pow(dz, 2)); - uE = (1.0/3.0)*D11E*inE*(comp*out_of_plane - 3)/pow(dx, 2); - uS = -D33S*pow(inS, 2)/pow(dz, 2); - uN = -D33N*pow(inN, 2)/pow(dz, 2); - vSW = ((1.0/3.0)*D11W*comp*inW*out_of_plane - D33S*inS)/(dx*dz); - vSE = (-1.0/3.0*D11E*comp*inE*out_of_plane + D33S*inS)/(dx*dz); - vNW = (-1.0/3.0*D11W*comp*inW*out_of_plane + D33N*inN)/(dx*dz); - vNE = ((1.0/3.0)*D11E*comp*inE*out_of_plane - D33N*inN)/(dx*dz); - + double inS=1.0, inN = 1.0, NeuS = 0.0, NeuN = 0.0, DirS = 0.0, DirN = 0.0; + if (mesh->BCg.type[ixyS] == 30 || mesh->BCu.type[iVxS] == 13) {NeuS = 1.0; inS = 0.0;} + if (mesh->BCg.type[ixyN] == 30 || mesh->BCu.type[iVxN] == 13) {NeuN = 1.0; inN = 0.0;} + if ( mesh->BCu.type[iVxS] == 11) {DirS = 1.0; inS = 0.0;} + if ( mesh->BCu.type[iVxN] == 11) {DirN = 1.0; inN = 0.0;} + + + // Simpler + uC = (-D11E*pow(dz, 2)*inE*(0.33333333333333331*inE*out_of_plane - 1) + pow(dx, 2)*(D33N*(2*DirS + inN) + D33S*(2*DirN + inS)))/(pow(dx, 2)*pow(dz, 2)); + uE = D11E*inE*(0.33333333333333331*inE*out_of_plane - 1)/pow(dx, 2); + uS = -D33S*inS/pow(dz, 2); + uN = -D33N*inN/pow(dz, 2); + vSW = -D33S/(dx*dz); + vSE = (-0.33333333333333331*D11E*pow(inE, 2)*out_of_plane + D33S)/(dx*dz); + vNW = D33N/(dx*dz); + vNE = (0.33333333333333331*D11E*pow(inE, 2)*out_of_plane - D33N)/(dx*dz); + pE = inE*one_dx; + // Stabilisation with density gradients - if (stab==1) { + if ( stab==1 ) { double drhodx = (mesh->rho_n[c2+1] - mesh->rho_n[c2])*one_dx; - double uC_corr = 1.00 * om * model.dt * mesh->gx[c1] * drhodx; + uC_corr = 1.00 * om * model.dt * mesh->gx[c1] * drhodx; // Importante trique, voire meme gigantesque! - if (uC+uC_corr>0.0) uC += uC_corr; - } - - // Add contribution from non-conforming Dirichlets - if ( mesh->BCu.type[iVxS] == 11 ) uC -= uS; - if ( mesh->BCu.type[iVxN] == 11 ) uC -= uN; - - // Pressure gradient - if ( mesh->BCp.type[iPrE] != 30 && mesh->BCp.type[iPrW] != 30 ) { - pW = -one_dx * 1e0; - pE = one_dx * 1e0; + if (uC+uC_corr<0.0) uC_corr = 0.0; + uC += uC_corr; } if ( Assemble == 1 ) { @@ -1080,248 +1062,38 @@ void Xmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars StokesB->b[eqn] *= celvol; //-------------------- // dsxx/dx - normal stencil - if (mesh->BCu.type[iVxS] != 30) { - if (mesh->BCu.type[iVxS] != 11) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxS], &(nnzc2A[ith]), uS*celvol, mesh->BCu.type[iVxS], mesh->BCu.val[iVxS], StokesA->bbc); - else AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxS], &(nnzc2A[ith]), uS*celvol, mesh->BCu.type[iVxS], 2.0*mesh->BCu.val[iVxS], StokesA->bbc); - } - if (mesh->BCu.type[iVxW] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxW], &(nnzc2A[ith]), uW*celvol, mesh->BCu.type[iVxW], mesh->BCu.val[iVxW], StokesA->bbc); + if (mesh->BCu.type[iVxS] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxS], &(nnzc2A[ith]), uS*celvol, mesh->BCu.type[iVxS], mesh->BCu.val[iVxS], StokesA->bbc); AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), uC*celvol, mesh->BCu.type[iVxC], mesh->BCu.val[iVxC], StokesA->bbc); - if (mesh->BCu.type[iVxE] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxE], &(nnzc2A[ith]), uE*celvol, mesh->BCu.type[iVxE], mesh->BCu.val[iVxE], StokesA->bbc); - if (mesh->BCu.type[c1+nx] != 30) { - if (mesh->BCu.type[iVxN] != 11) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxN], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[iVxN], mesh->BCu.val[iVxN], StokesA->bbc); - else AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxN], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[iVxN], 2*mesh->BCu.val[iVxN], StokesA->bbc); - } + if (mesh->BCu.type[iVxE] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxE], &(nnzc2A[ith]), uE*celvol, mesh->BCu.type[iVxE], mesh->BCu.val[iVxE], StokesA->bbc); + if (mesh->BCu.type[iVxN] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[iVxN], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[iVxN], mesh->BCu.val[iVxN], StokesA->bbc); //-------------------- if ( mesh->BCv.type[iVzSW] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSW], &(nnzc2A[ith]), vSW*celvol, mesh->BCv.type[iVzSW], mesh->BCv.val[iVzSW], StokesA->bbc); if ( mesh->BCv.type[iVzSE] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzSE], &(nnzc2A[ith]), vSE*celvol, mesh->BCv.type[iVzSE], mesh->BCv.val[iVzSE], StokesA->bbc); if ( mesh->BCv.type[iVzNW] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzNW], &(nnzc2A[ith]), vNW*celvol, mesh->BCv.type[iVzNW], mesh->BCv.val[iVzNW], StokesA->bbc); if ( mesh->BCv.type[iVzNE] != 30 ) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[iVzNE], &(nnzc2A[ith]), vNE*celvol, mesh->BCv.type[iVzNE], mesh->BCv.val[iVzNE], StokesA->bbc); - - //-------------------- - if ( mesh->BCp.type[iPrE] != 30 && mesh->BCp.type[iPrW] != 30 ) { - AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrW] - Stokes->neq_mom, &(nnzc2B[ith]), pW*celvol, mesh->BCp.type[iPrW], mesh->BCp.val[iPrW], StokesB->bbc); - AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrE] - Stokes->neq_mom, &(nnzc2B[ith]), pE*celvol, mesh->BCp.type[iPrE], mesh->BCp.val[iPrE], StokesB->bbc); - } - - } - else { + // printf("%d\n", Stokes->eqn_v[iVzSW]); + // printf("%d %d\n", eqn, Stokes->eqn_v[iVzSE]); + // printf("%d\n", Stokes->eqn_v[iVzNW]); + // printf("%d %d\n", eqn, Stokes->eqn_v[iVzNE]); - // Residual function - // StokesA->F[eqn] = 0.0; - // StokesA->F[eqn] += (mesh->sxxd[iPrE] - mesh->sxxd[iPrW])/dx; - // StokesA->F[eqn] -= (mesh->p_corr[iPrE] - mesh->p_corr[iPrW])/dx; - // StokesA->F[eqn] += (mesh->sxz[ixyN] - mesh->sxz[ixyS]) /dz; - // StokesA->F[eqn] *= -1.0; - // StokesA->F[eqn] -= (StokesA->b[eqn]); - // StokesA->F[eqn] *= celvol; - - // Residual function - StokesA->F[eqn] = uC*u[iVxC]; - if ( mesh->BCp.type[iPrE] != 30 && mesh->BCp.type[iPrW] != 30 ) { - StokesA->F[eqn] += pW*p[iPrW] + pE*p[iPrE]; - } - if ( mesh->BCv.type[iVzSE] != 30 ) StokesA->F[eqn] += vSE*v[iVzSE]; - if ( mesh->BCv.type[iVzSW] != 30 ) StokesA->F[eqn] += vSW*v[iVzSW]; - if ( mesh->BCv.type[iVzNE] != 30 ) StokesA->F[eqn] += vNE*v[iVzNE]; - if ( mesh->BCv.type[iVzNW] != 30 ) StokesA->F[eqn] += vNW*v[iVzNW]; - if ( mesh->BCu.type[iVxS] != 30 && mesh->BCu.type[iVxS] != 11 ) StokesA->F[eqn] += uS*u[iVxS]; - // if ( mesh->BCu.type[iVxS] == 11 ) StokesA->F[eqn] += -2.0*D33S*one_dz_dz*mesh->BCu.val[iVxS]; - if ( mesh->BCu.type[iVxS] == 11 ) StokesA->F[eqn] += 2.0*uS*mesh->BCu.val[iVxS]; - if ( mesh->BCu.type[iVxN] != 30 && mesh->BCu.type[iVxN] != 11 ) StokesA->F[eqn] += uN*u[iVxN]; - // if ( mesh->BCu.type[iVxN] == 11 ) StokesA->F[eqn] += -2.0*D33N*one_dz_dz*mesh->BCu.val[iVxN]; - if ( mesh->BCu.type[iVxN] == 11 ) StokesA->F[eqn] += 2.0*uN*mesh->BCu.val[iVxN]; - if ( mesh->BCu.type[iVxW] != 30 ) StokesA->F[eqn] += uW*u[iVxW]; - if ( mesh->BCu.type[iVxE] != 30 ) StokesA->F[eqn] += uE*u[iVxE]; - - // printf("before %2.2e %2.2e\n", StokesA->F[eqn], StokesA->b[eqn]); - StokesA->F[eqn] -= (StokesA->b[eqn]);// + Stokes->bbc[eqn]; - // printf("%2.2e\n", StokesA->b[eqn]); - // printf("after %2.2e %2.2e\n", StokesA->F[eqn], StokesA->b[eqn]); - StokesA->F[eqn] *= celvol; - if (isnan(StokesA->F[eqn]) ) printf("cuC = %2.2e u[iVxC] = %2.2e, D11W = %2.2e, D11E = %2.2e D33S = %2.2e D33N = %2.2e \n", uC, u[iVxC], D11W, D11E, D33S, D33N ); - } -} -#endif -/*--------------------------------------------------------------------------------------------------------------------*/ -/*------------------------------------------------------ M-Doodz -----------------------------------------------------*/ -/*--------------------------------------------------------------------------------------------------------------------*/ - -void Xmomentum_WestNeumannDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int lev, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int ncx, int nxvz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { - - double etaE, etaW, etaN, etaS; - etaE = mesh->eta_n[c2+1]; - etaW = etaE; - etaN = mesh->eta_s[c1]; - etaS = mesh->eta_s[c1-nx]; - - // double rhoVx = 0.5*(mesh->rho_s[c1] + mesh->rho_s[c1-nx]); - - double uS=0.0, uN=0.0, uW=0.0, uE=0.0, uC=0.0, vSW=0.0, vSE=0.0, vNW=0.0, vNE=0.0, pE=0.0, pW=0.0; - - StokesA->b[eqn] -= mesh->BCu.val[c1]*one_dx; - - // dsxx/dx - if ( mesh->BCu.type[c1+1] != 30 ) { - uC = -2.0*one_dx_dx * (etaE) + comp*2.0/3.0*one_dx_dx * etaE; - uE = 2.0*one_dx_dx * etaE - comp*2.0/3.0*one_dx_dx * etaE; - } - - // eta*dvx/dz - if ( mesh->BCu.type[c1-nx] != 30 && mesh->BCu.type[c1+nx] != 30 ) { - if ( mesh->BCu.type[c1-nx] != 13 ) {uS = one_dz_dz * etaS; uC += -one_dz_dz * etaS;} - if ( mesh->BCu.type[c1+nx] != 13 ) {uN = one_dz_dz * etaN; uC += -one_dz_dz * etaN;} - } - if ( mesh->BCu.type[c1-nx] == 30 && mesh->BCu.type[c1+nx] != 30 ) { - if ( mesh->BCu.type[c1+nx] != 13 ) uC += -one_dz_dz * etaN; - if ( mesh->BCu.type[c1+nx] != 13 ) uN = one_dz_dz * etaN; - } - if ( mesh->BCu.type[c1-nx] != 30 && mesh->BCu.type[c1+nx] == 30 ) { - if ( mesh->BCu.type[c1-nx] != 13 ) uC += -one_dz_dz * etaS; - if ( mesh->BCu.type[c1-nx] != 13 ) uS = one_dz_dz * etaS; - } - - if ( mesh->BCu.type[c1-nx] == 11 ) uC += -one_dz_dz * etaS; - if ( mesh->BCu.type[c1+nx] == 11 ) uC += -one_dz_dz * etaN; - - // eta*dvz/dx - if ( mesh->BCv.type[c3-nxvz] != 30 && mesh->BCv.type[c3-nxvz+1] != 30 ) { - vSE = -one_dx_dz * etaS + comp*2.0/3.0*one_dx_dz * etaE; - vSW = one_dx_dz * etaS - comp*2.0/3.0*one_dx_dz * etaW; - } - - if ( mesh->BCv.type[c3+1] != 30 && mesh->BCv.type[c3] != 30 ) { - vNE = one_dx_dz * etaN - comp*2.0/3.0*one_dx_dz * etaE; - vNW = -one_dx_dz * etaN + comp*2.0/3.0*one_dx_dz * etaW; - } - - // Pressure gradient - if ( mesh->BCp.type[c2+1] != 30 ) { - pW = one_dx; - pE = -pW; - } - - // // Inertia - // if ( model.isinertial == 1 || model.isinertial == 2 ) { - // uC -= sign * rhoVx/model.dt; - // } - // if ( model.isinertial == 2 ) { - // uN -= rhoVx/2*one_dz*mesh->VzVx[c1]; - // uS += rhoVx/2*one_dz*mesh->VzVx[c1]; - // uW += rhoVx/2*one_dx*mesh->u_in[c1]; - // uE -= rhoVx/2*one_dx*mesh->u_in[c1]; - // } - - uS=-uS, uN=-uN, uW=-uW, uE=-uE, uC=-uC, vSW=-vSW, vSE=-vSE, vNW=-vNW, vNE=-vNE, pE=-pE, pW=-pW; - - if ( Assemble == 1 ) { - StokesA->b[eqn] *= celvol; - StokesB->b[eqn] *= celvol; //-------------------- - // dsxx/dx - normal stencil - if (mesh->BCu.type[c1-nx] != 30) { - if (mesh->BCu.type[c1-nx] != 11) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1-nx], &(nnzc2A[ith]), uS*celvol, mesh->BCu.type[c1-nx], mesh->BCu.val[c1-nx], StokesA->bbc); - else AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1-nx], &(nnzc2A[ith]), uS*celvol, mesh->BCu.type[c1-nx], 2.0*mesh->BCu.val[c1-nx], StokesA->bbc); - } - AddCoeff3( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), uC*celvol, mesh->BCu.type[c1], mesh->BCu.val[c1], StokesA->bbc); - if (mesh->BCu.type[c1+1] != 30) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+1], &(nnzc2A[ith]), uE*celvol, mesh->BCu.type[c1+1], mesh->BCu.val[c1+1], StokesA->bbc); - if (mesh->BCu.type[c1+nx] != 30) { - if (mesh->BCu.type[c1+nx] != 11) AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[c1+nx], mesh->BCu.val[c1+nx], StokesA->bbc); - else AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx], &(nnzc2A[ith]), uN*celvol, mesh->BCu.type[c1+nx], 2*mesh->BCu.val[c1+nx], StokesA->bbc); - } - //-------------------- - - if ( mesh->BCv.type[c3-nxvz] != 30 && mesh->BCv.type[c3-nxvz+1] != 30 ) { - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-nxvz], &(nnzc2A[ith]), vSW*celvol, mesh->BCv.type[c3-nxvz], mesh->BCv.val[c3-nxvz], StokesA->bbc); - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-nxvz+1], &(nnzc2A[ith]), vSE*celvol, mesh->BCv.type[c3-nxvz+1], mesh->BCv.val[c3-nxvz+1], StokesA->bbc); - } - - if ( mesh->BCv.type[c3] != 30 && mesh->BCv.type[c3+1] != 30 ) { - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3], &(nnzc2A[ith]), vNW*celvol, mesh->BCv.type[c3], mesh->BCv.val[c3], StokesA->bbc); - AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+1], &(nnzc2A[ith]), vNE*celvol, mesh->BCv.type[c3+1], mesh->BCv.val[c3+1], StokesA->bbc); - } - //-------------------- - if ( mesh->BCp.type[c2+1] != 30 ) { - AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[c2+1] - Stokes->neq_mom, &(nnzc2B[ith]), pE*celvol, mesh->BCp.type[c2+1], mesh->BCp.val[c2+1], StokesB->bbc); - } + if ( mesh->BCp.type[iPrE] != 30) { + AddCoeff3( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[iPrE] - Stokes->neq_mom, &(nnzc2B[ith]), pE*celvol, mesh->BCp.type[iPrE], mesh->BCp.val[iPrE], StokesB->bbc); + } } else { - - // RHS - StokesA->F[eqn] = -(StokesA->b[eqn] + StokesA->bbc[eqn]/celvol); - StokesA->F[eqn] += uC*u[c1]; - if (mesh->BCu.type[c1-nx] != 30) { - if (mesh->BCu.type[c1-nx] != 11) StokesA->F[eqn] += uS*u[c1-nx]; - } - if (mesh->BCu.type[c1+1] != 30) StokesA->F[eqn] += uE*u[c1+1]; - if (mesh->BCu.type[c1+nx] != 30) { - if (mesh->BCu.type[c1+nx] != 11) StokesA->F[eqn] += uN*u[c1+nx]; - } - //-------------------- - if ( mesh->BCv.type[c3-nxvz] != 30 && mesh->BCv.type[c3-nxvz+1] != 30 ) { - if (mesh->BCv.type[c3-nxvz+1] != 11 && mesh->BCv.type[c3-nxvz+1] != 0) StokesA->F[eqn] += vSE*v[c3-nxvz+1]; - if (mesh->BCv.type[c3-nxvz ] != 11 && mesh->BCv.type[c3-nxvz ] != 0) StokesA->F[eqn] += vSW*v[c3-nxvz]; - } - if ( mesh->BCv.type[c3] != 30 && mesh->BCv.type[c3+1] != 30 ) { - if (mesh->BCv.type[c3+1] != 11 && mesh->BCv.type[c3+1] != 0) StokesA->F[eqn] += vNE*v[c3+1]; - if (mesh->BCv.type[c3 ] != 11 && mesh->BCv.type[c3 ] != 0) StokesA->F[eqn] += vNW*v[c3]; - } - //-------------------- - if ( mesh->BCp.type[c2+1] != 30 ) StokesA->F[eqn] += pE*p[c2+1]; + // Residual function + const double SxxE = mesh->sxxd[iPrE] - mesh->p_corr[iPrE]; + const double SxxW = 2*mesh->BCu.val[iVxC] - SxxE; + StokesA->F[eqn] = 0.0; + StokesA->F[eqn] += (SxxE - SxxW )/dx; + StokesA->F[eqn] += (mesh->sxz[ixyN] - mesh->sxz[ixyS])/dz*2; // dodgy factor 2 make it symmetric + StokesA->F[eqn] *= -1/2; + StokesA->F[eqn] -= StokesA->b[eqn] - uC_corr*u[iVxC]; // no body force StokesA->F[eqn] *= celvol; - - // printf("%2.4f %2.4f\n", u[c1-nx], uS); - // printf("%2.4f %2.4f\n", u[c1], uC); - // printf("%2.4f %2.4f\n", u[c1+1],uE); - //// printf("%2.2e\n", u[c1+nx]); - //// printf("%2.2e\n", v[c3-nxvz]); - // printf("%2.4f %2.4f\n", v[c3-nxvz+1], vSE); - //// printf("%2.2e\n", v[c3]); - //// printf("%2.2e\n", v[c3+1]); - // printf("%2.4f %2.4f\n", p[c2+1], pE); - // printf("F1 = %2.4f\n", u[c1-nx]*uS+u[c1]*uC+u[c1+1]*uE+p[c2+1]*pE+v[c3-nxvz+1]*vSE); - - - // // Residual function - // StokesA->F[eqn] = uC*u[c1]; - // if ( mesh->BCp.type[c2+1] != 30 ) { - // StokesA->F[eqn] += pE*p[c2+1]; - // } - // if ( mesh->BCv.type[c3-nxvz+1] == -1 || mesh->BCv.type[c3-nxvz+1] == 2 ) StokesA->F[eqn] += vSE*v[c3-nxvz+1]; - // if ( mesh->BCv.type[c3-nxvz] == -1 || mesh->BCv.type[c3-nxvz] == 2 ) StokesA->F[eqn] += vSW*v[c3-nxvz]; - // if ( mesh->BCv.type[c3+1] == -1 || mesh->BCv.type[c3+1] == 2 ) StokesA->F[eqn] += vNE*v[c3+1]; - // if ( mesh->BCv.type[c3] == -1 || mesh->BCv.type[c3] == 2 ) StokesA->F[eqn] += vNW*v[c3]; - // - // if ( mesh->BCu.type[c1-nx] != 30 && mesh->BCu.type[c1-nx] != 11 ) StokesA->F[eqn] += uS*u[c1-nx]; - // if ( mesh->BCu.type[c1+nx] != 30 && mesh->BCu.type[c1+nx] != 11 ) StokesA->F[eqn] += uN*u[c1+nx]; - // if ( mesh->BCu.type[c1+1] != 30 ) StokesA->F[eqn] += uE*u[c1+1]; - // if ( mesh->BCu.type[c1-nx] == 11 ) StokesA->F[eqn] += -2*etaS*one_dz_dz*mesh->BCu.val[c1-nx]; - // if ( mesh->BCu.type[c1+nx] == 11 ) StokesA->F[eqn] += -2*etaN*one_dz_dz*mesh->BCu.val[c1+nx]; - // - // StokesA->F[eqn] -= (StokesA->b[eqn]) + StokesA->bbc[eqn]; - // StokesA->F[eqn] *= celvol; - - - // printf("\n F = %2.4e %2.2e %2.2e %2.2e\n\n", StokesA->F[eqn], StokesA->b[eqn], StokesA->bbc[eqn], StokesA->b[eqn] + StokesA->bbc[eqn]); - // printf("%2.4e %2.4e %2.4e %2.4e %2.4e\n", u[c1-nx], u[c1], u[c1+1], v[c3-nxvz+1], p[c2+1] ); - // printf("F = %2.2e\n", (StokesA->b[eqn]+ StokesA->bbc[eqn] ) - (u[c1-nx]*uS + u[c1]*uC + u[c1+1]*uE + v[c3-nxvz+1]*vSE + p[c2+1]*pE)); - } - - //// if (mesh->BCu.type[c1+nx]==30) { - // printf("x momentum W -- etaE = %2.2e etaS = %2.2e etaN = %2.2e\n", etaE, etaS, etaN); - // printf("uC = %02d %2.2e\n", mesh->BCu.type[c1+0], uC ); - // // printf("uW = %02d %2.2e\n", mesh->BCu.type[c1-1], uW ); - // printf("uE = %02d %2.2e\n", mesh->BCu.type[c1+1], uE ); - // printf("uS = %02d %2.2e \n", mesh->BCu.type[c1-nx], uS ); - // printf("uN = %02d %2.2e \n", mesh->BCu.type[c1+nx], uN ); - // printf("vSE= %02d %2.2e \n", mesh->BCv.type[c3-nxvz+1], vSE ); - // printf("vNE= %02d %2.2e \n", mesh->BCv.type[c3+1], vNE ); - // printf("vSW= %02d %2.2e \n", mesh->BCv.type[c3-nxvz], vSW ); - // printf("vNW= %02d %2.2e \n", mesh->BCv.type[c3], vNW ); - //// // printf("pW = %02d %2.2e \n", mesh->BCp.type[c2], pW ); - // printf("pE = %02d %2.2e \n",mesh->BCp.type[c2+1], pE); - //// } + } } /*--------------------------------------------------------------------------------------------------------------------*/ @@ -1493,11 +1265,7 @@ void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars double D22N = mesh->D22_n[iPrN]; double D33E = mesh->D33_s[ixyE]; double D33W = mesh->D33_s[ixyW]; - // double D22S = 2*mesh->eta_n[iPrS]; - // double D22N = 2*mesh->eta_n[iPrN]; - // double D33E = mesh->eta_s[ixyE]; - // double D33W = mesh->eta_s[ixyW]; - comp = 1.0; + comp = 1; double vS=0.0, vN=0.0, vW=0.0, vE=0.0, vC=0.0, uSW=0.0, uSE=0.0, uNW=0.0, uNE=0.0, pN=0.0, pS=0.0; @@ -1506,23 +1274,27 @@ void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars if (mesh->BCp.type[iPrS] == -1) inS = 1.0; if (mesh->BCp.type[iPrN] == -1) inN = 1.0; - double inW=0.0, inE = 0.0; - if (mesh->BCg.type[ixyW] != 30 && mesh->BCv.type[iVzW] != 13) inW = 1.0; - if (mesh->BCg.type[ixyE] != 30 && mesh->BCv.type[iVzE] != 13) inE = 1.0; + double inW=1.0, inE = 1.0, NeuW = 0., NeuE = 0, DirW = 0., DirE = 0.; + if (mesh->BCg.type[ixyW] == 30 || mesh->BCv.type[iVzW] == 13) {NeuW = 1.0; inW = 0.0;} + if (mesh->BCg.type[ixyE] == 30 || mesh->BCv.type[iVzE] == 13) {NeuE = 1.0; inE = 0.0;} + if ( mesh->BCv.type[iVzW] == 11) {DirW = 1.0; inW = 0.0;} + if ( mesh->BCv.type[iVzE] == 11) {DirE = 1.0; inE = 0.0;} // Simpler vW = -D33W*inW/pow(dx, 2); - vC = (1.0/3.0)*(-pow(dx, 2)*(D22N*inN + D22S*inS)*(comp*out_of_plane - 3) + 3*pow(dz, 2)*(D33E*inE + D33W*inW))/(pow(dx, 2)*pow(dz, 2)); + vC = (-pow(dx, 2)*(D22N + D22S)*(0.33333333333333331*out_of_plane - 1) + pow(dz, 2)*(D33E*(2*DirE + inE) + D33W*(2*DirW + inW)))/(pow(dx, 2)*pow(dz, 2)); vE = -D33E*inE/pow(dx, 2); - vS = (1.0/3.0)*D22S*inS*(comp*out_of_plane - 3)/pow(dz, 2); - vN = (1.0/3.0)*D22N*inN*(comp*out_of_plane - 3)/pow(dz, 2); - uSW = ((1.0/3.0)*D22S*comp*inS*out_of_plane - D33W*inW)/(dx*dz); - uSE = (-1.0/3.0*D22S*comp*inS*out_of_plane + D33E*inE)/(dx*dz); - uNW = (-1.0/3.0*D22N*comp*inN*out_of_plane + D33W*inW)/(dx*dz); - uNE = ((1.0/3.0)*D22N*comp*inN*out_of_plane - D33E*inE)/(dx*dz); - pS = -inS*one_dz; - pN = inN*one_dz; - + vS = D22S*(0.33333333333333331*out_of_plane - 1)/pow(dz, 2); + vN = D22N*(0.33333333333333331*out_of_plane - 1)/pow(dz, 2); + uSW = (0.33333333333333331*D22S*out_of_plane - D33W)/(dx*dz); + uSE = (-0.33333333333333331*D22S*out_of_plane + D33E)/(dx*dz); + uNW = (-0.33333333333333331*D22N*out_of_plane + D33W)/(dx*dz); + uNE = (0.33333333333333331*D22N*out_of_plane - D33E)/(dx*dz); + pS = -one_dz; + pN = one_dz; + + // if (eqn==100) printf("uSW = %2.2lf\n", uSW); + // Stabilisation with density gradients if ( stab==1 ) { double drhodz = (mesh->rho_n[c2+ncx] - mesh->rho_n[c2])*one_dz; @@ -1531,10 +1303,6 @@ void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars if ((vC+vC_corr)<0.0 || (vC+vC_corr)BCv.type[iVzW] == 11 ) vC -= vW; - if ( mesh->BCv.type[iVzE] == 11 ) vC -= vE; if ( Assemble == 1 ) { @@ -1577,196 +1345,13 @@ void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, Spars // StokesA->F[eqn] += (inE*mesh->sxz[ixyE] - inW*mesh->sxz[ixyW]) /dx; StokesA->F[eqn] += (mesh->szzd[iPrN] - mesh->szzd[iPrS])/dz; StokesA->F[eqn] -= (mesh->p_corr[iPrN] - mesh->p_corr[iPrS])/dz; - StokesA->F[eqn] += (mesh->sxz[ixyE] - mesh->sxz[ixyW]) /dx; + StokesA->F[eqn] += (mesh->sxz[ixyE] - mesh->sxz[ixyW]) /dx; StokesA->F[eqn] *= -1.0; StokesA->F[eqn] -= StokesA->b[eqn] - vC_corr*v[iVzC]; StokesA->F[eqn] *= celvol; } } - -#if 0 -// MD6 -void Zmomentum_InnerNodesDecoupled( SparseMat *Stokes, SparseMat *StokesA, SparseMat *StokesB, int Assemble, int lev, int stab, int comp, double om, int sign, params model, double one_dx, double one_dz, double one_dx_dx, double one_dz_dz, double one_dx_dz, double celvol, grid* mesh, int ith, int c1, int c2, int c3, int nx, int ncx, int nxvz, int eqn, double* u, double* v, double* p, int **JtempA, double **AtempA, int *nnzc2A, int **JtempB, double **AtempB, int *nnzc2B ) { - - double uSW=0.0, uSE=0.0, uNW=0.0, uNE=0.0, vS=0.0, vW=0.0, vC=0.0, vE=0.0, vN=0.0, pN=0.0, pS=0.0; - - // Coefficients - double D22S = mesh->D22_n[c2]; - double D22N = mesh->D22_n[c2+ncx]; - double D33E = mesh->D33_s[c1]; - double D33W = mesh->D33_s[c1-1]; - - if (mesh->BCp.type[c2+ncx] == 31) { - StokesA->b[eqn] += mesh->BCv.val[c3]*one_dz; - } - - double inS=0.0,inN=0.0,inWu=0.0,inWv=0.0,inEu=0.0,inEv=0.0; - - // dsyy/dz - if ( mesh->BCv.type[c3-nxvz] != 30 && mesh->BCv.type[c3+nxvz] !=30 ) { - inS = 1.0; - inN = 1.0; - } - if ( mesh->BCv.type[c3-nxvz] == 30 && mesh->BCv.type[c3+nxvz] != 30 ) { - inN = 1.0; - } - if ( mesh->BCv.type[c3-nxvz] != 30 && mesh->BCv.type[c3+nxvz] == 30 ) { - inS = 1.0; - } - - if (mesh->BCp.type[c2 ] == -1) inS = 1.0; - if (mesh->BCp.type[c2+ncx] == -1) inN = 1.0; - - double inW=0.0, inE = 0.0; - if (mesh->BCg.type[c1-1] != 30 && mesh->BCv.type[c3-1] != 13 ) inW = 1.0; - if (mesh->BCg.type[c1 ] != 30 && mesh->BCv.type[c3+1] != 13 ) inE = 1.0; - - double dx = mesh->dx; - double dz = mesh->dz; - - comp = 1.0; - // simpler - vW = -D33W*inW/pow(dx, 2); - vC = -(D22N*inN*((1.0/3.0)*comp/dz - 1/dz) - D22S*inS*(-1.0/3.0*comp/dz + 1.0/dz))/dz - (-D33E*inE/dx - D33W*inW/dx)/dx; - vE = -D33E*inE/pow(dx, 2); - vS = D22S*inS*((1.0/3.0)*comp/dz - 1/dz)/dz; - vN = -D22N*inN*(-1.0/3.0*comp/dz + 1.0/dz)/dz; - uSW = (1.0/3.0)*D22S*comp*inS/(dx*dz) - D33W*inW/(dx*dz); - uSE = -1.0/3.0*D22S*comp*inS/(dx*dz) + D33E*inE/(dx*dz); - uNW = -1.0/3.0*D22N*comp*inN/(dx*dz) + D33W*inW/(dx*dz); - uNE = (1.0/3.0)*D22N*comp*inN/(dx*dz) - D33E*inE/(dx*dz); - - if ( mesh->BCv.type[c3-1] == 11 ) vC -= vW; - if ( mesh->BCv.type[c3+1] == 11 ) vC -= vE; - - - // Stabilisation with density gradients - if (stab==1) { - double drhodz = (mesh->rho_n[c2+ncx] - mesh->rho_n[c2])*one_dz; - double drhodx = (mesh->rho_s[c1-1 ] - mesh->rho_s[c1])*one_dx; - double vC_corr = 1.00 * om * model.dt * mesh->gz[c3] * drhodz; - double vy_corr = 0.25 * om * model.dt * mesh->gz[c3] * drhodx; - -// if (mesh->rho_n[c2+ncx]*1e27<1000) { -// drhodz = 3*drhodz; -//// printf("mesh->rho_n[c2+ncx] = %2.2e\n", mesh->rho_n[c2+ncx]*1e27); -//// printf("mesh->rho_n[c2] = %2.2e\n", mesh->rho_n[c2]*1e27); -// } - // Importante trique, voire meme gigantesque! - if (vC+vC_corr>0.0) { - vC += vC_corr; -// uSW += vy_corr; -// uSE += vy_corr; -// uNW += vy_corr; -// uNE += vy_corr; -// if (vC_corr>1e-6)printf("vC_corr = %2.2e gz=%2.2e\n",vC_corr, mesh->gz[c3]); -// printf(mesh->rho_n[c2+ncx]) - - } - } - - // Pressure gradient - if ( mesh->BCp.type[c2] != 30 && mesh->BCp.type[c2+ncx] != 30 ) { - pS = -one_dz * 1e0; - pN = one_dz * 1e0; - } - -// uSW=-uSW, uSE=-uSE, uNW=-uNW, uNE=-uNE, vS=-vS, vW=-vW, vC=-vC, vE=-vE, vN=-vN, pN=-pN, pS=-pS; - - if ( Assemble == 1 ) { - StokesA->b[eqn] *= celvol; - StokesB->b[eqn] *= celvol; - //-------------------- -// if ( mesh->BCu.type[c1-1] != 30 && mesh->BCu.type[c1] != 30 ) { -// AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1-1], &(nnzc2A[ith]), uSW*celvol, mesh->BCu.type[c1-1], mesh->BCu.val[c1-1], StokesA->bbc ); -// AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1], &(nnzc2A[ith]), uSE*celvol, mesh->BCu.type[c1], mesh->BCu.val[c1], StokesA->bbc ); -// } -// if ( mesh->BCu.type[c1+nx-1] != 30 && mesh->BCu.type[c1+nx] != 30) { -// AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx-1], &(nnzc2A[ith]), uNW*celvol, mesh->BCu.type[c1+nx-1], mesh->BCu.val[c1+nx-1], StokesA->bbc ); -// AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx], &(nnzc2A[ith]), uNE*celvol, mesh->BCu.type[c1+nx], mesh->BCu.val[c1+nx], StokesA->bbc ); -// } - -// if ( mesh->BCu.type[c1+nx] != 30 && mesh->BCu.type[c1] != 30 ) { -// AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx], &(nnzc2A[ith]), uNE*celvol, mesh->BCu.type[c1+nx], mesh->BCu.val[c1+nx], StokesA->bbc ); -// AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1], &(nnzc2A[ith]), uSE*celvol, mesh->BCu.type[c1], mesh->BCu.val[c1], StokesA->bbc ); -// } -// if ( mesh->BCu.type[c1+nx-1] != 30 && mesh->BCu.type[c1-1] != 30) { -// AddCoeff3( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx-1], &(nnzc2A[ith]), uNW*celvol, mesh->BCu.type[c1+nx-1], mesh->BCu.val[c1+nx-1], StokesA->bbc ); -// AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1-1], &(nnzc2A[ith]), uSW*celvol, mesh->BCu.type[c1-1], mesh->BCu.val[c1-1], StokesA->bbc ); -// } -// - - if ( mesh->BCu.type[c1+nx] != 30 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx], &(nnzc2A[ith]), uNE*celvol, mesh->BCu.type[c1+nx], mesh->BCu.val[c1+nx], StokesA->bbc ); - if ( mesh->BCu.type[c1] != 30 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1], &(nnzc2A[ith]), uSE*celvol, mesh->BCu.type[c1], mesh->BCu.val[c1], StokesA->bbc ); - if ( mesh->BCu.type[c1+nx-1] != 30 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1+nx-1], &(nnzc2A[ith]), uNW*celvol, mesh->BCu.type[c1+nx-1], mesh->BCu.val[c1+nx-1], StokesA->bbc ); - if ( mesh->BCu.type[c1-1] != 30 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_u[c1-1], &(nnzc2A[ith]), uSW*celvol, mesh->BCu.type[c1-1], mesh->BCu.val[c1-1], StokesA->bbc ); - - //-------------------- - if ( mesh->BCv.type[c3-nxvz] != 30 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-nxvz], &(nnzc2A[ith]), vS*celvol, mesh->BCv.type[c3-nxvz], mesh->BCv.val[c3-nxvz], StokesA->bbc ); - if ( mesh->BCv.type[c3-1] != 30 && mesh->BCv.type[c3-1] != -12 ) { - if( mesh->BCv.type[c3-1] != 11 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-1], &(nnzc2A[ith]), vW*celvol, mesh->BCv.type[c3-1], mesh->BCv.val[c3-1], StokesA->bbc ); - else AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-1], &(nnzc2A[ith]), vW*celvol, mesh->BCv.type[c3-1], 2.0*mesh->BCv.val[c3-1], StokesA->bbc ); - } - // Periodic - if ( mesh->BCv.type[c3-1] == -12 ) { - AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+nxvz-3], &(nnzc2A[ith]), vW*celvol, mesh->BCv.type[c3+nxvz-3], mesh->BCv.val[c3+nxvz-3], StokesA->bbc ); - } - AddCoeff2( JtempA[ith], AtempA[ith], eqn, eqn, &(nnzc2A[ith]), vC*celvol, mesh->BCv.type[c3], mesh->BCv.val[c3], StokesA->bbc ); - if ( mesh->BCv.type[c3+1] != 30 && mesh->BCv.type[c3+1] != -12 ) { - if ( mesh->BCv.type[c3+1] != 11 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+1], &(nnzc2A[ith]), vE*celvol, mesh->BCv.type[c3+1], mesh->BCv.val[c3+1], StokesA->bbc ); - else AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+1], &(nnzc2A[ith]), vE*celvol, mesh->BCv.type[c3+1], 2.0*mesh->BCv.val[c3+1], StokesA->bbc ); - } - // Periodic - if ( mesh->BCv.type[c3+1] == -12 ) { - AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3-nxvz+3], &(nnzc2A[ith]), vE*celvol, mesh->BCv.type[c3-nxvz+3], mesh->BCv.val[c3-nxvz+3], StokesA->bbc ); - } - - if ( mesh->BCv.type[c3+nxvz] != 30 ) AddCoeff2( JtempA[ith], AtempA[ith], eqn, Stokes->eqn_v[c3+nxvz], &(nnzc2A[ith]), vN*celvol, mesh->BCv.type[c3+nxvz], mesh->BCv.val[c3+nxvz], StokesA->bbc ); - //-------------------- - if ( mesh->BCp.type[c2] != 30 && mesh->BCp.type[c2+ncx] != 30) { - AddCoeff2( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[c2] - Stokes->neq_mom, &(nnzc2B[ith]), pS*celvol, mesh->BCp.type[c2], mesh->BCp.val[c2], StokesB->bbc ); - AddCoeff2( JtempB[ith], AtempB[ith], eqn, Stokes->eqn_p[c2+ncx] - Stokes->neq_mom, &(nnzc2B[ith]), pN*celvol, mesh->BCp.type[c2+ncx], mesh->BCp.val[c2+ncx], StokesB->bbc ); - } - - } - else { - - // Residual - StokesA->F[eqn] = vC*v[c3]; - if ( mesh->BCp.type[c2] != 30 && mesh->BCp.type[c2+ncx] != 30 ) { - StokesA->F[eqn] += pS*p[c2] + pN*p[c2+ncx]; - } - if ( mesh->BCv.type[c3-nxvz] != 30 ) StokesA->F[eqn] += vS*v[c3-nxvz]; - if ( mesh->BCv.type[c3-1] != 30 && mesh->BCv.type[c3-1] != 11 && mesh->BCv.type[c3-1] != -12 ) StokesA->F[eqn] += vW*v[c3-1]; - if ( mesh->BCv.type[c3+1] != 30 && mesh->BCv.type[c3+1] != 11 && mesh->BCv.type[c3+1] != -12 ) StokesA->F[eqn] += vE*v[c3+1]; - if ( mesh->BCv.type[c3+nxvz] != 30 ) StokesA->F[eqn] += vN*v[c3+nxvz]; -// if ( mesh->BCu.type[c1-1] != 30 && mesh->BCu.type[c1] != 30 ) { -// StokesA->F[eqn] += uSW*u[c1-1] + uSE*u[c1]; -// } -// if ( mesh->BCu.type[c1+nx-1] != 30 && mesh->BCu.type[c1+nx] != 30 ) { -// StokesA->F[eqn] += uNW*u[c1+nx-1] + uNE*u[c1+nx]; -// } - - if ( mesh->BCu.type[c1-1] != 30 ) StokesA->F[eqn] += uSW*u[c1-1]; - if ( mesh->BCu.type[c1] != 30 ) StokesA->F[eqn] += uSE*u[c1]; - - if ( mesh->BCu.type[c1+nx-1] != 30 ) StokesA->F[eqn] += uNW*u[c1+nx-1]; - if ( mesh->BCu.type[c1+nx] != 30 ) StokesA->F[eqn] += uNE*u[c1+nx]; - - - if ( mesh->BCv.type[c3-1] == 11 ) StokesA->F[eqn] += -2*D33W*one_dx_dx*mesh->BCv.val[c3-1]; - if ( mesh->BCv.type[c3+1] == 11 ) StokesA->F[eqn] += -2*D33E*one_dx_dx*mesh->BCv.val[c3+1]; - if ( mesh->BCv.type[c3-1] == -12 ) StokesA->F[eqn] += vW*v[c3+nxvz-3]; - if ( mesh->BCv.type[c3+1] == -12 ) StokesA->F[eqn] += vE*v[c3-nxvz+3]; - StokesA->F[eqn] -= (StokesA->b[eqn]); - StokesA->F[eqn] *= celvol; -// printf("StokesA->F[eqn] = %2.2e %2.2e %2.2e %2.2e %2.2e\n", StokesA->F[eqn], vC*v[c3], vS*v[c3-nxvz], vW*v[c3-1], StokesA->b[eqn]); - } -} -#endif - - /*--------------------------------------------------------------------------------------------------------------------*/ /*------------------------------------------------------ M-Doodz -----------------------------------------------------*/ /*--------------------------------------------------------------------------------------------------------------------*/ @@ -2225,7 +1810,6 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ // Matrix initialisation if ( Assemble == 1 ) { - nnzA = 9*((mesh->Nx-1) * mesh->Nz + (mesh->Nz-1) * mesh->Nx); nnzB = 5*((mesh->Nx-1) * (mesh->Nz-1)); nnzC = nnzB; @@ -2239,7 +1823,6 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ AllocMat( StokesB, nnzB ); AllocMat( StokesC, nnzC ); AllocMat( StokesD, nnzD ); - } // Build velocity block RHS @@ -2332,14 +1915,11 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ if ( k==nx-1 && mesh->BCu.type[c1]==2 ) { Xmomentum_EastNeumannDecoupled( Stokes, StokesA, StokesB, Assemble, lev, stab, comp, theta, sign, model, one_dx, one_dz, one_dx_dx, one_dz_dz, one_dx_dz, celvol, mesh, ith, c1, c2, c3, nx, ncx, nxvz, eqn, u, v, p_corr, JtempA, AtempA, nnzc2A, JtempB, AtempB, nnzc2B ); } - // SUM+=StokesA->F[eqn];//*StokesA->F[eqn]; - } } } - // printf("%2.8e\n", SUM); - + //-----------------------------------------------------------------// if ( Assemble == 1 ) { @@ -2391,7 +1971,6 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ //--------------------- INNER NODES ---------------------// if ( k>0 && k0 && lb[k] = StokesC->b[k]; } - - // Final index StokesA->Ic[StokesA->neq] = nnzcA; StokesB->Ic[StokesB->neq] = nnzcB; @@ -2650,7 +2225,6 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ //--------------------------------------// -#ifndef _VG_ if ( model.writer_debug == 1 ) { char *filename; @@ -2659,33 +2233,33 @@ void BuildStokesOperatorDecoupled( grid *mesh, params model, int lev, double *p_ // Fill in DD data structure OutputSparseMatrix OutputDDA, OutputDDB, OutputDDC, OutputDDD; - - OutputDDA.V = StokesA->A; - OutputDDA.Ic = StokesA->Ic; - OutputDDA.J = StokesA->J; - OutputDDA.b = StokesA->b; - OutputDDA.F = StokesA->F; - OutputDDB.V = StokesB->A; - OutputDDB.Ic = StokesB->Ic; - OutputDDB.J = StokesB->J; - OutputDDB.b = StokesB->b; - OutputDDC.V = StokesC->A; - OutputDDC.Ic = StokesC->Ic; - OutputDDC.J = StokesC->J; - OutputDDC.b = StokesC->b; - OutputDDC.F = StokesC->F; - OutputDDD.V = StokesD->A; - OutputDDD.Ic = StokesD->Ic; - OutputDDD.J = StokesD->J; - OutputDDD.b = StokesD->b; - OutputDDA.eta_cell = mesh->eta_n; + + OutputDDA.V = StokesA->A; + OutputDDA.Ic = StokesA->Ic; + OutputDDA.J = StokesA->J; + OutputDDA.b = StokesA->b; + OutputDDA.F = StokesA->F; + OutputDDB.V = StokesB->A; + OutputDDB.Ic = StokesB->Ic; + OutputDDB.J = StokesB->J; + OutputDDB.b = StokesB->b; + OutputDDC.V = StokesC->A; + OutputDDC.Ic = StokesC->Ic; + OutputDDC.J = StokesC->J; + OutputDDC.b = StokesC->b; + OutputDDC.F = StokesC->F; + OutputDDD.V = StokesD->A; + OutputDDD.Ic = StokesD->Ic; + OutputDDD.J = StokesD->J; + OutputDDD.b = StokesD->b; + OutputDDA.eta_cell = mesh->eta_n; OutputDDA.params[0] = nx; OutputDDA.params[1] = nz; OutputDDA.params[2] = mesh->dx; OutputDDA.params[3] = mesh->dz; - OutputDDA.eqn_u = DoodzMalloc(nx*nzvx*sizeof(int)); - OutputDDA.eqn_v = DoodzMalloc(nxvz*nz*sizeof(int)); - OutputDDA.eqn_p = DoodzMalloc(ncx*ncz*sizeof(int)); + OutputDDA.eqn_u = DoodzMalloc(nx*nzvx*sizeof(int)); + OutputDDA.eqn_v = DoodzMalloc(nxvz*nz*sizeof(int)); + OutputDDA.eqn_p = DoodzMalloc(ncx*ncz*sizeof(int)); for( l=0; lBCu.type[i+j*(nx)]); -// } -// printf("\n"); -// } - - // Vx Dirichlet for( i=0; iF, StokesC->F, dx, *model, mesh, scaling, Stokes ); if ( model->lin_solver == 0 ) DirectStokesDecoupled ( StokesA, StokesB, StokesC, StokesD, PardisoStokes, StokesA->F, StokesC->F, dx, *model, mesh, scaling, Stokes ); if ( model->lin_solver == 1 ) KSPStokesDecoupled ( StokesA, StokesB, StokesC, StokesD, PardisoStokes, StokesA->F, StokesC->F, dx, *model, mesh, scaling, Stokes, Stokes, JacobA, JacobB, JacobC ); if ( model->lin_solver == 2 ) KillerSolver ( StokesA, StokesB, StokesC, StokesD, PardisoStokes, StokesA->F, StokesC->F, dx, *model, mesh, scaling, Stokes, Stokes, JacobA, JacobB, JacobC ); diff --git a/SETS/ViscousInclusion.c b/SETS/ViscousInclusion.c index 8e7d350b..f04aad90 100755 --- a/SETS/ViscousInclusion.c +++ b/SETS/ViscousInclusion.c @@ -89,7 +89,15 @@ SetBC SetBCVx(MdoodzInput *instance, POSITION position, Coordinates coord) { const double mc = 1e3; double Vx, Vz, P, eta, sxx, szz, x, z; if (instance->model.shear_style == 0) { - if (position == W || position == E) { + if (position == W) { + x = coord.x; + z = coord.z; + eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 1, radius, mm, mc); + // bc.type = 0; + // bc.value = Vx; + bc.type = 2; + bc.value = sxx; + } else if (position == E) { x = coord.x; z = coord.z; eval_anal_Dani(&Vx, &Vz, &P, &eta, &sxx, &szz, x, z, 1, radius, mm, mc); diff --git a/SETS/ViscousInclusion.txt b/SETS/ViscousInclusion.txt index 0302ab9d..cb7b8c77 100755 --- a/SETS/ViscousInclusion.txt +++ b/SETS/ViscousInclusion.txt @@ -37,39 +37,39 @@ RK = 4 /**** LINEAR SOLVER ****/ noisy = 0 penalty = 1e5 -lin_solver = -1 +lin_solver = -1 diag_scaling = 0 lin_abs_div = 1e-8 lin_rel_div = 1e-8 /**** NON-LINEAR SOLVER ****/ -Newton = 0 -line_search = 0 -nit_max = 100 -rel_tol_KSP = 1e-4 -nonlin_abs_mom = 1e-8 -nonlin_abs_mom = 1e-8 -nonlin_rel_div = 1e-8 -nonlin_rel_div = 1e-8 +Newton = 0 +line_search = 0 +nit_max = 100 +rel_tol_KSP = 1e-4 +nonlin_abs_mom = 1e-8 +nonlin_abs_mom = 1e-8 +nonlin_rel_div = 1e-8 +nonlin_rel_div = 1e-8 /**** VISCOSITY CUT-OFF ****/ min_eta = 1e-3 max_eta = 1e6 /**** SWITCHES ****/ -elastic = 0 -free_surface = 0 -free_surface_stab = 0 -finite_strain = 1 +elastic = 0 +free_surface = 0 +free_surface_stab = 0 +finite_strain = 1 advection = 0 -gnuplot_log_res = 0 -eta_average = 0 +gnuplot_log_res = 0 +eta_average = 0 /**** SETUP-DEPENDANT ****/ shear_style = 0 -periodic_x = 0 -pure_shear_ALE = 1 -bkg_strain_rate = 1 +periodic_x = 0 +pure_shear_ALE = 1 +bkg_strain_rate = 1 user0 = 0 user1 = 0.15 / inclusion radius [m] user2 = 0 @@ -83,7 +83,7 @@ gz = 0.000 Nb_phases = 2 /**** PHASE 1 ****/ -ID = 0 +ID = 0 plast = 0 cstv = 1 / constant visc law eta0 = 1 @@ -91,7 +91,7 @@ npwl = 1.0 Qpwl = 0 /**** PHASE 2 ****/ -ID = 1 +ID = 1 plast = 0 cstv = 1 / constant visc law eta0 = 1000 diff --git a/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb b/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb index aa5271b8..ae5d7840 100644 --- a/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb +++ b/misc/Python_CodeGen/AssembleGeneralStiffness_MDOODZ_7.0_v4_StressBC.ipynb @@ -2,17 +2,17 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Txx = D11*Exx + D12*Ezz + D13*Gxz + D14*divV;\n", - "Tzz = D21*Exx + D22*Ezz + D23*Gxz + D24*divV;\n", - "Txz = D31*Exx + D32*Ezz + D33*Gxz + D34*divV;\n", - "P = D41*Exx + D42*Ezz + D43*Gxz + D44*divV;\n" + "Txx = D11*Exx;\n", + "Tzz = D22*Ezz;\n", + "Txz = D33*Gxz;\n", + "P = 0;\n" ] } ], @@ -48,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 48, "metadata": { "scrolled": false }, @@ -57,39 +57,42 @@ "name": "stdout", "output_type": "stream", "text": [ - "uW = (1.0/3.0)*(-0.75*D13W*dx*(pow(inNWv, 2) - pow(inSWv, 2)) + 0.25*dx*pow(inW, 2)*(-D31N*(comp*inW*out_of_plane - 3) + D31S*(comp*inW*out_of_plane - 3) - D32N*comp*inW*out_of_plane + D32S*comp*inW*out_of_plane) + dz*inW*(D11W*(comp*inW*out_of_plane - 3) + D12W*comp*inW*out_of_plane))/(pow(dx, 2)*dz);\n", - "uC = (dx*(dx*(D33N*inN + D33S*inS) + dz*(-D31N*(pow(inE, 2)*(0.083333333333333329*comp*inE*out_of_plane - 0.25) + pow(inW, 2)*(-0.083333333333333329*comp*inW*out_of_plane + 0.25)) + D31S*(pow(inE, 2)*(0.083333333333333329*comp*inE*out_of_plane - 0.25) + pow(inW, 2)*(-0.083333333333333329*comp*inW*out_of_plane + 0.25)) - 0.083333333333333329*D32N*comp*out_of_plane*(pow(inE, 3) - pow(inW, 3)) + 0.083333333333333329*D32S*comp*out_of_plane*(pow(inE, 3) - pow(inW, 3)))) - 1.0/3.0*dz*(0.75*dx*(-D13E + D13W)*(pow(inN, 2) - pow(inS, 2)) + dz*(D11E*inE*(comp*inE*out_of_plane - 3) + D11W*inW*(comp*inW*out_of_plane - 3) + D12E*comp*pow(inE, 2)*out_of_plane + D12W*comp*pow(inW, 2)*out_of_plane)))/(pow(dx, 2)*pow(dz, 2));\n", - "uE = (1.0/3.0)*(0.75*D13E*dx*(pow(inNEv, 2) - pow(inSEv, 2)) + 0.25*dx*pow(inE, 2)*(D31N*(comp*inE*out_of_plane - 3) - D31S*(comp*inE*out_of_plane - 3) + D32N*comp*inE*out_of_plane - D32S*comp*inE*out_of_plane) + dz*inE*(D11E*(comp*inE*out_of_plane - 3) + D12E*comp*inE*out_of_plane))/(pow(dx, 2)*dz);\n", - "uS = (-D33S*dx*inS + 0.25*dz*pow(inS, 2)*(D13E - D13W) + dz*(D31S*(pow(inSEc, 2)*(0.083333333333333329*comp*inSEc*out_of_plane - 0.25) + pow(inSWc, 2)*(-0.083333333333333329*comp*inSWc*out_of_plane + 0.25)) + 0.083333333333333329*D32S*comp*out_of_plane*(pow(inSEc, 3) - pow(inSWc, 3))))/(dx*pow(dz, 2));\n", - "uN = (-D33N*dx*inN + 0.25*dz*pow(inN, 2)*(-D13E + D13W) - dz*(D31N*(pow(inNEc, 2)*(0.083333333333333329*comp*inNEc*out_of_plane - 0.25) + pow(inNWc, 2)*(-0.083333333333333329*comp*inNWc*out_of_plane + 0.25)) + 0.083333333333333329*D32N*comp*out_of_plane*(pow(inNEc, 3) - pow(inNWc, 3))))/(dx*pow(dz, 2));\n", - "vSW = (-dx*(D33S*dz*inS + dx*(0.083333333333333329*D31N*comp*pow(inW, 3)*out_of_plane + 0.083333333333333329*D31S*comp*out_of_plane*(pow(inSWc, 3) - pow(inW, 3)) + 0.083333333333333329*D32N*pow(inW, 2)*(comp*inW*out_of_plane - 3) + D32S*(pow(inSWc, 2)*(0.083333333333333329*comp*inSWc*out_of_plane - 0.25) + pow(inW, 2)*(-0.083333333333333329*comp*inW*out_of_plane + 0.25)))) + (1.0/3.0)*dz*(dx*inW*(D11W*comp*inW*out_of_plane + D12W*(comp*inW*out_of_plane - 3)) + 0.75*dz*(D13E*pow(inS, 2) - D13W*(pow(inS, 2) - pow(inSWv, 2)))))/(pow(dx, 2)*pow(dz, 2));\n", - "vSE = (dx*(D33S*dz*inS + dx*(-0.083333333333333329*D31N*comp*pow(inE, 3)*out_of_plane + 0.083333333333333329*D31S*comp*out_of_plane*(pow(inE, 3) - pow(inSEc, 3)) - 0.083333333333333329*D32N*pow(inE, 2)*(comp*inE*out_of_plane - 3) + D32S*(pow(inE, 2)*(0.083333333333333329*comp*inE*out_of_plane - 0.25) + pow(inSEc, 2)*(-0.083333333333333329*comp*inSEc*out_of_plane + 0.25)))) - 1.0/3.0*dz*(dx*inE*(D11E*comp*inE*out_of_plane + D12E*(comp*inE*out_of_plane - 3)) + 0.75*dz*(D13E*(pow(inS, 2) - pow(inSEv, 2)) - D13W*pow(inS, 2))))/(pow(dx, 2)*pow(dz, 2));\n", - "vNW = (dx*(D33N*dz*inN - dx*(0.083333333333333329*D31N*comp*out_of_plane*(pow(inNWc, 3) - pow(inW, 3)) + 0.083333333333333329*D31S*comp*pow(inW, 3)*out_of_plane + D32N*(pow(inNWc, 2)*(0.083333333333333329*comp*inNWc*out_of_plane - 0.25) + pow(inW, 2)*(-0.083333333333333329*comp*inW*out_of_plane + 0.25)) + 0.083333333333333329*D32S*pow(inW, 2)*(comp*inW*out_of_plane - 3))) - 1.0/3.0*dz*(dx*inW*(D11W*comp*inW*out_of_plane + D12W*(comp*inW*out_of_plane - 3)) + 0.75*dz*(-D13E*pow(inN, 2) + D13W*(pow(inN, 2) - pow(inNWv, 2)))))/(pow(dx, 2)*pow(dz, 2));\n", - "vNE = (-dx*(D33N*dz*inN + dx*(-0.083333333333333329*D31N*comp*out_of_plane*(pow(inE, 3) - pow(inNEc, 3)) + 0.083333333333333329*D31S*comp*pow(inE, 3)*out_of_plane - D32N*(pow(inE, 2)*(0.083333333333333329*comp*inE*out_of_plane - 0.25) + pow(inNEc, 2)*(-0.083333333333333329*comp*inNEc*out_of_plane + 0.25)) + 0.083333333333333329*D32S*pow(inE, 2)*(comp*inE*out_of_plane - 3))) + (1.0/3.0)*dz*(dx*inE*(D11E*comp*inE*out_of_plane + D12E*(comp*inE*out_of_plane - 3)) + 0.75*dz*(-D13E*(pow(inN, 2) - pow(inNEv, 2)) + D13W*pow(inN, 2))))/(pow(dx, 2)*pow(dz, 2));\n", - "uSW = (-0.25*D13W*pow(inSWv, 2) + 0.083333333333333329*D31S*pow(inSWc, 2)*(comp*inSWc*out_of_plane - 3) + 0.083333333333333329*D32S*comp*pow(inSWc, 3)*out_of_plane)/(dx*dz);\n", - "uSE = (0.25*D13E*pow(inSEv, 2) - 0.083333333333333329*D31S*pow(inSEc, 2)*(comp*inSEc*out_of_plane - 3) - 0.083333333333333329*D32S*comp*pow(inSEc, 3)*out_of_plane)/(dx*dz);\n", - "uNW = (0.25*D13W*pow(inNWv, 2) - 0.083333333333333329*D31N*pow(inNWc, 2)*(comp*inNWc*out_of_plane - 3) - 0.083333333333333329*D32N*comp*pow(inNWc, 3)*out_of_plane)/(dx*dz);\n", - "uNE = (-0.25*D13E*pow(inNEv, 2) + 0.083333333333333329*D31N*pow(inNEc, 2)*(comp*inNEc*out_of_plane - 3) + 0.083333333333333329*D32N*comp*pow(inNEc, 3)*out_of_plane)/(dx*dz);\n", - "vSWW = -0.25*D13W*pow(inSWv, 2)/pow(dx, 2);\n", - "vSEE = -0.25*D13E*pow(inSEv, 2)/pow(dx, 2);\n", - "vNWW = -0.25*D13W*pow(inNWv, 2)/pow(dx, 2);\n", - "vNEE = -0.25*D13E*pow(inNEv, 2)/pow(dx, 2);\n", - "vSSW = 0.083333333333333329*pow(inSWc, 2)*(D31S*comp*inSWc*out_of_plane + D32S*(comp*inSWc*out_of_plane - 3))/pow(dz, 2);\n", - "vSSE = 0.083333333333333329*pow(inSEc, 2)*(D31S*comp*inSEc*out_of_plane + D32S*(comp*inSEc*out_of_plane - 3))/pow(dz, 2);\n", - "vNNW = 0.083333333333333329*pow(inNWc, 2)*(D31N*comp*inNWc*out_of_plane + D32N*(comp*inNWc*out_of_plane - 3))/pow(dz, 2);\n", - "vNNE = 0.083333333333333329*pow(inNEc, 2)*(D31N*comp*inNEc*out_of_plane + D32N*(comp*inNEc*out_of_plane - 3))/pow(dz, 2);\n", - "pW = (0.25*dx*inW*(-D34N + D34S) + dz*(D14W - 1))/(dx*dz);\n", - "pE = (0.25*dx*inE*(-D34N + D34S) + dz*(1 - D14E))/(dx*dz);\n", - "pSW = 0.25*D34S*inSWc/dz;\n", - "pSE = 0.25*D34S*inSEc/dz;\n", - "pNW = -0.25*D34N*inSWc/dz;\n", - "pNE = -0.25*D34N*inSEc/dz;\n" + "uW = D11W*inW*(0.33333333333333331*inW*out_of_plane - 1)/pow(dx, 2);\n", + "uC = (pow(dx, 2)*(D33N*(2*DirS + inN) + D33S*(2*DirN + inS)) - pow(dz, 2)*(D11E*inE*(0.33333333333333331*inE*out_of_plane - 1) + D11W*inW*(0.33333333333333331*inW*out_of_plane - 1)))/(pow(dx, 2)*pow(dz, 2));\n", + "uE = D11E*inE*(0.33333333333333331*inE*out_of_plane - 1)/pow(dx, 2);\n", + "uS = -D33S*inS/pow(dz, 2);\n", + "uN = -D33N*inN/pow(dz, 2);\n", + "vSW = (0.33333333333333331*D11W*pow(inW, 2)*out_of_plane - D33S)/(dx*dz);\n", + "vSE = (-0.33333333333333331*D11E*pow(inE, 2)*out_of_plane + D33S)/(dx*dz);\n", + "vNW = (-0.33333333333333331*D11W*pow(inW, 2)*out_of_plane + D33N)/(dx*dz);\n", + "vNE = (0.33333333333333331*D11E*pow(inE, 2)*out_of_plane - D33N)/(dx*dz);\n", + "uSW = 0;\n", + "uSE = 0;\n", + "uNW = 0;\n", + "uNE = 0;\n", + "vSWW = 0;\n", + "vSEE = 0;\n", + "vNWW = 0;\n", + "vNEE = 0;\n", + "vSSW = 0;\n", + "vSSE = 0;\n", + "vNNW = 0;\n", + "vNNE = 0;\n", + "pW = -1/dx;\n", + "pE = 1.0/dx;\n", + "pSW = 0;\n", + "pSE = 0;\n", + "pNW = 0;\n", + "pNE = 0;\n" ] } ], "source": [ "# Compressibility switch\n", "comp = symbols('comp')\n", + "comp = 1\n", + "\n", + "SxxBC = symbols('SxxBC') # boundary stress\n", "\n", "# Stencil x\n", "uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE = symbols('uW,uC,uE,uS,uN,vSW,vSE,vNW,vNE')\n", @@ -109,6 +112,14 @@ "inS,inN,inW,inE = symbols('inS,inN,inW,inE')\n", "inSEc,inSWc,inNEc,inNWc = symbols('inSEc,inSWc,inNEc,inNWc')\n", "inSEv,inSWv,inNEv,inNWv = symbols('inSEv,inSWv,inNEv,inNWv')\n", + "NeuN, NeuS, DirN, DirS = symbols('NeuN, NeuS, DirN, DirS')\n", + "\n", + "dVxdzN_in = (uN - uC )/dz\n", + "dVxdzS_in = (uC - uS )/dz\n", + "dVxdzN_Dir = ((-uC) - uC )/dz\n", + "dVxdzS_Dir = (uC - (-uC) )/dz\n", + "dVxdzN_Neu = ((uC) - uC )/dz\n", + "dVxdzS_Neu = (uC - (uC) )/dz\n", "\n", "# Divergences\n", "divW = inW *out_of_plane*( (uC - uW )/dx + (vNW - vSW )/dz ) \n", @@ -131,9 +142,8 @@ "GxzSW = inSWv*((uW - uSW)/dz + ( vSW - vSWW )/dx)\n", "\n", "# Shear strain rate\n", - "GxzN = inN *((uN - uC )/dz + (vNE -vNW )/dx) # clear\n", - "GxzS = inS *((uC - uS )/dz + (vSE -vSW )/dx) \n", - "\n", + "GxzN = ( inN*dVxdzN_in + DirS*dVxdzN_Dir + NeuS*dVxdzN_Neu + (vNE -vNW )/dx) \n", + "GxzS = ( inS*dVxdzS_in + DirN*dVxdzS_Dir + NeuN*dVxdzS_Neu + (vSE -vSW )/dx) \n", "ExxNE = inNEc*((uNE -uN )/dx - comp*1/3*divNE)\n", "ExxNW = inNWc*((uN -uNW)/dx - comp*1/3*divNW)\n", "ExxSE = inSEc*((uSE -uS )/dx - comp*1/3*divSE)\n", @@ -169,12 +179,16 @@ "TN = Dv*EdN\n", "TxzN = TN[2].subs({D31:'D31N', D32:'D32N', D33:'D33N', D34:'D34N'})\n", "#----------------------------------------#\n", - "# Fx = 1/dx*(TxxE - TxxW) + 1/dz*(TxzN - TxzS) - 1/dx*(pE - pW)\n", + "Fx = 1/dx*(TxxE - TxxW) + 1/dz*(TxzN - TxzS) - 1/dx*(pE - pW)\n", + "Fx = -Fx\n", "\n", "# Stress boundary\n", + "SxxE = TxxE - pE\n", "SxxW = 2*SxxBC - SxxE\n", - "Fx = 1/dx*(SxxE - SxxW) + 1/dz*(TxzN - TxzS)\n", - "Fx = -Fx\n", + "Fx = 1/dx*(SxxE - SxxW) + 1/dz*(TxzN - TxzS)*2 # dodgy factor 2 make it symmetric\n", + "Fx = -Fx/2\n", + "\n", + "\n", "for i in range(len(dofs)):\n", " cUc = Fx.diff(dofs[i])\n", " print(str(dofs[i]) + ' = ' + ccode(simplify(cUc)) + ';')" @@ -182,7 +196,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 46, "metadata": {}, "outputs": [ { @@ -190,14 +204,14 @@ "output_type": "stream", "text": [ "vW = -D33W*inW/pow(dx, 2);\n", - "vC = (1.0/3.0)*(-pow(dx, 2)*(D22N*inN + D22S*inS)*(comp*out_of_plane - 3) + 3*pow(dz, 2)*(D33E*inE + D33W*inW))/(pow(dx, 2)*pow(dz, 2));\n", + "vC = (-pow(dx, 2)*(D22N + D22S)*(0.33333333333333331*out_of_plane - 1) + pow(dz, 2)*(D33E*(2*DirE + inE) + D33W*(2*DirW + inW)))/(pow(dx, 2)*pow(dz, 2));\n", "vE = -D33E*inE/pow(dx, 2);\n", - "vS = (1.0/3.0)*D22S*inS*(comp*out_of_plane - 3)/pow(dz, 2);\n", - "vN = (1.0/3.0)*D22N*inN*(comp*out_of_plane - 3)/pow(dz, 2);\n", - "uSW = ((1.0/3.0)*D22S*comp*inS*out_of_plane - D33W*inW)/(dx*dz);\n", - "uSE = (-1.0/3.0*D22S*comp*inS*out_of_plane + D33E*inE)/(dx*dz);\n", - "uNW = (-1.0/3.0*D22N*comp*inN*out_of_plane + D33W*inW)/(dx*dz);\n", - "uNE = ((1.0/3.0)*D22N*comp*inN*out_of_plane - D33E*inE)/(dx*dz);\n", + "vS = D22S*(0.33333333333333331*out_of_plane - 1)/pow(dz, 2);\n", + "vN = D22N*(0.33333333333333331*out_of_plane - 1)/pow(dz, 2);\n", + "uSW = (0.33333333333333331*D22S*out_of_plane - D33W)/(dx*dz);\n", + "uSE = (-0.33333333333333331*D22S*out_of_plane + D33E)/(dx*dz);\n", + "uNW = (-0.33333333333333331*D22N*out_of_plane + D33W)/(dx*dz);\n", + "uNE = (0.33333333333333331*D22N*out_of_plane - D33E)/(dx*dz);\n", "vSW = 0;\n", "vSE = 0;\n", "vNW = 0;\n", @@ -210,8 +224,8 @@ "uSSE = 0;\n", "uNNW = 0;\n", "uNNE = 0;\n", - "pS = -inS/dz;\n", - "pN = inN/dz;\n", + "pS = -1/dz;\n", + "pN = 1.0/dz;\n", "pSW = 0;\n", "pSE = 0;\n", "pNW = 0;\n", @@ -222,6 +236,7 @@ "source": [ "# Compressibility switch\n", "comp = symbols('comp')\n", + "comp = 1\n", "\n", "# Stencil z\n", "vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE = symbols('vW,vC,vE,vS,vN,uSW,uSE,uNW,uNE')\n", @@ -241,6 +256,7 @@ "inS,inN,inW,inE = symbols('inS,inN,inW,inE')\n", "inSEc,inSWc,inNEc,inNWc = symbols('inSEc,inSWc,inNEc,inNWc')\n", "inSEv,inSWv,inNEv,inNWv = symbols('inSEv,inSWv,inNEv,inNWv')\n", + "DirE, DirW, NeuE, NeuW = symbols('DirE, DirW, NeuE, NeuW')\n", "\n", "# Divergences\n", "divS = out_of_plane*( (uSE-uSW)/dx + (vC -vS )/dz ) \n", @@ -263,8 +279,14 @@ "GxzNE = ((uNNE - uNE )/dz + (vNE - vN )/dx)\n", "\n", "# Shear strain rate\n", - "GxzE = (uNE-uSE)/dz + (vE - vC )/dx # clear\n", - "GxzW = (uNW-uSW)/dz + (vC - vW )/dx\n", + "dVzdxE_in = (vE - vC )/dx\n", + "dVzdxW_in = (vC - vW )/dx\n", + "dVzdxE_Dir = (-vC - vC )/dx\n", + "dVzdxW_Dir = (vC - -vC )/dx\n", + "dVzdxE_Neu = (vC - vC )/dx\n", + "dVzdxW_Neu = (vC - vC )/dx\n", + "GxzE = (uNE-uSE)/dz + inE*dVzdxE_in + DirE*dVzdxE_Dir + NeuE*dVzdxE_Neu\n", + "GxzW = (uNW-uSW)/dz + inW*dVzdxW_in + DirW*dVzdxW_Dir + NeuW*dVzdxW_Neu\n", "\n", "# Additional missing stress tensor components for interpolation\n", "ExxNE = ((uNEE-uNE)/dx - comp*1/3*(divNE))\n", @@ -303,7 +325,7 @@ "TE = Dv*EdE\n", "TxzE = TE[2].subs({D31:'D31E', D32:'D32E', D33:'D33E', D34:'D34E'})\n", "#----------------------------------------#\n", - "Fz = 1/dz*(inN*TzzN - inS*TzzS) + 1/dx*(inE*TxzE - inW*TxzW) - (inN*pN - inS*pS)/dz \n", + "Fz = 1/dz*(TzzN - TzzS) + 1/dx*(TxzE - TxzW) - (pN - pS)/dz \n", "Fz = -Fz\n", "for i in range(len(dofs)):\n", " cVc = Fz.diff(dofs[i])\n",