diff --git a/Solver/TESTCASES/Diffusion_supra.lua b/Solver/TESTCASES/Diffusion_supra.lua new file mode 100644 index 0000000000000000000000000000000000000000..189d2439be3fc0865048744a355518b1146dd5c8 --- /dev/null +++ b/Solver/TESTCASES/Diffusion_supra.lua @@ -0,0 +1,77 @@ +model = GModel () +model:load ('plaque_supra.geo') +model:load ('plaque_supra.msh') +dg = dgSystemOfEquations (model) +dg:setOrder(1) + +Ec=1e-7 +Jc=100 +mu0=4*math.pi*1e-10 +c=mu0*Jc/Ec + +Bmax = 1e-6/Ec +t=0 +dt=1e-6 + + +function diffusivity(f ) + for i=0,f:size1()-1 do + f:set (i, 0, 1/c ) + end +end + +function valueRight(f) + for i=0,f:size1()-1 do + f:set(i,0,Bmax) + end +end + +function valueLeft(f) + for i=0,f:size1()-1 do + f:set(i,0,Bmax) + end +end + + + +nu=fullMatrix(1,1); +nu:set(0,0,1/c) + +-- conservation law +-- advection speed +law = dgConservationLawAdvectionDiffusion.diffusionLaw(functionLua(1,'diffusivity',{})) +dg:setConservationLaw(law) + +law:addBoundaryCondition('Top',law:new0FluxBoundary()) +law:addBoundaryCondition('Bottom',law:new0FluxBoundary()) +leftFunction=functionLua(1,'valueLeft',{}) +rightFunction=functionLua(1,'valueRight',{}) + + +law:addBoundaryCondition('Right',law:newNeumannBoundary(rightFunction)) +law:addBoundaryCondition('Left',law:newNeumannBoundary(leftFunction)) +dg:setup() + +-- initial condition +function initial_condition( xyz , f ) + for i=0,xyz:size1()-1 do + x = xyz:get(i,0) + y = xyz:get(i,1) + z = xyz:get(i,2) + f:set (i, 0, 0*math.exp(-100*((x-0.2)^2 +(y-0.3)^2))) + end +end + +dg:L2Projection(functionLua(1,'initial_condition',{functionCoordinates.get()})) +dg:exportSolution('output/DiffusionSupra_00000') + +-- main loop +for i=1,1000 do + -- norm = dg:RK44(0.01) + norm = dg:RK44_TransformNodalValue_limiter(0.00001) + if (i % 50 == 0) then + print('iter',i,norm) + dg:exportSolution(string.format("output/DiffusionSupra-%05d", i)) + end +t=t+dt +end diff --git a/Solver/TESTCASES/plaque_supra.geo b/Solver/TESTCASES/plaque_supra.geo new file mode 100644 index 0000000000000000000000000000000000000000..a4e076ba1cc308556a7c87fd47be0308ea5bc0c2 --- /dev/null +++ b/Solver/TESTCASES/plaque_supra.geo @@ -0,0 +1,23 @@ +lc=0.25; +a=5; +b=0.5; + +Point(1) = {-a, -b, 0, lc}; +Point(2) = {-a, b, 0, lc}; +Point(3) = {a, b, 0, lc}; +Point(4) = {a, -b, 0, lc}; +Line(1) = {4, 3}; +Line(2) = {3, 2}; +Line(3) = {2, 1}; +Line(4) = {1, 4}; +Line Loop(5) = {2, 3, 4, 1}; +Plane Surface(6) = {5}; + +Physical Line("Right")={1}; +Physical Line("Top")={2}; +Physical Line("Left")={3}; +Physical Line("Bottom")={4}; +Physical Surface("Inside") = {6}; +//Mesh.CharacteristicLengthExtendFromBoundary=1; +//Transfinite Line {1, 2, 4, 3} = 2 Using Progression 1; +//Transfinite Surface {6};