Thanks for the update.
I came up with another approach, but you will only see a slightly performance benefit for larger sizes.
Your code seems to be alright for this workload:
def fun1(a, b):
    c = torch.cat((a, b), dim=0)
    idx = c.nonzero()
    res = torch.zeros(c.size()[1:])
    res[idx[:, 1:].split(1, dim=1)] = c[idx.split(1, dim=1)]
    return res

def fun2(a, b):
    d = (a!=0).int()
    e = (b!=0).int()
    f = d + e
    final = (a + b)/f
    res = torch.where(final == final, final, torch.tensor(0.)) # to remove nan values
    return res

a = torch.zeros(1, 1, 100, 100)
a[:, :, torch.randperm(100)[:10], torch.randperm(100)[:10]] = torch.randn(1)
b = torch.zeros(1, 1, 100, 100)
b[:, :, torch.randperm(100)[:10], torch.randperm(100)[:10]] = torch.randn(1)

res1 = fun1(a, b)
res2 = fun2(a, b)
print((res1 == res2).all())
> tensor(True)

%timeit fun1(a, b)
%timeit fun2(a, b)